From 3b536d43328a6dc05f7da5f3e9b961fc5c1b19cf Mon Sep 17 00:00:00 2001 From: pruliere Date: Mon, 6 Jul 2026 09:36:33 +0200 Subject: [PATCH 01/14] Add MultiMesh postprocessing and FDH5 dataset support --- .../01-simple/spherical_shell_compression.py | 21 +- fedoo/__init__.py | 1 + fedoo/core/__init__.py | 2 + fedoo/core/base.py | 4 +- fedoo/core/dataset.py | 740 +++++++++++- fedoo/core/mesh.py | 511 ++++++++- fedoo/core/output.py | 228 +++- fedoo/core/problem.py | 6 +- fedoo/util/fdh5.py | 1007 +++++++++++++++++ fedoo/util/viewer.py | 210 +++- fedoo/util/xdmf_writer.py | 267 +++++ 11 files changed, 2845 insertions(+), 152 deletions(-) create mode 100644 fedoo/util/fdh5.py create mode 100644 fedoo/util/xdmf_writer.py diff --git a/examples/01-simple/spherical_shell_compression.py b/examples/01-simple/spherical_shell_compression.py index 125b9950..b1faad20 100644 --- a/examples/01-simple/spherical_shell_compression.py +++ b/examples/01-simple/spherical_shell_compression.py @@ -62,12 +62,25 @@ # Define a linear analysis and solve the problem. # # .. note:: -# Here we don't need to add other boundary conditions. The rigid body -# displacements and rotations of the sphere aren't constrained but the solver -# find a solution that is unique in terms of strain and stress (but not -# for displacements or rotations) +# To improve numerical stability, a few displacement boundary conditions are +# added to remove rigid-body motions. These constraints do not affect the +# strain/stress solution because they only suppress the nullspace modes. +# Without them, the unconstrained problem is singular; some solvers may still +# return a usable strain/stress field, but the displacement field can contain +# arbitrary rigid-body motion. +# pb = fd.problem.Linear(assembly) + +nodes = mesh.nodes +node_a = int(np.argmin(nodes[:, 0])) +node_b = int(np.argmax(nodes[:, 0])) +node_c = int(np.argmax(nodes[:, 1])) + +pb.bc.add("Dirichlet", node_a, "Disp", 0) +pb.bc.add("Dirichlet", node_b, ["DispY", "DispZ"], 0) +pb.bc.add("Dirichlet", node_c, "DispZ", 0) + pb.solve() ############################################################################### diff --git a/fedoo/__init__.py b/fedoo/__init__.py index 25017e6e..1177adcb 100644 --- a/fedoo/__init__.py +++ b/fedoo/__init__.py @@ -20,6 +20,7 @@ Mesh, ModelingSpace, MultiFrameDataSet, + MultiMeshData, Problem, WeakForm, WeakFormSum, diff --git a/fedoo/core/__init__.py b/fedoo/core/__init__.py index 662d1356..5223bd78 100644 --- a/fedoo/core/__init__.py +++ b/fedoo/core/__init__.py @@ -9,6 +9,7 @@ MultiFrameDataSet, read_data, ) +from .multimeshdata import MultiMeshData from .boundary_conditions import BoundaryCondition, MPC, ListBC from .problem import Problem @@ -23,6 +24,7 @@ "ModelingSpace", "DataSet", "MultiFrameDataSet", + "MultiMeshData", "read_data", "BoundaryCondition", "MPC", diff --git a/fedoo/core/base.py b/fedoo/core/base.py index 523c24e7..13fcb5c0 100644 --- a/fedoo/core/base.py +++ b/fedoo/core/base.py @@ -88,11 +88,11 @@ class MeshBase: __dic: dict[str, "MeshBase"] = {} - def __init__(self, name=""): + def __init__(self, name="", register=True): assert isinstance(name, str), "name must be a string" self.__name = name - if name != "": + if register and name != "": MeshBase.__dic[self.__name] = self def __class_getitem__(cls, item): diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 44eeed7d..46fa6569 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -4,8 +4,9 @@ import numpy as np import os -from zipfile import ZipFile, Path -from fedoo.core.mesh import Mesh +from zipfile import ZipFile, Path as ZipPath +from fedoo.core.mesh import Mesh, MultiMesh +from fedoo.core.multimeshdata import MultiMeshData, copy_data_value from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList try: @@ -36,6 +37,59 @@ USE_PANDAS = False +def _as_3d_points(points: np.ndarray) -> np.ndarray: + """Return point coordinates padded or truncated to 3D.""" + if points.shape[1] == 3: + return points + if points.shape[1] < 3: + return np.column_stack( + (points, np.zeros((points.shape[0], 3 - points.shape[1]))) + ) + return points[:, :3] + + +def _array_to_pyvista_data(data: np.ndarray) -> np.ndarray: + """Return data in the component-last shape expected by pyvista.""" + data = np.asarray(data) + if data.ndim == 1: + return data + return data.T + + +def _component_from_array(field: str, data: np.ndarray, component=None): + """Extract one component from a Fedoo data array.""" + if ( + component is None + or np.isscalar(data) + or not hasattr(data, "shape") + or len(data.shape) <= 1 + ): + return data + + if isinstance(component, str): + component = { + "X": 0, + "Y": 1, + "Z": 2, + "XX": 0, + "YY": 1, + "ZZ": 2, + "XY": 3, + "XZ": 4, + "YZ": 5, + }.get(component, component) + + if component == "norm": + return np.linalg.norm(data, axis=0) + + if isinstance(component, str): + if field == "Stress": + data = StressTensorList(data) + elif field == "Strain": + data = StrainTensorList(data) + return data[component] + + class DataSet: """ Object to store, save, load and plot data associated to a mesh. @@ -95,6 +149,14 @@ def __init__( self.element_data = {} self.gausspoint_data = {} self.scalar_data = {} + self.active_submesh = 0 + """Active submesh used by MultiMesh data access. + + For ``MultiMesh`` datasets, element and Gauss point fields are returned + as ``MultiMeshData`` objects. Plain array access to these objects points + to ``active_submesh``. Node fields stay plain arrays; by default, + ``get_data`` fills nodes unused by the active submesh with zero. + """ if isinstance(data, dict): data_type = data_type.lower() @@ -115,6 +177,114 @@ def __init__( self.meshplot = None self.meshplot_gp = None # a mesh with discontinuity between each element to plot gauss points field + def _is_multimesh(self) -> bool: + return isinstance(self.mesh, MultiMesh) + + def _resolve_submesh_indices(self, selector=None) -> list[int]: + if not self._is_multimesh(): + return [0] + if selector is None: + return list(range(len(self.mesh.submeshes))) + if isinstance(selector, (list, tuple, set, np.ndarray)): + indices = [] + for item in selector: + indices.extend(self._resolve_submesh_indices(item)) + return list(dict.fromkeys(indices)) + if isinstance(selector, int): + return [selector] + if isinstance(selector, str): + name_matches = [ + i for i, mesh in enumerate(self.mesh.submeshes) if mesh.name == selector + ] + if name_matches: + return name_matches + type_matches = [ + i + for i, mesh in enumerate(self.mesh.submeshes) + if mesh.elm_type == selector + ] + if type_matches: + return type_matches + raise KeyError(selector) + + def _active_submesh_index(self) -> int: + return self._resolve_submesh_indices(self.active_submesh)[0] + + def _selected_multimesh(self, selected_submeshes=None): + if not self._is_multimesh(): + return self.mesh, [0] + indices = self._resolve_submesh_indices(selected_submeshes) + mesh = MultiMesh.from_mesh_list( + [self.mesh[i] for i in indices], + name=self.mesh.name, + node_sets=self.mesh.node_sets, + register_name=False, + ) + return mesh, indices + + def _as_multimesh_data(self, data) -> MultiMeshData: + return MultiMeshData( + self.mesh, + data, + active_submesh=self._active_submesh_index(), + ) + + def _convert_multimesh_data(self, data, convert_from: str, convert_to: str): + """Convert MultiMesh data block-by-block on each submesh.""" + active_submesh = self._active_submesh_index() + + if convert_from == convert_to: + if convert_to in ["Element", "GaussPoint"]: + return self._as_multimesh_data(data) + return data + + if isinstance(data, MultiMeshData): + if convert_to == "Node": + block = data.submesh(active_submesh) + if block is None: + raise NameError("Field data not found on the active submesh.") + return self.mesh[active_submesh].convert_data( + np.asarray(block), + convert_from=convert_from, + convert_to=convert_to, + ) + + converted = { + submesh_id: self.mesh[submesh_id].convert_data( + np.asarray(block), + convert_from=convert_from, + convert_to=convert_to, + ) + for submesh_id, block in data.items() + if block is not None + } + return MultiMeshData( + self.mesh, + converted, + active_submesh=active_submesh, + ) + + if convert_from == "Node" and convert_to in ["Element", "GaussPoint"]: + converted = { + submesh_id: submesh.convert_data( + data, + convert_from=convert_from, + convert_to=convert_to, + ) + for submesh_id, submesh in enumerate(self.mesh.submeshes) + } + return MultiMeshData( + self.mesh, + converted, + active_submesh=active_submesh, + ) + + return self.mesh[active_submesh].convert_data( + data, + convert_from=convert_from, + convert_to=convert_to, + ) + def __getitem__(self, items): if isinstance(items, tuple): return self.get_data(*items) @@ -165,6 +335,7 @@ def plot( multiplot: bool | None = None, element_set: str | np.ndarray[int] | None = None, element_set_invert: bool = False, + selected_submeshes=None, clip_args: tuple | None = None, lock_view: bool = False, iteration: int | None = None, @@ -290,6 +461,10 @@ def plot( element_set_invert : bool, optional Used only if element_set is defined. Invert element set. + selected_submeshes : int, str or list, optional + Only used with MultiMesh objects. Restrict the plot to the selected + submesh ids, submesh names or element types. + clip_args : dict, optional Dictionary of arguments to pass to the pyvista clip filter in order to clip the current plot. @@ -341,6 +516,36 @@ def plot( "scalars", field ) # kargs scalars can be used instead of field + if self._is_multimesh(): + return self._plot_multimesh( + field=field, + component=component, + data_type=data_type, + scale=scale, + show=show, + show_edges=show_edges, + clim=clim, + node_labels=node_labels, + element_labels=element_labels, + show_nodes=show_nodes, + plotter=plotter, + screenshot=screenshot, + azimuth=azimuth, + elevation=elevation, + roll=roll, + title=title, + title_size=title_size, + window_size=window_size, + multiplot=multiplot, + element_set=element_set, + element_set_invert=element_set_invert, + selected_submeshes=selected_submeshes, + clip_args=clip_args, + lock_view=lock_view, + return_cpos=kargs.pop("return_cpos", False), + **kargs, + ) + if field is not None: data, data_type = self.get_data(field, component, data_type, True) else: @@ -351,6 +556,7 @@ def plot( return_cpos = kargs.pop("return_cpos", False) cmap = kargs.pop("cmap", "jet") # if cmap not defined, default to "jet" + extra_cell_data = kargs.pop("_extra_cell_data", None) if data_type == "GaussPoint": if self.meshplot_gp is None: @@ -404,6 +610,10 @@ def plot( center = 0.5 * (meshplot.points.min(axis=0) + meshplot.points.max(axis=0)) + if extra_cell_data: + for key, value in extra_cell_data.items(): + meshplot.cell_data[key] = np.asarray(value) + backgroundplotter = True if USE_PYVISTA_QT and (plotter is None or plotter == "qt"): # use pyvistaqt plotter @@ -647,7 +857,307 @@ def plot( return pl - def get_data(self, field, component=None, data_type=None, return_data_type=False): + def _multimesh_block_for_submesh( + self, + data, + submesh_id: int, + n_items: int, + fill_missing: bool = True, + ): + """Return one submesh data block without copying when possible.""" + multimesh_data = self._as_multimesh_data(data) + block = multimesh_data.submesh(submesh_id) + if block is not None or not fill_missing: + return block + + template = next( + (value for _, value in multimesh_data.items() if value is not None), + None, + ) + if template is None: + return None + + template = np.asarray(template) + if template.ndim > 1: + shape = template.shape[:-1] + (n_items,) + else: + shape = (n_items,) + return np.full(shape, np.nan) + + def _submesh_dataset(self, submesh_id: int) -> "DataSet": + """Build a lightweight DataSet view attached to one submesh.""" + submesh = self.mesh[submesh_id] + dataset = DataSet(submesh) + dataset.node_data = self.node_data + dataset.scalar_data = self.scalar_data + dataset.element_data = { + field: block + for field, value in self.element_data.items() + if ( + block := self._multimesh_block_for_submesh( + value, + submesh_id, + submesh.n_elements, + ) + ) + is not None + } + dataset.gausspoint_data = { + field: block + for field, value in self.gausspoint_data.items() + if ( + block := self._multimesh_block_for_submesh( + value, + submesh_id, + submesh.n_elements, + fill_missing=False, + ) + ) + is not None + } + return dataset + + def _submesh_element_set( + self, + element_set, + submesh_id: int, + offsets: dict[int, int], + global_element_set: bool, + element_set_invert: bool, + ): + """Map a MultiMesh element selection to one submesh.""" + if element_set is None: + return None, True + + submesh = self.mesh[submesh_id] + if isinstance(element_set, str): + if element_set in submesh.element_sets: + return element_set, True + return None, bool(element_set_invert) + + element_ids = np.asarray(element_set, dtype=int) + if not global_element_set: + return element_ids, True + + offset = offsets[submesh_id] + local_ids = element_ids[ + (element_ids >= offset) & (element_ids < offset + submesh.n_elements) + ] - offset + if len(local_ids) == 0 and not element_set_invert: + return None, False + if len(local_ids) == 0 and element_set_invert: + return None, True + return local_ids, True + + def _plot_multimesh( + self, + *, + field, + component, + data_type, + scale, + show, + show_edges, + clim, + node_labels, + element_labels, + show_nodes, + plotter, + screenshot, + azimuth, + elevation, + roll, + title, + title_size, + window_size, + multiplot, + element_set, + element_set_invert, + selected_submeshes, + clip_args, + lock_view, + return_cpos, + **kargs, + ): + """Plot a MultiMesh by reusing the single-mesh plotting path.""" + global_element_set = kargs.pop("global_element_set", False) + submesh_indices = self._resolve_submesh_indices(selected_submeshes) + if ( + element_set is not None + and not isinstance(element_set, str) + and not global_element_set + ): + submesh_indices = [self._active_submesh_index()] + + resolved_data_type = None + if field is not None: + _, resolved_data_type = self.get_data( + field, + component, + data_type, + True, + fill_unused_nodes=None, + ) + + offsets = {} + offset = 0 + for i, submesh in enumerate(self.mesh.submeshes): + offsets[i] = offset + offset += submesh.n_elements + + if screenshot and (plotter is None or plotter == "pv"): + pl = pv.Plotter(off_screen=True, window_size=window_size) + backgroundplotter = False + else: + pl = plotter + backgroundplotter = not (plotter is None or plotter == "pv") + if USE_PYVISTA_QT and (plotter is None or plotter == "qt"): + backgroundplotter = True + + base_name = kargs.pop("name", None) + scalar_bar_added = False + plotted = False + + for submesh_id in submesh_indices: + local_element_set, should_plot = self._submesh_element_set( + element_set, + submesh_id, + offsets, + global_element_set, + element_set_invert, + ) + if not should_plot: + continue + + subdataset = self._submesh_dataset(submesh_id) + if ( + field is not None + and resolved_data_type == "GaussPoint" + and field not in subdataset.gausspoint_data + ): + continue + sub_kargs = dict(kargs) + if base_name is not None: + sub_kargs["name"] = f"{base_name}_{submesh_id}" + if field is not None and scalar_bar_added: + sub_kargs["show_scalar_bar"] = False + sub_kargs["_extra_cell_data"] = { + "_fedoo_global_cell_ids": ( + np.arange(subdataset.mesh.n_elements, dtype=int) + + offsets[submesh_id] + ), + "_fedoo_submesh_id": np.full( + subdataset.mesh.n_elements, + submesh_id, + dtype=int, + ), + } + if clip_args is not None: + sub_clip_args = dict(clip_args) + else: + sub_clip_args = None + + pl = subdataset.plot( + field=field, + component=component, + data_type=data_type, + scale=scale, + show=False, + show_edges=show_edges, + clim=clim, + node_labels=False, + element_labels=element_labels, + show_nodes=False, + plotter=pl, + screenshot=None, + azimuth=azimuth, + elevation=elevation, + roll=roll, + title="", + title_size=title_size, + window_size=window_size, + multiplot=multiplot, + element_set=local_element_set, + element_set_invert=element_set_invert, + clip_args=sub_clip_args, + lock_view=True, + return_cpos=return_cpos, + **sub_kargs, + ) + plotted = True + if field is not None: + scalar_bar_added = True + + if pl is None: + pl = pv.Plotter(window_size=window_size) + backgroundplotter = False + + if title is None: + title = "" if field is None else f"{field}_{component}" + pl.add_text(title, name="name", color="Black", font_size=title_size) + + if not lock_view: + pl.set_background("White") + points = self.mesh.nodes + if "Disp" in self.node_data and scale != 0: + points = points + scale * np.asarray(self.node_data["Disp"]).T + points = _as_3d_points(points) + center = 0.5 * (points.min(axis=0) + points.max(axis=0)) + length = np.linalg.norm(points.max(axis=0) - points.min(axis=0)) + if length == 0: + length = 1 + pl.camera.SetFocalPoint(center) + pl.camera.position = tuple(center + np.array([0, 0, 2 * length])) + pl.camera.up = tuple([0, 1, 0]) + if roll != 0: + pl.camera.Roll(roll) + if self.mesh.ndim == 3: + pl.camera.Azimuth(azimuth) + pl.camera.Elevation(elevation) + pl.add_axes(color="Black", interactive=True) + + if plotted and (node_labels or show_nodes): + crd_labels = self.mesh.nodes + if "Disp" in self.node_data and scale != 0: + crd_labels = crd_labels + scale * np.asarray(self.node_data["Disp"]).T + crd_labels = _as_3d_points(crd_labels) + if node_labels: + if node_labels is True: + node_labels = list(range(self.mesh.n_nodes)) + pl.add_point_labels(crd_labels, node_labels) + if show_nodes: + if show_nodes is True: + show_nodes = 5 + pl.add_points( + crd_labels, + render_points_as_spheres=True, + point_size=show_nodes, + ) + + pl.renderer.ResetCameraClippingRange() + + if screenshot: + ext = os.path.splitext(screenshot)[1].lower() + if ext in [".pdf", ".svg", ".eps", ".ps", ".tex"]: + pl.save_graphic(screenshot) + else: + pl.screenshot(screenshot) + return pl + + if not backgroundplotter and show: + return pl.show(return_cpos=return_cpos) + + return pl + + def get_data( + self, + field, + component=None, + data_type=None, + return_data_type=False, + *, + fill_unused_nodes=0.0, + ): """Retrieve data from the DataSet for a given field. This method is equivalent to the `DataSet.__getitem__` magic method. @@ -670,11 +1180,17 @@ def get_data(self, field, component=None, data_type=None, return_data_type=False If True, the method returns a tuple `(data, data_type)` where `data_type` is the type of the returned data. If False, only the data is returned. + fill_unused_nodes : scalar, optional + Value used for nodes that are not used by the active submesh when a + MultiMesh node field is restricted to one submesh. Default is 0. + If ``None``, the original node array is returned unchanged. Returns ------- - data : np.ndarray - The requested data, possibly converted to the specified type. + data : np.ndarray or MultiMeshData + The requested data, possibly converted to the specified type. For a + ``MultiMesh`` dataset, element and Gauss point fields are returned + as ``MultiMeshData`` objects. Node fields remain NumPy arrays. data_type : str, optional Returned only if `return_data_type` is True. Indicates the type of the returned data. @@ -683,6 +1199,10 @@ def get_data(self, field, component=None, data_type=None, return_data_type=False ----- This method supports automatic conversion between node, element, and Gauss point data types when applicable. + + With a ``MultiMesh``, element and Gauss point fields may be stored as a + dictionary keyed by submesh id, submesh name, or unique element type. + The active submesh is controlled by ``dataset.active_submesh``. """ if data_type is None: # search if field exist somewhere @@ -702,38 +1222,45 @@ def get_data(self, field, component=None, data_type=None, return_data_type=False data = self.dict_data[data_type][field] else: # if field is not present whith the given data_type search if it exist elsewhere and convert it data, current_data_type = self.get_data( - field, component, return_data_type=True + field, + component, + return_data_type=True, + fill_unused_nodes=fill_unused_nodes, ) if current_data_type != "Scalar": - data = self.mesh.convert_data( - data, convert_from=current_data_type, convert_to=data_type - ) + if self._is_multimesh(): + data = self._convert_multimesh_data( + data, + convert_from=current_data_type, + convert_to=data_type, + ) + else: + data = self.mesh.convert_data( + data, + convert_from=current_data_type, + convert_to=data_type, + ) + + if self._is_multimesh() and data_type in ["Element", "GaussPoint"]: + data = self._as_multimesh_data(data) + + if isinstance(data, MultiMeshData): + data = data.map(lambda val: _component_from_array(field, val, component)) + else: + data = _component_from_array(field, data, component) if ( - component is not None and not (np.isscalar(data)) and len(data.shape) > 1 - ): # if data is scalar or 1d array, component ignored - if isinstance(component, str): - component = { - "X": 0, - "Y": 1, - "Z": 2, - "XX": 0, - "YY": 1, - "ZZ": 2, - "XY": 3, - "XZ": 4, - "YZ": 5, - }.get(component, component) - - if component == "norm": - data = np.linalg.norm(data, axis=0) - else: - if isinstance(component, str): - if field == "Stress": - data = StressTensorList(data) - elif field == "Strain": - data = StrainTensorList(data) - data = data[component] + self._is_multimesh() + and data_type == "Node" + and fill_unused_nodes is not None + and field in self.node_data + ): + active_mesh = self.mesh[self._active_submesh_index()] + used_nodes = np.unique(active_mesh.elements) + unused_nodes = np.setdiff1d(np.arange(self.mesh.n_nodes), used_nodes) + if len(unused_nodes): + data = np.array(data, copy=True) + data[..., unused_nodes] = fill_unused_nodes if return_data_type: return data, data_type @@ -745,7 +1272,7 @@ def field_names(self): set( list(self.gausspoint_data.keys()) + list(self.node_data.keys()) - + list(self.node_data.keys()) + + list(self.element_data.keys()) ) ) @@ -774,6 +1301,8 @@ def save( pandas installed). The mesh is not included and may be saved beside in a vtk file. * 'xlsx': Same as csv but with the excel format. + * 'fdh5': HDF5 format for Fedoo meshes and results, including + MultiMesh element and Gauss point fields. Parameters ---------- @@ -788,7 +1317,7 @@ def save( ext = os.path.splitext(filename)[1] ext = ext.lower() if ext == "": - ext = ".fdz" + ext = ".fdh5" filename = filename + ext if ext == ".vtk": self.to_vtk(filename) @@ -807,6 +1336,8 @@ def save( self.to_fdz( filename, save_mesh=True, compressed=compressed ) # create a new file and add the mesh + elif ext == ".fdh5": + self.to_fdh5(filename, iteration=0, overwrite=True) def save_mesh(self, filename: str): """Save the mesh using a vtk file. The extension of filename is ignored and modified to '.vtk'.""" @@ -833,7 +1364,7 @@ def load(self, data: object, load_mesh: bool = False, iteration: int = 0): copy. * **str** : Path to a data file. Supported file extensions are - ``'vtk'``, ``'msh'``, ``'fdz'``, and ``'npz'``. + ``'vtk'``, ``'msh'``, ``'fdz'``, ``'fdh5'`` and ``'npz'``. load_mesh : bool, optional If ``True``, the mesh is loaded from the file (when the file @@ -858,7 +1389,7 @@ def load(self, data: object, load_mesh: bool = False, iteration: int = 0): self.element_data = {k: v.T for k, v in data.cell_data.items()} if load_mesh: self.mesh = Mesh.from_pyvista(data) - elif isinstance(data, Path): + elif isinstance(data, ZipPath): # used to load one iteration in fdz file data = np.load(data.open("rb")) self.load_dict(data) @@ -876,6 +1407,10 @@ def load(self, data: object, load_mesh: bool = False, iteration: int = 0): DataSet.load(self, pv.read(filename)) elif ext == ".msh": return NotImplemented + elif ext == ".fdh5": + from fedoo.util.fdh5 import load_dataset_iteration + + load_dataset_iteration(self, filename, iteration) elif ext in [".npz", ".fdz"]: if ext == ".fdz": file = ZipFile(filename, "r") @@ -1030,20 +1565,45 @@ def to_vtk( write_vtk(self, filename, gp_data_to_node) - def to_pyvista(self, gp_data_to_node: bool = True): + def to_pyvista(self, gp_data_to_node: bool = True, selected_submeshes=None): + """Convert the dataset to a PyVista unstructured grid. + + Parameters + ---------- + gp_data_to_node : bool, default=True + For a single ``Mesh``, convert Gauss point data to node data before + export. For a ``MultiMesh``, Gauss point conversion is currently + skipped. + selected_submeshes : int, str or sequence, optional + Only used with ``MultiMesh`` datasets. Restrict the exported grid + and element data to the selected submesh ids, names, or element + types. + + Returns + ------- + pyvista.UnstructuredGrid + Mesh with point, cell and field data attached. + """ if self.mesh is not None: - pv_data = self.mesh.to_pyvista() + if self._is_multimesh(): + mesh, submesh_indices = self._selected_multimesh(selected_submeshes) + else: + mesh, submesh_indices = self.mesh, None + + pv_data = mesh.to_pyvista() for key, val in self.node_data.items(): pv_data.point_data[key] = val.T for key, val in self.element_data.items(): - pv_data.cell_data[key] = val.T + if self._is_multimesh(): + val = self._as_multimesh_data(val).to_global(submesh_indices) + pv_data.cell_data[key] = _array_to_pyvista_data(val) for key, val in self.scalar_data.items(): pv_data.field_data[key] = np.array(val) - if gp_data_to_node: + if gp_data_to_node and not self._is_multimesh(): for key in self.gausspoint_data: pv_data.point_data[key] = self.get_data(key, data_type="Node").T @@ -1108,6 +1668,35 @@ def to_fdz( os.remove("_mesh_.vtk") file.close() + def to_fdh5( + self, + filename: str, + iteration: int = 0, + overwrite: bool = False, + ) -> None: + """Write the dataset to a FDH5 file. + + Parameters + ---------- + filename : str + Name of the FDH5 file. If no extension is provided, ``.fdh5`` is + appended. + iteration : int, default=0 + Result iteration id written under ``results/iter_``. + overwrite : bool, default=False + If True, an existing file is removed before writing. If False, the + mesh is kept and only the requested iteration is added or replaced. + + Notes + ----- + Element and Gauss point fields attached to a ``MultiMesh`` are written + under their matching ``submesh_X`` groups. Single ``Mesh`` datasets are + written under ``submesh_0``. + """ + from fedoo.util.fdh5 import write_dataset + + write_dataset(self, filename, iteration=iteration, overwrite=overwrite) + def savez(self, filename: str, save_mesh: bool = False) -> None: """Write a npz file using the numpy savez function. @@ -1139,12 +1728,12 @@ def savez_compressed(self, filename: str, save_mesh: bool = False) -> None: self.save_mesh(filename) @staticmethod - def read(filename: str, file_format: str = "fdz") -> DataSet | MultiFrameDataSet: + def read(filename: str, file_format: str = "fdh5") -> DataSet | MultiFrameDataSet: """Read a file from disk. Same as :py:func:`fedoo.read_data`. """ - return read_data(filename, file_format="fdz") + return read_data(filename, file_format=file_format) @property def dict_data(self) -> dict: @@ -1172,6 +1761,7 @@ def copy(self): copy.element_data = dict(self.element_data) copy.gausspoint_data = dict(self.gausspoint_data) copy.scalar_data = dict(self.scalar_data) + copy.active_submesh = self.active_submesh return copy def deepcopy(self): @@ -1183,17 +1773,23 @@ def deepcopy(self): """ copy = DataSet() copy.mesh = self.mesh.deepcopy() - copy.node_data = {key: value.copy() for key, value in self.node_data.items()} + copy.node_data = { + key: copy_data_value(value, deep=True) + for key, value in self.node_data.items() + } copy.element_data = { - key: value.copy() for key, value in self.element_data.items() + key: copy_data_value(value, deep=True) + for key, value in self.element_data.items() } copy.gausspoint_data = { - key: value.copy() for key, value in self.gausspoint_data.items() + key: copy_data_value(value, deep=True) + for key, value in self.gausspoint_data.items() } copy.scalar_data = { key: value if np.isscalar(value) else np.array(value).copy() for key, value in self.scalar_data.items() } + copy.active_submesh = self.active_submesh return copy @@ -1215,20 +1811,21 @@ def __getitem__(self, items): return DataSet.__getitem__(self, items) def save_all( - self, filename: str, file_format: str = "fdz", compressed: bool = False + self, filename: str, file_format: str = "fdh5", compressed: bool = False ): """Save all data from MultiFrameDataSet. If filename has no extension, the format is given in the parameter file_format - (default = 'fdz'). - If format is not 'fdz', the data files are saved using the given filename and format + (default = 'fdh5'). + If format is not 'fdz' or 'fdh5', the data files are saved using the given filename and format simply adding the iteration number to the file name. The mesh is also saved in vtk format in the same directory. + For 'fdh5', all iterations are written into one HDF5 file. """ dirname = os.path.dirname(filename) extension = os.path.splitext(filename)[1] if extension == "": file_format = file_format.lower() - if file_format != "fdz": + if file_format not in ["fdz", "fdh5"]: dirname = filename + "/" filename = dirname + os.path.basename(filename) else: @@ -1246,6 +1843,12 @@ def save_all( for i in range(1, len(self.list_data)): self.load(i) self.to_fdz(filename, False, i, compressed) + elif file_format == "fdh5": + if os.path.splitext(filename)[1] == "": + filename += ".fdh5" + for i in range(len(self.list_data)): + self.load(i) + self.to_fdh5(filename, iteration=i, overwrite=(i == 0)) else: for i in range(len(self.list_data)): self.load(i) @@ -1267,9 +1870,17 @@ def load(self, data=-1, load_mesh=False): return if iteration > len(self.list_data) or iteration < 0: raise NameError("Number of iteration out of bounds") - DataSet.load(self, self.list_data[iteration]) + data_ref = self.list_data[iteration] + if isinstance(data_ref, tuple) and len(data_ref) == 3 and data_ref[0] == "fdh5": + DataSet.load(self, data_ref[1], load_mesh, iteration=data_ref[2]) + else: + DataSet.load(self, data_ref) self.loaded_iter = iteration + elif isinstance(data, tuple) and len(data) == 3 and data[0] == "fdh5": + DataSet.load(self, data[1], load_mesh, iteration=data[2]) + self.loaded_iter = data[2] + elif data: DataSet.load(self, data, load_mesh) @@ -1506,6 +2117,8 @@ def get_history( self.load(it) for i, field in enumerate(list_fields): data = self.get_data(field, component[i], data_type[i]) + if isinstance(data, MultiMeshData): + data = data.to_global() if list_indices[i] is None or np.isscalar( data ): # modify for allowing scalar_data @@ -1604,9 +2217,11 @@ def get_all_frame_lim(self, field, component=0, data_type=None, scale=1): self.load(i) if field is not None: data = self.get_data(field, component, data_type) + if isinstance(data, MultiMeshData): + data = data.to_global() clim = [ - np.min([data.min(), clim[0]]), - np.max([data.max(), clim[1]]), + np.nanmin([np.nanmin(data), clim[0]]), + np.nanmax([np.nanmax(data), clim[1]]), ] if "Disp" in self.node_data: @@ -1641,6 +2256,7 @@ def copy(self): The copied MultiFrameDataSet object. """ copy = MultiFrameDataSet(self.mesh.copy(), list(self.list_data)) + copy.active_submesh = self.active_submesh if self.loaded_iter: copy.load(self.loaded_iter) return copy @@ -1656,6 +2272,7 @@ def deepcopy(self): The copied MultiFrameDataSet object. """ copy = MultiFrameDataSet(self.mesh.deepcopy(), list(self.list_data)) + copy.active_submesh = self.active_submesh if self.loaded_iter: copy.load(self.loaded_iter) return copy @@ -1665,16 +2282,17 @@ def n_iter(self): return len(self.list_data) -def read_data(filename: str, file_format: str = "fdz"): +def read_data(filename: str, file_format: str = "fdh5"): """Read a file from disk. The file may be a single file or a directory containing files from several iterations. The file format may be specified either by the filename extension or by the ``file_format`` parameter (default is - ``"fdz"``) when the filename has no extension. + ``"fdh5"``) when the filename has no extension. - Supported file formats are ``"fdz"``, ``"vtk"``, and ``"npz"``. - For ``"npz"`` files, a VTK mesh with the same base name is also searched. + Supported file formats are ``"fdz"``, ``"fdh5"``, ``"vtk"``, and + ``"npz"``. For ``"npz"`` files, a VTK mesh with the same base name is + also searched. Parameters ---------- @@ -1682,7 +2300,7 @@ def read_data(filename: str, file_format: str = "fdz"): Path to the file or directory to read. file_format : str, optional File format identifier to use when the filename has no extension. - Default is ``"fdz"``. + Default is ``"fdh5"``. Returns ------- @@ -1696,6 +2314,10 @@ def read_data(filename: str, file_format: str = "fdz"): if file_format == "fdz": return read_fdz(filename) + if file_format == "fdh5": + from fedoo.util.fdh5 import read_fdh5 + + return read_fdh5(filename) dirname = os.path.dirname(filename) if extension == "": @@ -1752,7 +2374,7 @@ def read_fdz(filename: str): dataset = MultiFrameDataSet(mesh) i = 0 while "iter_" + str(i) + ".npz" in list_iter: - dataset.list_data.append(Path(filename, "iter_" + str(i) + ".npz")) + dataset.list_data.append(ZipPath(filename, "iter_" + str(i) + ".npz")) i += 1 return dataset diff --git a/fedoo/core/mesh.py b/fedoo/core/mesh.py index f434ee75..0227af94 100644 --- a/fedoo/core/mesh.py +++ b/fedoo/core/mesh.py @@ -11,6 +11,7 @@ import re from os.path import splitext +from pathlib import Path try: import pyvista as pv @@ -27,6 +28,36 @@ USE_PYVISTA_QT = False +_PYVISTA_CELL_TYPES = { + "spring": (3, 2), + "lin2": (3, 2), + "tri3": (5, 3), + "quad4": (9, 4), + "tet4": (10, 4), + "hex8": (12, 8), + "wed6": (13, 6), + "pyr5": (14, 5), + "lin3": (21, 3), + "tri6": (22, 6), + "quad8": (23, 8), + "tet10": (24, 10), + "hex20": (25, 20), + "wed15": (26, 15), + "pyr13": (27, 13), + "quad9": (28, 9), + "hex27": (29, 27), + "wed18": (32, 18), +} + + +def _get_pyvista_cell_info(elm_type: str) -> tuple[int, int]: + """Return the VTK cell type and number of element nodes for pyvista.""" + cell_info = _PYVISTA_CELL_TYPES.get(elm_type, None) + if cell_info is None: + raise NameError("Element Type " + str(elm_type) + " not available in pyvista") + return cell_info + + class Mesh(MeshBase): """Fedoo Mesh object. @@ -76,8 +107,9 @@ def __init__( element_sets: dict | None = None, ndim: int | None = None, name: str = "", + register_name: bool = True, ) -> None: - MeshBase.__init__(self, name) + MeshBase.__init__(self, name, register=register_name) self.nodes = nodes # node coordinates """List of nodes coordinates: nodes[i] gives the coordinates of the ith node.""" self.elements = elements # element table @@ -922,30 +954,7 @@ def to_pyvista(self, include_sets=False): pvmesh.field_data["ndim"] = self.ndim return pvmesh if USE_PYVISTA: - cell_type, n_elm_nodes = { - "spring": (3, 2), - "lin2": (3, 2), - "tri3": (5, 3), - "quad4": (9, 4), - "tet4": (10, 4), - "hex8": (12, 8), - "wed6": (13, 6), - "pyr5": (14, 5), - "lin3": (21, 3), - "tri6": (22, 6), - "quad8": (23, 8), - "tet10": (24, 10), - "hex20": (25, 20), - "wed15": (26, 15), - "pyr13": (27, 13), - "quad9": (28, 9), - "hex27": (29, 27), - "wed18": (32, 18), - }.get(self.elm_type, None) - if cell_type is None: - raise NameError( - "Element Type " + str(self.elm_type) + " not available in pyvista" - ) + cell_type, n_elm_nodes = _get_pyvista_cell_info(self.elm_type) # elm = np.empty((self.elements.shape[0], self.elements.shape[1]+1), dtype=int) elm = np.empty((self.elements.shape[0], n_elm_nodes + 1), dtype=int) @@ -973,9 +982,10 @@ def save( ) -> None: """Save the mesh object to file. - This function use the save function of the pyvista UnstructuredGrid - object. If using a 'vtk' format, The node and element sets are saved - as field data with a prefix "_nset_" or "_eset_". + This function uses the save function of the pyvista UnstructuredGrid + object, except for FDH5 files which are written directly as raw mesh + data. If using a 'vtk' format, The node and element sets are saved as + field data with a prefix "_nset_" or "_eset_". Parameters ---------- @@ -989,9 +999,32 @@ def save( extension = splitext(filename)[1] if extension == "": filename = filename + ".vtk" + extension = ".vtk" + + if extension.lower() == ".fdh5": + self._save_fdh5(filename) + return self.to_pyvista(include_sets=True).save(filename, binary=binary) + def _save_fdh5(self, filename: str) -> None: + """Save the mesh definition to a raw FDH5 file.""" + from fedoo.util.fdh5 import FDH5Writer + + path = Path(filename) + if path.exists(): + path.unlink() + + writer = FDH5Writer(path) + writer.write_mesh(self.nodes, node_sets=self.node_sets, overwrite=True) + writer.add_submesh( + self.elm_type, + self.elements, + element_sets=self.element_sets, + name=self.name or None, + submesh_id="submesh_0", + ) + def plot(self, show_edges: bool = True, **kargs) -> None: """Simple plot function using pyvista. @@ -1642,15 +1675,75 @@ def bounding_box(self) -> "BoundingBox": class MultiMesh(Mesh): + """Fedoo mesh made of several ordered submeshes sharing one node table. + + A ``MultiMesh`` is useful when a mesh contains several element blocks, + possibly with different element types or with several blocks using the + same element type. All submeshes share ``nodes``. Internally, submeshes + are stored in order; this order defines integer indexing and FDH5 + ``submesh_X`` ids. + + A ``MultiMesh`` can be created in three ways: + + * With a simple element-type dictionary, when each element type appears + once: + + >>> mesh = MultiMesh(nodes, {"tri3": tri_elements, "quad4": quad_elements}) + + * With an explicit submesh-name dictionary, when several submeshes may + share the same element type: + + >>> mesh = MultiMesh( + ... nodes, + ... { + ... "part_a": ("tri3", tri_elements_a), + ... "part_b": ("tri3", tri_elements_b), + ... }, + ... ) + + The dictionary key is stored as the submesh name, but internal + submeshes are not registered in ``MeshBase.get_all()``. + + * From an existing list of ``Mesh`` objects: + + >>> mesh = MultiMesh.from_mesh_list([tri_mesh, quad_mesh]) + + Submeshes can be retrieved by integer id, submesh name, or element type. + Integer indexing always returns a concrete ``Mesh``. String indexing first + tries exact submesh names, then element types. Element-type indexing + returns a ``Mesh`` if only one submesh matches, or a view-like + ``MultiMesh`` if several submeshes share that element type. + + Parameters + ---------- + nodes : numpy.ndarray + Shared node coordinates. + elements_dict : dict, optional + Submesh definitions. Values may be connectivity arrays, ``Mesh`` + objects, ``(elm_type, elements)`` tuples, or + ``(elm_type, elements, element_sets)`` tuples. + node_sets : dict, optional + Node sets associated with the shared node table. + ndim : int, optional + Dimension of the mesh. By default, deduced from ``nodes.shape[1]``. + name : str, optional + Name of the ``MultiMesh``. If non-empty and ``register_name`` is True, + the ``MultiMesh`` is registered in ``MeshBase.get_all()``. + register_name : bool, default=True + If False, the ``MultiMesh`` name is kept locally but not registered in + ``MeshBase.get_all()``. This is used for internal filtered views. + """ + def __init__( self, nodes: np.ndarray[float], - elements_dict: dict = None, + elements_dict: dict | None = None, node_sets: dict | None = None, ndim: int | None = None, name: str = "", + register_name: bool = True, ) -> None: - MeshBase.__init__(self, name) + MeshBase.__init__(self, name, register=register_name) self.nodes = nodes # node coordinates """ List of nodes coordinates: nodes[i] gives the coordinates of the ith node.""" @@ -1669,34 +1762,366 @@ def __init__( else: self.crd_name = ("X", "Y", "Z") - self.mesh_dict = {} + self._mesh_list = [] + """Ordered list of submeshes. The list order defines submesh ids.""" + if elements_dict is not None: - for elm_type, elements in elements_dict.items(): - self.mesh_dict[elm_type] = Mesh( - self.nodes, elements, elm_type, ndim=ndim + for key, value in elements_dict.items(): + submesh_name = "" + element_sets = None + if isinstance(value, Mesh): + elm_type = value.elm_type + elements = value.elements + element_sets = value.element_sets + submesh_name = value.name + elif isinstance(value, tuple): + if len(value) == 2: + elm_type, elements = value + elif len(value) == 3: + elm_type, elements, element_sets = value + else: + raise ValueError( + "MultiMesh element entries should be arrays, " + "(elm_type, elements), or " + "(elm_type, elements, element_sets)." + ) + submesh_name = str(key) + else: + elm_type = str(key) + elements = value + + self._mesh_list.append( + Mesh( + self.nodes, + elements, + elm_type, + element_sets=element_sets, + ndim=ndim, + name=submesh_name, + register_name=False, + ) ) - self.node_sets = {} + if node_sets is None: + node_sets = {} + self.node_sets = node_sets """Dict containing node sets associated to the mesh""" - def __getitem__(self, item: str) -> Mesh: - return self.mesh_dict[item] + def __getitem__(self, item: int | str) -> Mesh | "MultiMesh": + if isinstance(item, int): + return self._mesh_list[item] + elif isinstance(item, str): + mesh_list = [mesh for mesh in self._mesh_list if mesh.name == item] + if len(mesh_list) == 0: + mesh_list = [mesh for mesh in self._mesh_list if mesh.elm_type == item] + if len(mesh_list) == 0: + raise KeyError(item) + elif len(mesh_list) == 1: + return mesh_list[0] + else: + return MultiMesh.from_mesh_list( + mesh_list, + name=self.name, + node_sets=self.node_sets, + register_name=False, + ) + else: + raise TypeError( + "MultiMesh indices should be int, submesh name str, " + "or element type str." + ) def __repr__(self): return ( f"{object.__repr__(self)}\n\n" - f"elm_types: {tuple(self.mesh_dict.keys())}\n" + f"elm_types: {self.elm_types}\n" f"n_nodes: {self.n_nodes}\n" + f"n_submeshes: {len(self._mesh_list)}\n" f"n node_sets: {len(self.node_sets)}" ) @staticmethod - def from_mesh_list(mesh_list: list[Mesh], name: str = "") -> "MultiMesh": - multi_mesh = MultiMesh(mesh_list[0].nodes, name=name) - for mesh in mesh_list: - multi_mesh.mesh_dict[mesh.elm_type] = mesh + def from_mesh_list( + mesh_list: list[Mesh], + name: str = "", + node_sets: dict | None = None, + register_name: bool = True, + ) -> "MultiMesh": + """Build a MultiMesh from an ordered list of Mesh objects. + + The returned MultiMesh shares the node coordinates and element arrays + of the input meshes. Submesh names are copied from the input meshes but + are not registered globally. + """ + if len(mesh_list) == 0: + raise ValueError("Can't build a MultiMesh from an empty mesh list.") + if node_sets is None: + node_sets = mesh_list[0].node_sets + + return MultiMesh( + mesh_list[0].nodes, + {i: mesh for i, mesh in enumerate(mesh_list)}, + node_sets=node_sets, + name=name, + register_name=register_name, + ) + + @property + def submeshes(self) -> tuple[Mesh, ...]: + """Ordered view of the submeshes. + + The tuple order is the submesh id order used by integer indexing and + FDH5 export. + """ + return tuple(self._mesh_list) + + @property + def mesh_dict(self) -> dict[str, Mesh | "MultiMesh"]: + """Element-type indexed view kept for backward compatibility. + + If an element type appears once, the value is the matching ``Mesh``. + If an element type appears several times, the value is a filtered + ``MultiMesh`` containing all matching submeshes. + """ + mesh_dict = {} + for elm_type in dict.fromkeys(self.elm_types): + mesh_list = [ + mesh for mesh in self._mesh_list if mesh.elm_type == elm_type + ] + if len(mesh_list) == 1: + mesh_dict[elm_type] = mesh_list[0] + else: + mesh_dict[elm_type] = MultiMesh.from_mesh_list( + mesh_list, + name=self.name, + node_sets=self.node_sets, + register_name=False, + ) + return mesh_dict + + @property + def elm_types(self) -> tuple[str, ...]: + """Element type of each submesh, in submesh id order.""" + return tuple(mesh.elm_type for mesh in self._mesh_list) + + @property + def n_elements(self) -> int: + """Total number of elements in all submeshes.""" + return sum(mesh.n_elements for mesh in self._mesh_list) + + def copy(self, name: str = "") -> "MultiMesh": + """Make a shallow copy of the MultiMesh. + + The nodes, element tables, node sets and element sets are shared with + the original object. Submeshes are rebuilt as unregistered mesh views + that share the same node table. + """ + return MultiMesh.from_mesh_list( + self._mesh_list, + name=name, + node_sets=self.node_sets, + ) + + def deepcopy(self, name: str = "") -> "MultiMesh": + """Make a deep copy of the node and element arrays.""" + nodes = self.nodes.copy() + mesh_list = [ + Mesh( + nodes, + mesh.elements.copy(), + mesh.elm_type, + element_sets=mesh.element_sets, + ndim=self.ndim, + name=mesh.name, + register_name=False, + ) + for mesh in self._mesh_list + ] + return MultiMesh.from_mesh_list( + mesh_list, + name=name, + node_sets=self.node_sets, + ) + + def as_ndim(self, ndim, inplace: bool = False) -> "MultiMesh": + """Return a MultiMesh view with node coordinates in ``ndim`` dimensions.""" + if self.ndim == ndim: + return self + if self.ndim < ndim: + nodes = np.column_stack( + (self.nodes, np.zeros((self.n_nodes, ndim - self.ndim))) + ) + else: + nodes = self.nodes[:, :ndim] + + if inplace: + self.nodes = nodes + for mesh in self._mesh_list: + mesh.nodes = self.nodes + return self + + mesh_list = [ + Mesh( + nodes, + mesh.elements, + mesh.elm_type, + element_sets=mesh.element_sets, + ndim=ndim, + name=mesh.name, + register_name=False, + ) + for mesh in self._mesh_list + ] + return MultiMesh.from_mesh_list( + mesh_list, + name=self.name, + node_sets=self.node_sets, + register_name=False, + ) + + def as_2d(self, inplace: bool = False) -> "MultiMesh": + """Return a view of the current MultiMesh in 2D.""" + return self.as_ndim(2, inplace) + + def as_3d(self, inplace: bool = False) -> "MultiMesh": + """Return a view of the current MultiMesh in 3D.""" + return self.as_ndim(3, inplace) - return multi_mesh + def find_elements( + self, + selection_criterion: str, + value: float = 0, + tol: float = 1e-6, + select_by: str = "centers", + all_nodes: bool = True, + name: str | None = None, + ): + """Return global element ids matching a selection criterion.""" + element_ids = [] + offset = 0 + for mesh in self._mesh_list: + local_ids = mesh.find_elements( + selection_criterion, + value, + tol, + select_by, + all_nodes, + ) + element_ids.extend(np.asarray(local_ids, dtype=int) + offset) + offset += mesh.n_elements + + element_ids = np.asarray(element_ids, dtype=int) + if name: + raise NotImplementedError( + "Saving a new element set from MultiMesh.find_elements is not " + "implemented. Add the set to the relevant submesh instead." + ) + return element_ids + + @property + def element_sets(self) -> dict[str, np.ndarray]: + """Element sets using the global MultiMesh cell numbering. + + The returned indices follow the same cell order as + :meth:`to_pyvista`: all elements of submesh 0, then all elements of + submesh 1, and so on. If several submeshes define an element set with + the same name, the shifted indices are concatenated. + """ + element_sets = {} + offset = 0 + for mesh in self._mesh_list: + for name, values in mesh.element_sets.items(): + element_sets.setdefault(name, []).append( + np.asarray(values, dtype=int) + offset + ) + offset += mesh.n_elements + + return { + name: np.concatenate(values) + for name, values in element_sets.items() + } + + def to_pyvista(self, include_sets=False): + """Wrap the current multimesh to a pyvista UnstructuredGrid. + + If include_sets is True (default = False), the node and element sets + are included in the pyvista object as field data. Element set indices + are shifted to match the cell numbering of the generated mixed-element + grid. + """ + if not USE_PYVISTA: + raise NameError("Pyvista not installed.") + + if self.ndim < 3: + nodes = np.column_stack( + (self.nodes, np.zeros((self.n_nodes, 3 - self.ndim))) + ) + else: + nodes = self.nodes[:, :3] + + cells = [] + cell_types = [] + element_sets = {} + element_offset = 0 + + for mesh in self._mesh_list: + elm_type = mesh.elm_type + cell_type, n_elm_nodes = _get_pyvista_cell_info(elm_type) + + elm = np.empty((mesh.elements.shape[0], n_elm_nodes + 1), dtype=int) + elm[:, 0] = n_elm_nodes + elm[:, 1:] = mesh.elements[:, :n_elm_nodes] + cells.append(elm.ravel()) + cell_types.append(np.full(mesh.elements.shape[0], cell_type, dtype=int)) + + if include_sets: + for set_name, set_elements in mesh.element_sets.items(): + set_key = "_eset_" + set_name + element_sets.setdefault(set_key, []).append( + np.asarray(set_elements, dtype=int) + element_offset + ) + + element_offset += mesh.elements.shape[0] + + if not cells: + raise NotImplementedError("This mesh contains no compatible element.") + + pvmesh = pv.UnstructuredGrid( + np.concatenate(cells), + np.concatenate(cell_types), + nodes, + ) + + if self.ndim != 3: + pvmesh.field_data["ndim"] = self.ndim + + if include_sets: + for nset in self.node_sets: + pvmesh.field_data["_nset_" + nset] = self.node_sets[nset] + for set_key, set_elements in element_sets.items(): + pvmesh.field_data[set_key] = np.concatenate(set_elements) + + return pvmesh + + def _save_fdh5(self, filename: str) -> None: + """Save the multimesh definition to a raw FDH5 file.""" + from fedoo.util.fdh5 import FDH5Writer + + path = Path(filename) + if path.exists(): + path.unlink() + + writer = FDH5Writer(path) + writer.write_mesh(self.nodes, node_sets=self.node_sets, overwrite=True) + + for i, mesh in enumerate(self._mesh_list): + writer.add_submesh( + mesh.elm_type, + mesh.elements, + element_sets=mesh.element_sets, + name=mesh.name or mesh.elm_type, + submesh_id=f"submesh_{i}", + ) class BoundingBox(list): diff --git a/fedoo/core/output.py b/fedoo/core/output.py index d64f9d17..782224e5 100644 --- a/fedoo/core/output.py +++ b/fedoo/core/output.py @@ -2,6 +2,7 @@ # from fedoo.core.mesh import * from fedoo.core.base import AssemblyBase +from fedoo.core.mesh import MultiMesh # from fedoo.util.ExportData import ExportData from fedoo.core.dataset import DataSet, MultiFrameDataSet @@ -66,9 +67,198 @@ "npz", "csv", "xlsx", + "fdh5", ] +def _unique_assembly_mesh_entries(assemb, element_set=None): + """Return unique mesh entries for an AssemblySum. + + The returned list keeps the first occurrence order of meshes in + ``assemb.list_assembly``. If several elementary assemblies are attached to + the same mesh, they are grouped in the same entry and keep their original + assembly order. This order is later used to select the first assembly that + can provide a requested output field. + """ + entries = [] + mesh_id_to_entry = {} + + for sub_assemb in assemb.list_assembly: + mesh = sub_assemb.mesh + + if element_set is not None and isinstance(element_set, str): + if element_set not in mesh.element_sets: + continue + + mesh_id = id(mesh) + if mesh_id not in mesh_id_to_entry: + if element_set is None: + output_mesh = mesh + else: + output_mesh = mesh.extract_elements(element_set) + + mesh_id_to_entry[mesh_id] = len(entries) + entries.append( + { + "mesh": output_mesh, + "assemblies": [], + } + ) + + entries[mesh_id_to_entry[mesh_id]]["assemblies"].append(sub_assemb) + + return entries + + +def _output_mesh_for_assembly(assemb, element_set=None): + """Return the mesh associated with an output request. + + For a regular assembly, this is its mesh, optionally restricted to + ``element_set``. For an ``AssemblySum`` without ``assembly_output``, this + builds a ``MultiMesh`` from the unique elementary assembly meshes, in + ``list_assembly`` order. If all elementary assemblies share one mesh, the + mesh itself is returned. + """ + if hasattr(assemb, "list_assembly") and assemb.assembly_output is None: + entries = _unique_assembly_mesh_entries(assemb, element_set) + meshes = [entry["mesh"] for entry in entries] + if not meshes: + raise NameError("No mesh available for the requested element set.") + if len(meshes) == 1: + return meshes[0] + return MultiMesh.from_mesh_list( + meshes, + name=getattr(assemb, "name", ""), + register_name=False, + ) + + if hasattr(assemb, "list_assembly"): + assemb = assemb.assembly_output + + if element_set is None: + return assemb.mesh + return assemb.mesh.extract_elements(element_set) + + +def _dataset_field(data_set, field): + """Return a field and its data type from a DataSet, or (None, None).""" + if field in data_set.node_data: + return data_set.node_data[field], "Node" + if field in data_set.element_data: + return data_set.element_data[field], "Element" + if field in data_set.gausspoint_data: + return data_set.gausspoint_data[field], "GaussPoint" + if field in data_set.scalar_data: + return data_set.scalar_data[field], "Scalar" + return None, None + + +def _store_assemblysum_field(result, field, data, data_type, submesh_id, multimesh): + """Store one field block extracted from an AssemblySum. + + Element and Gauss point fields are stored per submesh when the AssemblySum + output mesh is a ``MultiMesh``. Node and scalar fields are shared by the + node list or by the problem and therefore keep the first available value. + """ + if data_type == "Node": + if field not in result.node_data: + result.node_data[field] = data + elif data_type == "Element": + if multimesh: + result.element_data.setdefault(field, {})[submesh_id] = data + else: + result.element_data[field] = data + elif data_type == "GaussPoint": + if multimesh: + result.gausspoint_data.setdefault(field, {})[submesh_id] = data + else: + result.gausspoint_data[field] = data + elif data_type == "Scalar": + if field not in result.scalar_data: + result.scalar_data[field] = data + + +def _get_assemblysum_results( + pb, + assemb, + output_list, + output_type=None, + position=1, + element_set=None, + include_mesh=True, +): + """Collect output data from an AssemblySum. + + Each requested field is searched independently on the elementary + assemblies. For each unique mesh, assemblies attached to that mesh are + tried in their ``list_assembly`` order and the first assembly providing + the field is used. If no assembly provides a field for a mesh, that mesh is + simply left without data for the field. If no assembly provides the field + at all, the field is omitted from the returned dataset. + + When several unique meshes are present, the returned ``DataSet`` is + associated with a ``MultiMesh`` and element/Gauss point fields are stored + as per-submesh dictionaries consumable by ``MultiMeshData``. + """ + entries = _unique_assembly_mesh_entries(assemb, element_set) + multimesh = len(entries) > 1 + + if include_mesh: + mesh = _output_mesh_for_assembly(assemb, element_set) + result = DataSet(mesh) + else: + result = DataSet() + + for res in output_list: + field_type = None + found = False + + for submesh_id, entry in enumerate(entries): + for sub_assemb in entry["assemblies"]: + try: + sub_result = _get_results( + pb, + sub_assemb, + [res], + output_type, + position, + element_set, + False, + ) + except NameError: + continue + + data, data_type = _dataset_field(sub_result, res) + if data_type is None: + continue + if field_type is None: + field_type = data_type + elif data_type != field_type: + raise NameError( + f'Field "{res}" has inconsistent output types in ' + "the AssemblySum." + ) + + _store_assemblysum_field( + result, + res, + data, + data_type, + submesh_id, + multimesh, + ) + found = True + break + + if not found: + continue + + if hasattr(pb, "time"): + result.scalar_data["Time"] = pb.time + + return result + + def _get_results( pb, assemb, @@ -96,6 +286,17 @@ def _get_results( if isinstance(assemb, str): assemb = AssemblyBase.get_all()[assemb] + if hasattr(assemb, "list_assembly") and assemb.assembly_output is None: + return _get_assemblysum_results( + pb, + assemb, + output_list, + output_type, + position, + element_set, + include_mesh, + ) + # for i, res in enumerate(output_list): # if ( # res not in _available_output @@ -112,20 +313,14 @@ def _get_results( data_sav = {} # dict to keep data in memory that may be used more that one time if hasattr(assemb, "list_assembly"): # AssemblySum object - if assemb.assembly_output is None: - raise NameError("AssemblySum objects can't be used to extract outputs") - else: - assemb = assemb.assembly_output + assemb = assemb.assembly_output sv = assemb.sv # state variables associated to the assembly if include_mesh: - if element_set is None: - result = DataSet(assemb.mesh) - else: - result = DataSet(assemb.mesh.extract_elements(element_set)) - if isinstance(element_set, str): - element_set = assemb.mesh.element_sets[element_set] + result = DataSet(_output_mesh_for_assembly(assemb, element_set)) + if element_set is not None and isinstance(element_set, str): + element_set = assemb.mesh.element_sets[element_set] else: result = DataSet() @@ -340,7 +535,7 @@ def add_output( assemb, output_list, output_type=None, - file_format="fdz", + file_format="fdh5", compressed=False, position=1, element_set=None, @@ -351,7 +546,7 @@ def add_output( extension = os.path.splitext(filename)[1] if extension == "": file_format = file_format.lower() - if file_format != "fdz": + if file_format not in ["fdz", "fdh5"]: # if no extention -> create a new dir using filename as dirname dirname = filename + "/" filename = dirname + os.path.basename(filename) @@ -386,10 +581,7 @@ def add_output( if isinstance(assemb, str): assemb = AssemblyBase.get_all()[assemb] - if element_set is None: - mesh = assemb.mesh - else: - mesh = assemb.mesh.extract_elements(element_set) + mesh = _output_mesh_for_assembly(assemb, element_set) if not (os.path.isdir(dirname)) and dirname != "": os.mkdir(dirname) @@ -415,7 +607,7 @@ def add_output( file.write("_mesh_.vtk") # add '_mesh_.vtk' to the zip archive os.remove("_mesh_.vtk") file.close() - elif save_mesh and (file_format not in ["vtk", "msh"]): + elif save_mesh and (file_format not in ["vtk", "msh", "fdh5"]): mesh.save(filename) res = MultiFrameDataSet(mesh, []) @@ -462,7 +654,7 @@ def save_results(self, pb, comp_output=None): list_file_format.append(file_format) list_compressed.append(compressed) - out = DataSet(assemb.mesh) + out = DataSet(_output_mesh_for_assembly(assemb, element_set)) list_data.append(out) else: # else, the same file is used diff --git a/fedoo/core/problem.py b/fedoo/core/problem.py index e4500e57..7e58a4b2 100644 --- a/fedoo/core/problem.py +++ b/fedoo/core/problem.py @@ -129,7 +129,7 @@ def add_output( assembly, output_list=None, output_type=None, - file_format="fdz", + file_format="fdh5", compressed=False, position=1, element_set=None, @@ -159,9 +159,9 @@ def add_output( Type of results. If None, the type of output is not converted. Scalar results are not concerned by this parameter. - file_format : "fdz", "vtk", "msh", "npz", "csv", "xlsx" + file_format : "fdh5", "fdz", "vtk", "msh", "npz", "csv", "xlsx" file format used to save the results. The default file format - and recommanding one is "fdz". + and recommanding one is "fdh5". compressed : bool, default = False if True, the fdz data are compressed. diff --git a/fedoo/util/fdh5.py b/fedoo/util/fdh5.py new file mode 100644 index 00000000..89bb97a0 --- /dev/null +++ b/fedoo/util/fdh5.py @@ -0,0 +1,1007 @@ +""" +The Fedoo FDH5 format is a HDF5 file that is structured +to store multiple data computed from fedoo. +The structure of the file is the following: + +mesh/ +├── nodes +├── node_sets/ +├── submesh_0/ +│ ├── elements +│ ├── element_sets/ +│ └── metadata/ +│ └── element_type = "tri3" +├── submesh_1/ +│ ├── elements +│ ├── element_sets/ +│ └── metadata/ +│ └── element_type = "tri3" +├── submesh_2/ +│ ├── elements +│ ├── element_sets/ +│ └── metadata/ +│ └── element_type = "quad4" +└── ... + +results/ +├── iter_0/ +│ ├── node_data/ +│ ├── element_data/ +│ │ ├── submesh_0/ +│ │ │ ├── stress +│ │ │ ├── plasticity +│ │ │ └── ... +│ │ ├── submesh_1/ +│ │ └── ... +│ ├── gausspoint_data/ +│ │ ├── submesh_0/ +│ │ ├── submesh_1/ +│ │ └── ... +│ ├── scalars/ +│ └── metadata/ +├── iter_1/ +└── ... +""" + + +from __future__ import annotations + +from pathlib import Path +from typing import Dict, Mapping, Optional, Union, Any, Literal, List, Iterator +from dataclasses import dataclass +from datetime import datetime, timezone +import contextlib + +import numpy as np +import h5py + +PathLike = Union[str, Path] +NDArray = np.ndarray + + +def submesh_id(index: int) -> str: + """Return the canonical FDH5 submesh id for a Fedoo submesh index.""" + return f"submesh_{index}" + + +def iteration_name(iteration: int) -> str: + """Return the canonical FDH5 group name for a result iteration.""" + return f"iter_{iteration}" + +@dataclass(frozen=True) +class CompressionConfig: + compression: Optional[Literal["gzip", "lzf"]] = "gzip" + compression_opts: Optional[int] = 4 + chunks: Optional[Union[bool, tuple[int, ...]]] = True + + +class FDH5Writer: + """ + Writer for Finite Element HDF5 files (FDH5). + + This writer implements the following structure: + + mesh/ + nodes + node_sets/ + submesh_X/ + elements + element_sets/ + metadata/ + + results/ + iter_0/ + node_data/ + element_data/submesh_X/ + gausspoint_data/submesh_X/ + scalars/ + metadata/ + """ + + FILE_VERSION = "1.0" + + def __init__( + self, + file_path: PathLike, + *, + compression: CompressionConfig = CompressionConfig(), + validate: bool = True, + create_parents: bool = True, + ) -> None: + self.path = Path(file_path) + if create_parents: + self.path.parent.mkdir(parents=True, exist_ok=True) + + self.compression = compression + self.validate = validate + + with h5py.File(self.path, "a") as f: + if "format" not in f.attrs: + f.attrs["format"] = "FDH5" + f.attrs["version"] = self.FILE_VERSION + f.attrs["created_utc"] = datetime.now(timezone.utc).isoformat() + + # ------------------------------------------------------------------ + # Internal helpers + # ------------------------------------------------------------------ + def _create_dataset( + self, + parent: h5py.Group, + name: str, + data: NDArray, + *, + overwrite: bool = False, + ) -> h5py.Dataset: + if name in parent: + if overwrite: + del parent[name] + else: + raise ValueError(f"Dataset '{parent.name}/{name}' already exists") + + kwargs: Dict[str, Any] = {} + is_scalar = np.asarray(data).shape == () + if self.compression.compression and not is_scalar: + kwargs["compression"] = self.compression.compression + if self.compression.compression == "gzip": + kwargs["compression_opts"] = self.compression.compression_opts + if self.compression.chunks is not None and not is_scalar: + kwargs["chunks"] = self.compression.chunks + + return parent.create_dataset(name, data=data, **kwargs) + + def _next_submesh_id(self, f: h5py.File) -> str: + existing = [k for k in f["mesh"].keys() if k.startswith("submesh_")] + if not existing: + return "submesh_0" + nums = [int(k.split("_")[1]) for k in existing] + return f"submesh_{max(nums) + 1}" + + # ------------------------------------------------------------------ + # Mesh writing + # ------------------------------------------------------------------ + def write_mesh( + self, + nodes: NDArray, + *, + node_sets: Optional[Mapping[str, NDArray]] = None, + overwrite: bool = False, + ) -> None: + nodes = np.asarray(nodes) + + if self.validate: + if nodes.ndim != 2: + raise ValueError("nodes must have shape (n_nodes, dim)") + if not np.issubdtype(nodes.dtype, np.floating): + raise TypeError("nodes must be floating point") + + with h5py.File(self.path, "a") as f: + mesh = f.require_group("mesh") + + if "nodes" in mesh and overwrite: + del mesh["nodes"] + + self._create_dataset(mesh, "nodes", nodes, overwrite=overwrite) + + if node_sets: + ns = mesh.require_group("node_sets") + if overwrite: + for k in list(ns.keys()): + del ns[k] + + for name, idx in node_sets.items(): + arr = np.asarray(idx) + self._create_dataset(ns, name, arr, overwrite=overwrite) + + def add_submesh( + self, + element_type: str, + elements: NDArray, + *, + element_sets: Optional[Mapping[str, NDArray]] = None, + name: Optional[str] = None, + submesh_id: Optional[str] = None, + overwrite: bool = False, + ) -> str: + elements = np.asarray(elements) + + if self.validate: + if elements.ndim != 2: + raise ValueError("elements must be 2D (n_elements, nodes_per_element)") + if not np.issubdtype(elements.dtype, np.integer): + raise TypeError("elements must be integer indices (0-based)") + + with h5py.File(self.path, "a") as f: + mesh = f.require_group("mesh") + sid = submesh_id or self._next_submesh_id(f) + + if sid in mesh: + if overwrite: + del mesh[sid] + else: + raise ValueError(f"Submesh '{sid}' already exists") + + sm = mesh.create_group(sid) + self._create_dataset(sm, "elements", elements) + + # Metadata + meta = sm.create_group("metadata") + meta.attrs["element_type"] = element_type + meta.attrs["n_elements"] = elements.shape[0] + meta.attrs["nodes_per_element"] = elements.shape[1] + if name: + meta.attrs["name"] = name + + # Element sets + if element_sets: + es = sm.create_group("element_sets") + for set_name, idx in element_sets.items(): + self._create_dataset(es, set_name, np.asarray(idx)) + + return sid + + # ------------------------------------------------------------------ + # Iteration writing + # ------------------------------------------------------------------ + def write_iteration( + self, + iteration: int, + *, + node_data: Optional[Mapping[str, NDArray]] = None, + element_data: Optional[Mapping[str, Mapping[str, NDArray]]] = None, + gausspoint_data: Optional[Mapping[str, Mapping[str, NDArray]]] = None, + scalars: Optional[Mapping[str, Union[int, float, NDArray]]] = None, + time: Optional[float] = None, + dt: Optional[float] = None, + overwrite: bool = False, + ) -> str: + iter_name = iteration_name(iteration) + + with h5py.File(self.path, "a") as f: + results = f.require_group("results") + + if iter_name in results: + if overwrite: + del results[iter_name] + else: + raise ValueError(f"Iteration '{iter_name}' already exists") + + it = results.create_group(iter_name) + + # Metadata + md = it.create_group("metadata") + md.attrs["iteration"] = iteration + if time is not None: + md.attrs["time"] = float(time) + if dt is not None: + md.attrs["dt"] = float(dt) + md.attrs["created_utc"] = datetime.now(timezone.utc).isoformat() + + # Node data + if node_data: + nd = it.create_group("node_data") + for name, arr in node_data.items(): + self._create_dataset(nd, name, np.asarray(arr)) + + # Element data + if element_data: + ed = it.create_group("element_data") + for sid, fields in element_data.items(): + smg = ed.create_group(sid) + for fname, arr in fields.items(): + self._create_dataset(smg, fname, np.asarray(arr)) + + # Gauss-point data (GP-major flattened) + if gausspoint_data: + gd = it.create_group("gausspoint_data") + for sid, fields in gausspoint_data.items(): + smg = gd.create_group(sid) + n_elements = f[f"mesh/{sid}/elements"].shape[0] + for fname, arr in fields.items(): + arr = np.asarray(arr) + ds = self._create_dataset(smg, fname, arr) + ds.attrs["gp_order"] = "gp-major" + ds.attrs["n_elements"] = n_elements + if arr.shape[0] % n_elements == 0: + ds.attrs["n_gauss_points"] = arr.shape[0] // n_elements + if arr.ndim == 1: + ds.attrs["n_components"] = 1 + else: + ds.attrs["n_components"] = arr.shape[1] + + # Scalars + if scalars: + sc = it.create_group("scalars") + for name, val in scalars.items(): + self._create_dataset(sc, name, np.asarray(val)) + + return it.name + + + + + +class FDH5Reader: + """ + Reader for FDH5 (Finite Element HDF5) files. + + Supports: + - eager reads (NumPy arrays) + - lazy reads (h5py.Dataset) + + Lazy access MUST be used inside `with reader.open():`. + """ + + def __init__(self, file_path: PathLike) -> None: + self.path = Path(file_path) + if not self.path.exists(): + raise FileNotFoundError(self.path) + + # ------------------------------------------------------------------ + # File handling + # ------------------------------------------------------------------ + @contextlib.contextmanager + def open(self) -> Iterator[h5py.File]: + """ + Open the HDF5 file. + + Required for lazy access (datasets remain valid only + while the file is open). + """ + f = h5py.File(self.path, "r") + try: + yield f + finally: + f.close() + + def _open(self) -> h5py.File: + """Internal helper for eager reads.""" + return h5py.File(self.path, "r") + + @staticmethod + def _decode(value: Any) -> Any: + if isinstance(value, (bytes, bytearray)): + return value.decode("utf-8") + return value + + def _require_lazy_file(self, file: h5py.File | None) -> h5py.File: + if file is None: + raise ValueError( + "Lazy reads require an open HDF5 file: " + "use `with reader.open() as f:` and pass `file=f`." + ) + return file + + # ------------------------------------------------------------------ + # Mesh + # ------------------------------------------------------------------ + def read_nodes( + self, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Union[NDArray, h5py.Dataset]: + if not lazy: + with self._open() as f: + return f["mesh/nodes"][...] + + f = self._require_lazy_file(file) + return f["mesh/nodes"] + + def read_node_sets( + self, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + if not lazy: + out: Dict[str, NDArray] = {} + with self._open() as f: + grp = f.get("mesh/node_sets") + if grp is None: + return out + for name, ds in grp.items(): + out[name] = ds[...] + return out + + f = self._require_lazy_file(file) + grp = f.get("mesh/node_sets") + return {} if grp is None else {name: ds for name, ds in grp.items()} + + def list_submeshes(self, *, file: h5py.File | None = None) -> List[str]: + if file is None: + with self._open() as f: + return self.list_submeshes(file=f) + + mesh = file.get("mesh") + if mesh is None: + return [] + return sorted(name for name in mesh if name.startswith("submesh_")) + + def read_submesh( + self, + submesh_id: str, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Any]: + """ + Read a submesh definition. + """ + if not lazy: + with self._open() as f: + sm = f[f"mesh/{submesh_id}"] + meta = sm["metadata"] + + data: Dict[str, Any] = { + "id": submesh_id, + "element_type": self._decode(meta.attrs.get("element_type")), + "elements": sm["elements"][...], + "element_sets": {}, + } + + if "name" in meta.attrs: + data["name"] = self._decode(meta.attrs["name"]) + + eset_grp = sm.get("element_sets") + if eset_grp is not None: + for name, ds in eset_grp.items(): + data["element_sets"][name] = ds[...] + + return data + + f = self._require_lazy_file(file) + sm = f[f"mesh/{submesh_id}"] + meta = sm["metadata"] + + return { + "id": submesh_id, + "element_type": self._decode(meta.attrs.get("element_type")), + "name": self._decode(meta.attrs.get("name")) if "name" in meta.attrs else None, + "elements": sm["elements"], + "element_sets": ( + {name: ds for name, ds in sm["element_sets"].items()} + if "element_sets" in sm + else {} + ), + } + + def read_mesh( + self, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Any]: + """ + Read the full mesh (nodes + submeshes). + """ + mesh: Dict[str, Any] = {} + mesh["nodes"] = self.read_nodes(lazy=lazy, file=file) + mesh["node_sets"] = self.read_node_sets(lazy=lazy, file=file) + + subs: Dict[str, Any] = {} + for sid in self.list_submeshes(file=file): + subs[sid] = self.read_submesh(sid, lazy=lazy, file=file) + + mesh["submeshes"] = subs + return mesh + + # ------------------------------------------------------------------ + # Iterations + # ------------------------------------------------------------------ + def list_iterations(self, *, file: h5py.File | None = None) -> List[int]: + if file is None: + with self._open() as f: + return self.list_iterations(file=f) + + grp = file.get("results") + if grp is None: + return [] + return sorted( + int(name.split("_")[1]) + for name in grp + if name.startswith("iter_") + ) + + def read_iteration_metadata(self, iteration: int) -> Dict[str, Any]: + with self._open() as f: + md = f[f"results/{iteration_name(iteration)}/metadata"] + return {k: self._decode(v) for k, v in md.attrs.items()} + + def read_iteration( + self, + iteration: int, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Any]: + """Read all data stored for one result iteration.""" + if lazy: + self._require_lazy_file(file) + + return { + "metadata": self.read_iteration_metadata(iteration), + "node_data": self.read_node_data(iteration, lazy=lazy, file=file), + "element_data": self.read_element_data(iteration, lazy=lazy, file=file), + "gausspoint_data": self.read_gausspoint_data( + iteration, lazy=lazy, file=file + ), + "scalars": self.read_scalars(iteration, lazy=lazy, file=file), + } + + def read_scalars( + self, + iteration: int, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + if not lazy: + out: Dict[str, NDArray] = {} + with self._open() as f: + grp = f.get(f"results/{iteration_name(iteration)}/scalars") + if grp is None: + return out + for name, ds in grp.items(): + out[name] = ds[...] + return out + + f = self._require_lazy_file(file) + grp = f.get(f"results/{iteration_name(iteration)}/scalars") + return {} if grp is None else {name: ds for name, ds in grp.items()} + + def read_node_data( + self, + iteration: int, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + if not lazy: + out: Dict[str, NDArray] = {} + with self._open() as f: + grp = f.get(f"results/{iteration_name(iteration)}/node_data") + if grp is None: + return out + for name, ds in grp.items(): + out[name] = ds[...] + return out + + f = self._require_lazy_file(file) + grp = f.get(f"results/{iteration_name(iteration)}/node_data") + return {} if grp is None else {name: ds for name, ds in grp.items()} + + def read_element_data( + self, + iteration: int, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Dict[str, Union[NDArray, h5py.Dataset]]]: + if not lazy: + out: Dict[str, Dict[str, NDArray]] = {} + with self._open() as f: + grp = f.get(f"results/{iteration_name(iteration)}/element_data") + if grp is None: + return out + for sid, sm_grp in grp.items(): + out[sid] = {name: ds[...] for name, ds in sm_grp.items()} + return out + + f = self._require_lazy_file(file) + grp = f.get(f"results/{iteration_name(iteration)}/element_data") + if grp is None: + return {} + return {sid: {name: ds for name, ds in sm_grp.items()} for sid, sm_grp in grp.items()} + + def read_gausspoint_data( + self, + iteration: int, + *, + lazy: bool = False, + file: h5py.File | None = None, + ) -> Dict[str, Dict[str, Union[NDArray, h5py.Dataset]]]: + if not lazy: + out: Dict[str, Dict[str, NDArray]] = {} + with self._open() as f: + grp = f.get(f"results/{iteration_name(iteration)}/gausspoint_data") + if grp is None: + return out + for sid, sm_grp in grp.items(): + out[sid] = {name: ds[...] for name, ds in sm_grp.items()} + return out + + f = self._require_lazy_file(file) + grp = f.get(f"results/{iteration_name(iteration)}/gausspoint_data") + if grp is None: + return {} + return {sid: {name: ds for name, ds in sm_grp.items()} for sid, sm_grp in grp.items()} + + +def mesh_to_fedoo(mesh_data: dict): + """Build a Fedoo Mesh or MultiMesh from FDH5 reader mesh data.""" + from fedoo.core.mesh import Mesh, MultiMesh + + nodes = mesh_data["nodes"] + node_sets = mesh_data.get("node_sets", {}) + submeshes = mesh_data.get("submeshes", {}) + + if len(submeshes) == 0: + return Mesh(nodes, node_sets=node_sets) + + ordered_submeshes = [ + submeshes[sid] + for sid in sorted(submeshes, key=lambda name: int(name.split("_")[1])) + ] + + if len(ordered_submeshes) == 1: + submesh = ordered_submeshes[0] + return Mesh( + nodes, + submesh["elements"], + submesh["element_type"], + node_sets=node_sets, + element_sets=submesh.get("element_sets", {}), + name=submesh.get("name", ""), + register_name=False, + ) + + elements_dict = {} + for i, submesh in enumerate(ordered_submeshes): + name = submesh.get("name") or submesh_id(i) + elements_dict[name] = ( + submesh["element_type"], + submesh["elements"], + submesh.get("element_sets", {}), + ) + + return MultiMesh( + nodes, + elements_dict, + node_sets=node_sets, + register_name=False, + ) + + +def fields_to_submesh_dict(dataset, fields: dict) -> dict: + """Convert Fedoo field dictionaries to the FDH5 submesh-first layout.""" + out = {} + + if dataset._is_multimesh(): + for field, value in fields.items(): + data = dataset._as_multimesh_data(value) + for sid, block in data.items(): + out.setdefault(submesh_id(sid), {})[field] = block + else: + for field, value in fields.items(): + out.setdefault("submesh_0", {})[field] = value + + return out + + +def fields_from_submesh_dict(mesh, fields: dict) -> dict: + """Convert FDH5 submesh-first data to Fedoo field dictionaries.""" + from fedoo.core.mesh import MultiMesh + + if not isinstance(mesh, MultiMesh): + if "submesh_0" in fields: + submesh_fields = fields["submesh_0"] + elif fields: + first_key = sorted(fields)[0] + submesh_fields = fields[first_key] + else: + submesh_fields = {} + return dict(submesh_fields) + + out = {} + for sid, submesh_fields in fields.items(): + submesh_index = int(sid.split("_")[1]) + for field, value in submesh_fields.items(): + out.setdefault(field, {})[submesh_index] = value + return out + + +def scalar_value(value): + """Return a scalar value as a Python scalar when possible.""" + arr = np.asarray(value) + if arr.shape == (): + return arr.item() + return arr + + +def load_dataset_iteration(dataset, filename: str, iteration: int = 0) -> None: + """Load one FDH5 iteration into an existing DataSet object.""" + reader = FDH5Reader(filename) + if dataset.mesh is None: + dataset.mesh = mesh_to_fedoo(reader.read_mesh()) + + iter_data = reader.read_iteration(iteration) + dataset.node_data = iter_data["node_data"] + dataset.element_data = fields_from_submesh_dict( + dataset.mesh, + iter_data["element_data"], + ) + dataset.gausspoint_data = fields_from_submesh_dict( + dataset.mesh, + iter_data["gausspoint_data"], + ) + dataset.scalar_data = { + key: scalar_value(value) + for key, value in iter_data["scalars"].items() + } + + +def write_dataset(dataset, filename: str, iteration: int = 0, overwrite: bool = False): + """Write a Fedoo DataSet iteration to a FDH5 file.""" + from fedoo.core.mesh import MultiMesh + + if dataset.mesh is None: + raise TypeError("Mesh should be defined before writing a FDH5 file.") + + path = Path(filename) + if path.suffix == "": + path = path.with_suffix(".fdh5") + + file_exists = path.exists() + if overwrite and file_exists: + path.unlink() + file_exists = False + + writer = FDH5Writer(path) + if not file_exists: + writer.write_mesh( + dataset.mesh.nodes, + node_sets=dataset.mesh.node_sets, + overwrite=True, + ) + + submeshes = ( + dataset.mesh.submeshes + if isinstance(dataset.mesh, MultiMesh) + else (dataset.mesh,) + ) + for i, mesh in enumerate(submeshes): + writer.add_submesh( + mesh.elm_type, + mesh.elements, + element_sets=mesh.element_sets, + name=mesh.name or mesh.elm_type, + submesh_id=submesh_id(i), + overwrite=True, + ) + + time = dataset.scalar_data.get("Time", None) + writer.write_iteration( + iteration, + node_data=dataset.node_data, + element_data=fields_to_submesh_dict(dataset, dataset.element_data), + gausspoint_data=fields_to_submesh_dict(dataset, dataset.gausspoint_data), + scalars=dataset.scalar_data, + time=time, + overwrite=True, + ) + + +def read_fdh5(filename: str): + """Read a FDH5 file as a DataSet or MultiFrameDataSet.""" + from fedoo.core.dataset import DataSet, MultiFrameDataSet + + path = Path(filename) + if path.suffix == "": + path = path.with_suffix(".fdh5") + + assert path.is_file(), "File not found" + + reader = FDH5Reader(path) + mesh = mesh_to_fedoo(reader.read_mesh()) + iterations = reader.list_iterations() + + if len(iterations) == 0: + return DataSet(mesh) + + if len(iterations) == 1: + dataset = DataSet(mesh) + load_dataset_iteration(dataset, str(path), iterations[0]) + return dataset + + dataset = MultiFrameDataSet(mesh) + dataset.list_data = [("fdh5", str(path), iteration) for iteration in iterations] + return dataset + + +# class FDH5File: +# # Class that allow both to write and read data +# # Don't know if this may be usefull +# def __init__(self, path, mode="r"): +# self.path = path +# self.mode = mode + +# if mode == "r": +# self.reader = FDH5Reader(path) +# self.writer = None +# elif mode in ("a", "w"): +# self.writer = FDH5Writer(path) +# self.reader = FDH5Reader(path) +# else: +# raise ValueError("mode must be 'r', 'a', or 'w'") + +# def write_iteration(self, *args, **kwargs): +# if self.writer is None: +# raise RuntimeError("File opened read-only") +# return self.writer.write_iteration(*args, **kwargs) + +# def read_node_data(self, *args, **kwargs): +# return self.reader.read_node_data(*args, **kwargs) + +# def open(self): +# return self.reader.open() + +if __name__ == "__main__": + # example of use + import numpy as np + from pathlib import Path + + # Import your classes + # from fdh5 import FDH5Writer, FDH5Reader + + # -------------------------------------------------- + # File path + # -------------------------------------------------- + path = Path("example.fdh5") + + writer = FDH5Writer(path, validate=True) + + # -------------------------------------------------- + # Mesh definition + # -------------------------------------------------- + + # 4 nodes, 2D + nodes = np.array([ + [0.0, 0.0], + [1.0, 0.0], + [1.0, 1.0], + [0.0, 1.0], + ], dtype=float) + + # Write mesh + node sets + writer.write_mesh( + nodes, + node_sets={ + "boundary": np.array([0, 1, 2, 3], dtype=int), + }, + overwrite=True, + ) + + # Two TRI3 elements + elements = np.array([ + [0, 1, 2], + [0, 2, 3], + ], dtype=int) + + submesh_id = writer.add_submesh( + element_type="tri3", + elements=elements, + name="square_triangles", + ) + + # -------------------------------------------------- + # Iteration 0 + # -------------------------------------------------- + + # Node field: displacement (n_nodes, 2) + node_data = { + "displacement": np.array([ + [0.0, 0.0], + [0.1, 0.0], + [0.1, 0.1], + [0.0, 0.1], + ], dtype=float) + } + + # Element field: von Mises stress (n_elements,) + element_data = { + submesh_id: { + "stress_vm": np.array([100.0, 120.0], dtype=float), + } + } + + # Gauss‑point field (GP‑major, flattened) + # Here: 2 elements, 2 Gauss points each + # Order: + # GP0: elem0, elem1 + # GP1: elem0, elem1 + gausspoint_data = { + submesh_id: { + "strain_eq": np.array([ + 0.01, 0.02, # GP0 + 0.015, 0.025 # GP1 + ], dtype=float).reshape(-1, 1) # (n_elem * n_gp, n_comp) + } + } + + # Scalars + scalars = { + "time": 0.0, + "total_energy": 42.0, + } + + # Write iteration + writer.write_iteration( + iteration=0, + node_data=node_data, + element_data=element_data, + gausspoint_data=gausspoint_data, + scalars=scalars, + time=0.0, + dt=0.1, + ) + + print("✅ File written:", path) + + + + + # ------------------------------------------------------------------------- + # Read written file + # ------------------------------------------------------------------------- + reader = FDH5Reader(path) + + # Read mesh + mesh = reader.read_mesh() + print("Nodes:\n", mesh["nodes"]) + print("Submeshes:", list(mesh["submeshes"].keys())) + + # Read iteration 0 + it0 = reader.read_iteration(0) + + u = it0["node_data"]["displacement"] + stress = it0["element_data"]["submesh_0"]["stress_vm"] + strain_gp = it0["gausspoint_data"]["submesh_0"]["strain_eq"] + + print("\nDisplacement:\n", u) + print("\nElement stress:\n", stress) + print("\nGauss-point strain (flattened):\n", strain_gp) + + + # ------------------------------------------------------------------------- + # Read written file lazily + # ------------------------------------------------------------------------- + reader = FDH5Reader(path) + + with reader.open() as f: + # Lazy node data + node_data = reader.read_node_data(0, lazy=True, file=f) + u_ds = node_data["displacement"] # h5py.Dataset + + # Only read first 2 nodes + print("First two displacements:\n", u_ds[:2]) + + # Lazy element data + elem_data = reader.read_element_data(0, lazy=True, file=f) + stress_ds = elem_data["submesh_0"]["stress_vm"] + + print("First element stress:", stress_ds[0]) + + # Lazy Gauss-point data + gp_data = reader.read_gausspoint_data(0, lazy=True, file=f) + eps_ds = gp_data["submesh_0"]["strain_eq"] + + # Infer GP layout manually + n_elems = f["mesh/submesh_0/elements"].shape[0] + + # GP0 for all elements + gp0 = eps_ds[0:n_elems] + print("GP0 strain:\n", gp0) + + # ------------------------------------------------------------------------- + # Read written file lazily + # ------------------------------------------------------------------------- + + from xdmf_writer import XDMFExporter # your exporter class + + exporter = XDMFExporter(Path("example.fdh5")) + exporter.export() + + import pyvista as pv + mesh = pv.read("example.xdmf") + mesh.plot(show_edges=True) + diff --git a/fedoo/util/viewer.py b/fedoo/util/viewer.py index 50ea4a5a..6217954b 100644 --- a/fedoo/util/viewer.py +++ b/fedoo/util/viewer.py @@ -26,9 +26,91 @@ import os import re +from fedoo.core.multimeshdata import MultiMeshData +from fedoo.core.mesh import MultiMesh + USE_PYVISTA_QT = True +def _is_multimesh(mesh): + return isinstance(mesh, MultiMesh) + + +def _global_data_array(data): + if isinstance(data, MultiMeshData): + return data.to_global() + return data + + +def _global_element_nodes(mesh, element_id): + """Return node ids for a global element id.""" + if not _is_multimesh(mesh): + return mesh.elements[element_id] + + submesh_id, local_id = _global_element_location(mesh, element_id) + return mesh[submesh_id].elements[local_id] + + +def _global_element_location(mesh, element_id): + """Return ``(submesh_id, local_element_id)`` for a global element id.""" + if not _is_multimesh(mesh): + return 0, element_id + + offset = 0 + for submesh_id, submesh in enumerate(mesh.submeshes): + stop = offset + submesh.n_elements + if offset <= element_id < stop: + return submesh_id, element_id - offset + offset = stop + raise IndexError(element_id) + + +def _element_data_value(data, field, component, data_type, element_id): + values = data.get_data(field, component, data_type) + values = _global_data_array(values) + return values[..., element_id] + + +def _gausspoint_values_for_element(data, field, component, element_id): + values = data.get_data(field, component, "GaussPoint") + if isinstance(values, MultiMeshData): + submesh_id, local_id = _global_element_location(data.mesh, element_id) + block = values.submesh(submesh_id) + if block is None: + return np.array([]) + block = np.asarray(block) + n_elements = data.mesh[submesh_id].n_elements + return block.reshape(-1, n_elements)[:, local_id] + + values = np.asarray(values) + return values.reshape(-1, data.mesh.n_elements)[:, element_id] + + +def _gausspoint_global_indices_for_element(data, field, component, element_id): + values = data.get_data(field, component, "GaussPoint") + if isinstance(values, MultiMeshData): + submesh_id, local_id = _global_element_location(data.mesh, element_id) + offset = 0 + for i in range(submesh_id): + block = values.submesh(i) + if block is None: + offset += data.mesh[i].n_elements + else: + offset += np.asarray(block).shape[-1] + + block = values.submesh(submesh_id) + if block is None: + return [] + block = np.asarray(block) + n_elements = data.mesh[submesh_id].n_elements + n_gp = block.shape[-1] // n_elements + return [offset + local_id + i * n_elements for i in range(n_gp)] + + values = np.asarray(values) + n_gp = values.shape[-1] // data.mesh.n_elements + return [element_id + i * data.mesh.n_elements for i in range(n_gp)] + + class DockTitleBar(QtWidgets.QWidget): clicked = Signal() @@ -212,9 +294,41 @@ def __init__(self, data, title, parent=None, opts=None): @property def pv_mesh(self, id_mesh=1): mesh_name = "data" + str(id_mesh) - if mesh_name not in self.plotter.actors: + if mesh_name in self.plotter.actors: + return pv.wrap(self.plotter.actors[mesh_name].GetMapper().GetInput()) + + def actor_sort_key(name): + try: + return (0, int(name.rsplit("_", 1)[1])) + except (IndexError, ValueError): + return (1, name) + + multimesh_actor_names = sorted( + (name for name in self.plotter.actors if name.startswith("data_")), + key=actor_sort_key, + ) + if not multimesh_actor_names: return None - return pv.wrap(self.plotter.actors[mesh_name].GetMapper().GetInput()) + + meshes = [ + pv.wrap(self.plotter.actors[name].GetMapper().GetInput()) + for name in multimesh_actor_names + ] + if len(meshes) == 1: + return meshes[0] + return pv.MultiBlock(meshes).combine() + + def picked_actor_mesh(self, actor): + """Return the actor name and PyVista mesh for a picked VTK actor.""" + for name, candidate in self.plotter.actors.items(): + if candidate is actor: + return name, pv.wrap(candidate.GetMapper().GetInput()) + try: + if candidate.GetAddressAsString("") == actor.GetAddressAsString(""): + return name, pv.wrap(candidate.GetMapper().GetInput()) + except AttributeError: + pass + return None, None def get_components(self, field): if field == "": @@ -274,6 +388,10 @@ def update_plot(self, val=None, iteration=None, lock_view=True, plotter=None): } # plotter.clear() # not compatible with pbr ??? plotter.renderer.clear_actors() + multimesh_kargs = {} + if _is_multimesh(self.data.mesh): + multimesh_kargs["global_element_set"] = True + multimesh_kargs["name"] = "data" self.data.plot( field=self.current_field, @@ -299,6 +417,7 @@ def update_plot(self, val=None, iteration=None, lock_view=True, plotter=None): element_set=self.opts["element_set"], # element_set_invert = self.opts["element_set_invert"], cmap=self.opts["cmap"], + **multimesh_kargs, ) if self.parent()._plane_widget_enabled: @@ -1261,7 +1380,7 @@ def open_file(self): self, "Open file", "", - "Fedoo files (*.fdz) ;;VTK Files (*.vtk);;CSV Files (*.csv) ;; All Files (*)", + "Fedoo files (*.fdz *.fdh5) ;;VTK Files (*.vtk);;CSV Files (*.csv) ;; All Files (*)", ) if fname: data = fd.read_data(fname) @@ -1937,7 +2056,10 @@ def _picked(point): pid = self.data.mesh.nearest_node(point) else: pid = mesh.find_closest_point(point) - if self.current_data_type == "GaussPoint": + if ( + self.current_data_type == "GaussPoint" + and not _is_multimesh(self.data.mesh) + ): pid = self.data.mesh.elements.ravel()[pid] except Exception: pid = None @@ -1948,7 +2070,7 @@ def _picked(point): msg = "" msg += f"x={point[0]:.3g}, y={point[1]:.3g}, z={point[2]:.3g}" name = self.current_field - if name: + if name and pid is not None: msg += ( f" | {name}={self.data[name, self.current_component, 'Node'][pid]}" ) @@ -1974,7 +2096,10 @@ def _on_click(obj, event): cell_id = picker.GetCellId() if cell_id < 0: return - mesh = self.active_dock.pv_mesh + display_cell_id = cell_id + actor_name, mesh = self.active_dock.picked_actor_mesh(picker.GetActor()) + if mesh is None: + mesh = self.active_dock.pv_mesh if "vtkOriginalCellIds" in mesh.cell_data: save_original_cell_ids = mesh.cell_data["vtkOriginalCellIds"] cell = mesh.extract_cells(cell_id) @@ -1992,17 +2117,35 @@ def _on_click(obj, event): elif "vtkOriginalCellIds" in mesh.cell_data: # don't work !!!! cell_id = mesh.cell_data["vtkOriginalCellIds"][cell_id] + if "_fedoo_global_cell_ids" in mesh.cell_data: + cell_id = mesh.cell_data["_fedoo_global_cell_ids"][display_cell_id] + cell_id = int(cell_id) + msg = "" if mesh is not None and cell_id is not None: msg = f"Element id={cell_id} | " # msg += f"x={point[0]:.3g}, y={point[1]:.3g}, z={point[2]:.3g}" name = self.current_field - data_gp = self.data[name, self.current_component, "GaussPoint"].reshape( - -1, self.data.mesh.n_elements - )[:, cell_id] if name: - msg += f"{name}_{self.current_component}={self.data[name, self.current_component, 'Element'][cell_id]}" - msg += f" | gp vals: {data_gp}" + try: + elem_value = _element_data_value( + self.data, + name, + self.current_component, + "Element", + cell_id, + ) + msg += f"{name}_{self.current_component}={elem_value}" + except Exception: + pass + + try: + data_gp = _gausspoint_values_for_element( + self.data, name, self.current_component, cell_id + ) + msg += f" | gp vals: {data_gp}" + except Exception: + pass self.statusBar().showMessage(msg) @@ -2379,7 +2522,8 @@ def _update_clim_values(self): parent = self.parent() if self.rb_current.isChecked(): data = parent.get_current_data() - clim = [data.min(), data.max()] + data = _global_data_array(data) + clim = [np.nanmin(data), np.nanmax(data)] elif self.rb_all.isChecked(): if hasattr(parent.active_dock.data, "get_all_frame_lim"): clim = parent.active_dock.data.get_all_frame_lim( @@ -2389,7 +2533,8 @@ def _update_clim_values(self): )[2] else: data = parent.get_current_data() - clim = [data.min(), data.max()] + data = _global_data_array(data) + clim = [np.nanmin(data), np.nanmax(data)] self.vmin_spin.setValue(float(clim[0])) self.vmax_spin.setValue(float(clim[1])) @@ -3128,7 +3273,9 @@ def _on_rectangle(rect_selection): dtype=bool, ) picked_ids = np.nonzero(inside_mask)[0].astype(int) - if "vtkOriginalCellIds" in pvmesh.cell_data: + if "_fedoo_global_cell_ids" in pvmesh.cell_data: + picked_ids = pvmesh.cell_data["_fedoo_global_cell_ids"][picked_ids] + elif "vtkOriginalCellIds" in pvmesh.cell_data: picked_ids = pvmesh.cell_data["vtkOriginalCellIds"][picked_ids] picked_ids = set(picked_ids) @@ -3186,6 +3333,9 @@ def _apply_expression_to_selection(self): if "Disp" in data.node_data and self.use_def_mesh.isChecked(): mesh = data.mesh.copy() mesh.nodes = mesh.nodes + data.node_data["Disp"].T + if _is_multimesh(mesh): + for submesh in mesh.submeshes: + submesh.nodes = mesh.nodes else: mesh = data.mesh ids = set(mesh.find_elements(expr)) @@ -3612,7 +3762,10 @@ def _picked(point): pid = mainwin.data.mesh.nearest_node(point) else: pid = mesh.find_closest_point(point) - if mainwin.current_data_type == "GaussPoint": + if ( + mainwin.current_data_type == "GaussPoint" + and not _is_multimesh(mainwin.data.mesh) + ): pid = mainwin.data.mesh.elements.ravel()[pid] if mesh is not None and pid is not None and pid >= 0: @@ -3637,7 +3790,12 @@ def _on_click(obj, event): cell_id = picker.GetCellId() if cell_id < 0: return - mesh = mainwin.active_dock.pv_mesh + display_cell_id = cell_id + actor_name, mesh = mainwin.active_dock.picked_actor_mesh( + picker.GetActor() + ) + if mesh is None: + mesh = mainwin.active_dock.pv_mesh cell = mesh.extract_cells(cell_id) mainwin.plotter.remove_actor("_picked_cell") mainwin.plotter.add_mesh( @@ -3645,8 +3803,12 @@ def _on_click(obj, event): ) if mainwin.opts["clip_args"]: cell_id = mesh.cell_data["cell_ids"][cell_id] + elif "vtkOriginalCellIds" in mesh.cell_data: + cell_id = mesh.cell_data["vtkOriginalCellIds"][cell_id] + if "_fedoo_global_cell_ids" in mesh.cell_data: + cell_id = mesh.cell_data["_fedoo_global_cell_ids"][display_cell_id] if mesh is not None and cell_id is not None: - self.id_spin.setValue(cell_id) + self.id_spin.setValue(int(cell_id)) # Remove observer after pick if hasattr(self, "_lbp_tag") and self._lbp_tag is not None: mainwin.plotter.interactor.RemoveObserver(self._lbp_tag) @@ -3790,7 +3952,9 @@ def check_valide(self, field, comp, data_type, indice): if indice >= data.mesh.n_elements: return False else: # GaussPoint - if indice >= data[field, comp, "GaussPoint"].shape[-1]: + gp_data = data.get_data(field, comp, "GaussPoint") + gp_data = _global_data_array(gp_data) + if indice >= gp_data.shape[-1]: return False return True @@ -3815,7 +3979,7 @@ def add_y_data(self, **kargs): field = self.field_combo.currentText() comp = self.comp_combo.currentText() dtype = self.data_type_combo.currentText() - if not indice: + if indice is None: idx = self.id_spin.value() else: idx = indice @@ -3827,10 +3991,10 @@ def add_y_data(self, **kargs): else: if indice is None and dtype == "GaussPoint": # add all gauss points values - n_elements = data.mesh.n_elements - n_gp = data[field, comp, "GaussPoint"].shape[-1] // n_elements - for i in range(n_gp): - self.add_y_data(indice=idx + i * n_elements) + for gp_id in _gausspoint_global_indices_for_element( + data, field, comp, idx + ): + self.add_y_data(indice=gp_id) return label = f"{field}_{comp} ({dtype}, ID={idx})" # Add to list and table diff --git a/fedoo/util/xdmf_writer.py b/fedoo/util/xdmf_writer.py new file mode 100644 index 00000000..3c339292 --- /dev/null +++ b/fedoo/util/xdmf_writer.py @@ -0,0 +1,267 @@ + +from __future__ import annotations + +from pathlib import Path +from typing import Iterable, Optional, Any +from lxml import etree + +import h5py + +from .fdh5 import iteration_name + + +class XDMFExporter: + """ + Export FDH5 files to XDMF (XDMF 3.0, HDF5-backed). + + - One XDMF file for the full time series + - Multiple submeshes supported + - Node, element and Gauss-point data + - Gauss-point data exported in flattened GP-major form + """ + + XDMF_VERSION = "3.0" + + # FE element type -> (XDMF topology type, nodes per element) + ELEMENT_TYPE_MAP = { + "tri3": ("Triangle", 3), + "quad4": ("Quadrilateral", 4), + "tet4": ("Tetrahedron", 4), + "hex8": ("Hexahedron", 8), + } + + def __init__( + self, + h5_path: Path, + xdmf_path: Optional[Path] = None, + ) -> None: + self.h5_path = Path(h5_path) + self.xdmf_path = xdmf_path or self.h5_path.with_suffix(".xdmf") + + # ------------------------------------------------------------------ + # Public API + # ------------------------------------------------------------------ + def export(self) -> Path: + """ + Generate the XDMF file for all iterations (time series). + """ + root = etree.Element("Xdmf", Version=self.XDMF_VERSION) + domain = etree.SubElement(root, "Domain") + + temporal = etree.SubElement( + domain, + "Grid", + Name="TimeSeries", + GridType="Collection", + CollectionType="Temporal", + ) + + with h5py.File(self.h5_path, "r") as f: + for it in self._list_iterations(f): + temporal.append(self._build_iteration_grid(f, it)) + + tree = etree.ElementTree(root) + tree.write( + str(self.xdmf_path), + pretty_print=True, + xml_declaration=True, + encoding="UTF-8", + ) + + return self.xdmf_path + + # ------------------------------------------------------------------ + # Iterations + # ------------------------------------------------------------------ + def _list_iterations(self, f: h5py.File) -> Iterable[int]: + grp = f.get("results") + if grp is None: + return [] + return sorted( + int(name.split("_")[1]) + for name in grp + if name.startswith("iter_") + ) + + def _build_iteration_grid(self, f: h5py.File, iteration: int) -> etree.Element: + it_name = iteration_name(iteration) + it_grp = f[f"results/{it_name}"] + + time_value = it_grp["metadata"].attrs.get("time", iteration) + + it_grid = etree.Element( + "Grid", + Name=it_name, + GridType="Collection", + CollectionType="Spatial", + ) + etree.SubElement(it_grid, "Time", Value=str(time_value)) + + for submesh_id in self._list_submeshes(f): + it_grid.append(self._build_submesh_grid(f, it_grp, submesh_id)) + + return it_grid + + # ------------------------------------------------------------------ + # Submeshes + # ------------------------------------------------------------------ + def _list_submeshes(self, f: h5py.File) -> Iterable[str]: + mesh = f.get("mesh") + if mesh is None: + return [] + return sorted(name for name in mesh if name.startswith("submesh_")) + + def _build_submesh_grid( + self, + f: h5py.File, + it_grp: h5py.Group, + submesh_id: str, + ) -> etree.Element: + sm = f[f"mesh/{submesh_id}"] + meta = sm["metadata"] + + element_type = meta.attrs.get("element_type") + if isinstance(element_type, (bytes, bytearray)): + element_type = element_type.decode("utf-8") + + if element_type not in self.ELEMENT_TYPE_MAP: + raise ValueError(f"Unsupported element type '{element_type}'") + + topo_type, nodes_per_elem = self.ELEMENT_TYPE_MAP[element_type] + + elements = sm["elements"] + n_elems = elements.shape[0] + + grid = etree.Element("Grid", Name=submesh_id, GridType="Uniform") + + # ------------------------------ + # Topology + # ------------------------------ + topo = etree.SubElement( + grid, + "Topology", + TopologyType=topo_type, + NumberOfElements=str(n_elems), + ) + + etree.SubElement( + topo, + "DataItem", + Dimensions=f"{n_elems} {nodes_per_elem}", + NumberType="Int", + Format="HDF", + ).text = f"{self.h5_path.name}:/mesh/{submesh_id}/elements" + + # ------------------------------ + # Geometry (global nodes) + # ------------------------------ + nodes = f["mesh/nodes"] + n_nodes, dim = nodes.shape + geom_type = "XYZ" if dim == 3 else "XY" + + geom = etree.SubElement( + grid, + "Geometry", + GeometryType=geom_type, + ) + + etree.SubElement( + geom, + "DataItem", + Dimensions=f"{n_nodes} {dim}", + NumberType="Float", + Precision="8", + Format="HDF", + ).text = f"{self.h5_path.name}:/mesh/nodes" + + # ------------------------------ + # Node data (global) + # ------------------------------ + node_grp = it_grp.get("node_data") + if node_grp is not None: + for name, ds in node_grp.items(): + self._add_attribute( + grid, + name=name, + center="Node", + dataset_path=f"/results/{it_grp.name.split('/')[-1]}/node_data/{name}", + shape=ds.shape, + ) + + # ------------------------------ + # Element data (per submesh) + # ------------------------------ + elem_grp = it_grp.get(f"element_data/{submesh_id}") + if elem_grp is not None: + for name, ds in elem_grp.items(): + self._add_attribute( + grid, + name=name, + center="Cell", + dataset_path=( + f"/results/{it_grp.name.split('/')[-1]}/" + f"element_data/{submesh_id}/{name}" + ), + shape=ds.shape, + ) + + # ------------------------------ + # Gauss-point data (flattened GP-major) + # Exported as Cell-centered attributes + # ------------------------------ + gp_grp = it_grp.get(f"gausspoint_data/{submesh_id}") + if gp_grp is not None: + for name, ds in gp_grp.items(): + n_gauss_points = int(ds.attrs.get("n_gauss_points", 1)) + if n_gauss_points != 1: + continue + self._add_attribute( + grid, + name=name, + center="Cell", + dataset_path=( + f"/results/{it_grp.name.split('/')[-1]}/" + f"gausspoint_data/{submesh_id}/{name}" + ), + shape=ds.shape, + ) + + return grid + + # ------------------------------------------------------------------ + # Attributes + # ------------------------------------------------------------------ + def _add_attribute( + self, + grid: etree.Element, + *, + name: str, + center: str, + dataset_path: str, + shape: Any, + ) -> None: + """ + Add an Attribute element referencing an HDF5 dataset. + """ + attr_type = "Scalar" + if len(shape) == 2 and shape[1] == 3: + attr_type = "Vector" + elif len(shape) == 2 and shape[1] == 6: + attr_type = "Tensor6" + + attr = etree.SubElement( + grid, + "Attribute", + Name=name, + AttributeType=attr_type, + Center=center, + ) + + etree.SubElement( + attr, + "DataItem", + Dimensions=" ".join(str(s) for s in shape), + NumberType="Float", + Precision="8", + Format="HDF", + ).text = f"{self.h5_path.name}:{dataset_path}" From 010fcdfb8776ac2681ec1e90f79394a29beb2623 Mon Sep 17 00:00:00 2001 From: pruliere Date: Mon, 6 Jul 2026 09:38:11 +0200 Subject: [PATCH 02/14] ruff format --- .../01-simple/spherical_shell_compression.py | 2 +- fedoo/core/dataset.py | 15 +- fedoo/core/mesh.py | 9 +- fedoo/util/fdh5.py | 142 ++++++++++-------- fedoo/util/viewer.py | 14 +- fedoo/util/xdmf_writer.py | 5 +- 6 files changed, 97 insertions(+), 90 deletions(-) diff --git a/examples/01-simple/spherical_shell_compression.py b/examples/01-simple/spherical_shell_compression.py index b1faad20..3195dab3 100644 --- a/examples/01-simple/spherical_shell_compression.py +++ b/examples/01-simple/spherical_shell_compression.py @@ -68,7 +68,7 @@ # Without them, the unconstrained problem is singular; some solvers may still # return a usable strain/stress field, but the displacement field can contain # arbitrary rigid-body motion. -# +# pb = fd.problem.Linear(assembly) diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 46fa6569..2a040b93 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -940,9 +940,12 @@ def _submesh_element_set( return element_ids, True offset = offsets[submesh_id] - local_ids = element_ids[ - (element_ids >= offset) & (element_ids < offset + submesh.n_elements) - ] - offset + local_ids = ( + element_ids[ + (element_ids >= offset) & (element_ids < offset + submesh.n_elements) + ] + - offset + ) if len(local_ids) == 0 and not element_set_invert: return None, False if len(local_ids) == 0 and element_set_invert: @@ -1871,7 +1874,11 @@ def load(self, data=-1, load_mesh=False): if iteration > len(self.list_data) or iteration < 0: raise NameError("Number of iteration out of bounds") data_ref = self.list_data[iteration] - if isinstance(data_ref, tuple) and len(data_ref) == 3 and data_ref[0] == "fdh5": + if ( + isinstance(data_ref, tuple) + and len(data_ref) == 3 + and data_ref[0] == "fdh5" + ): DataSet.load(self, data_ref[1], load_mesh, iteration=data_ref[2]) else: DataSet.load(self, data_ref) diff --git a/fedoo/core/mesh.py b/fedoo/core/mesh.py index 0227af94..06beb64e 100644 --- a/fedoo/core/mesh.py +++ b/fedoo/core/mesh.py @@ -1885,9 +1885,7 @@ def mesh_dict(self) -> dict[str, Mesh | "MultiMesh"]: """ mesh_dict = {} for elm_type in dict.fromkeys(self.elm_types): - mesh_list = [ - mesh for mesh in self._mesh_list if mesh.elm_type == elm_type - ] + mesh_list = [mesh for mesh in self._mesh_list if mesh.elm_type == elm_type] if len(mesh_list) == 1: mesh_dict[elm_type] = mesh_list[0] else: @@ -2036,10 +2034,7 @@ def element_sets(self) -> dict[str, np.ndarray]: ) offset += mesh.n_elements - return { - name: np.concatenate(values) - for name, values in element_sets.items() - } + return {name: np.concatenate(values) for name, values in element_sets.items()} def to_pyvista(self, include_sets=False): """Wrap the current multimesh to a pyvista UnstructuredGrid. diff --git a/fedoo/util/fdh5.py b/fedoo/util/fdh5.py index 89bb97a0..812a3db4 100644 --- a/fedoo/util/fdh5.py +++ b/fedoo/util/fdh5.py @@ -43,7 +43,6 @@ └── ... """ - from __future__ import annotations from pathlib import Path @@ -68,6 +67,7 @@ def iteration_name(iteration: int) -> str: """Return the canonical FDH5 group name for a result iteration.""" return f"iter_{iteration}" + @dataclass(frozen=True) class CompressionConfig: compression: Optional[Literal["gzip", "lzf"]] = "gzip" @@ -317,9 +317,6 @@ def write_iteration( return it.name - - - class FDH5Reader: """ Reader for FDH5 (Finite Element HDF5) files. @@ -456,7 +453,9 @@ def read_submesh( return { "id": submesh_id, "element_type": self._decode(meta.attrs.get("element_type")), - "name": self._decode(meta.attrs.get("name")) if "name" in meta.attrs else None, + "name": self._decode(meta.attrs.get("name")) + if "name" in meta.attrs + else None, "elements": sm["elements"], "element_sets": ( {name: ds for name, ds in sm["element_sets"].items()} @@ -497,9 +496,7 @@ def list_iterations(self, *, file: h5py.File | None = None) -> List[int]: if grp is None: return [] return sorted( - int(name.split("_")[1]) - for name in grp - if name.startswith("iter_") + int(name.split("_")[1]) for name in grp if name.startswith("iter_") ) def read_iteration_metadata(self, iteration: int) -> Dict[str, Any]: @@ -591,7 +588,10 @@ def read_element_data( grp = f.get(f"results/{iteration_name(iteration)}/element_data") if grp is None: return {} - return {sid: {name: ds for name, ds in sm_grp.items()} for sid, sm_grp in grp.items()} + return { + sid: {name: ds for name, ds in sm_grp.items()} + for sid, sm_grp in grp.items() + } def read_gausspoint_data( self, @@ -614,7 +614,10 @@ def read_gausspoint_data( grp = f.get(f"results/{iteration_name(iteration)}/gausspoint_data") if grp is None: return {} - return {sid: {name: ds for name, ds in sm_grp.items()} for sid, sm_grp in grp.items()} + return { + sid: {name: ds for name, ds in sm_grp.items()} + for sid, sm_grp in grp.items() + } def mesh_to_fedoo(mesh_data: dict): @@ -725,8 +728,7 @@ def load_dataset_iteration(dataset, filename: str, iteration: int = 0) -> None: iter_data["gausspoint_data"], ) dataset.scalar_data = { - key: scalar_value(value) - for key, value in iter_data["scalars"].items() + key: scalar_value(value) for key, value in iter_data["scalars"].items() } @@ -839,29 +841,32 @@ def read_fdh5(filename: str): # example of use import numpy as np from pathlib import Path - + # Import your classes # from fdh5 import FDH5Writer, FDH5Reader - + # -------------------------------------------------- # File path # -------------------------------------------------- path = Path("example.fdh5") - + writer = FDH5Writer(path, validate=True) - + # -------------------------------------------------- # Mesh definition # -------------------------------------------------- - + # 4 nodes, 2D - nodes = np.array([ - [0.0, 0.0], - [1.0, 0.0], - [1.0, 1.0], - [0.0, 1.0], - ], dtype=float) - + nodes = np.array( + [ + [0.0, 0.0], + [1.0, 0.0], + [1.0, 1.0], + [0.0, 1.0], + ], + dtype=float, + ) + # Write mesh + node sets writer.write_mesh( nodes, @@ -870,40 +875,46 @@ def read_fdh5(filename: str): }, overwrite=True, ) - + # Two TRI3 elements - elements = np.array([ - [0, 1, 2], - [0, 2, 3], - ], dtype=int) - + elements = np.array( + [ + [0, 1, 2], + [0, 2, 3], + ], + dtype=int, + ) + submesh_id = writer.add_submesh( element_type="tri3", elements=elements, name="square_triangles", ) - + # -------------------------------------------------- # Iteration 0 # -------------------------------------------------- - + # Node field: displacement (n_nodes, 2) node_data = { - "displacement": np.array([ - [0.0, 0.0], - [0.1, 0.0], - [0.1, 0.1], - [0.0, 0.1], - ], dtype=float) + "displacement": np.array( + [ + [0.0, 0.0], + [0.1, 0.0], + [0.1, 0.1], + [0.0, 0.1], + ], + dtype=float, + ) } - + # Element field: von Mises stress (n_elements,) element_data = { submesh_id: { "stress_vm": np.array([100.0, 120.0], dtype=float), } } - + # Gauss‑point field (GP‑major, flattened) # Here: 2 elements, 2 Gauss points each # Order: @@ -911,19 +922,24 @@ def read_fdh5(filename: str): # GP1: elem0, elem1 gausspoint_data = { submesh_id: { - "strain_eq": np.array([ - 0.01, 0.02, # GP0 - 0.015, 0.025 # GP1 - ], dtype=float).reshape(-1, 1) # (n_elem * n_gp, n_comp) + "strain_eq": np.array( + [ + 0.01, + 0.02, # GP0 + 0.015, + 0.025, # GP1 + ], + dtype=float, + ).reshape(-1, 1) # (n_elem * n_gp, n_comp) } } - + # Scalars scalars = { "time": 0.0, "total_energy": 42.0, } - + # Write iteration writer.write_iteration( iteration=0, @@ -934,11 +950,8 @@ def read_fdh5(filename: str): time=0.0, dt=0.1, ) - - print("✅ File written:", path) - - + print("✅ File written:", path) # ------------------------------------------------------------------------- # Read written file @@ -949,19 +962,18 @@ def read_fdh5(filename: str): mesh = reader.read_mesh() print("Nodes:\n", mesh["nodes"]) print("Submeshes:", list(mesh["submeshes"].keys())) - + # Read iteration 0 it0 = reader.read_iteration(0) - + u = it0["node_data"]["displacement"] stress = it0["element_data"]["submesh_0"]["stress_vm"] strain_gp = it0["gausspoint_data"]["submesh_0"]["strain_eq"] - + print("\nDisplacement:\n", u) print("\nElement stress:\n", stress) print("\nGauss-point strain (flattened):\n", strain_gp) - - + # ------------------------------------------------------------------------- # Read written file lazily # ------------------------------------------------------------------------- @@ -970,24 +982,24 @@ def read_fdh5(filename: str): with reader.open() as f: # Lazy node data node_data = reader.read_node_data(0, lazy=True, file=f) - u_ds = node_data["displacement"] # h5py.Dataset - + u_ds = node_data["displacement"] # h5py.Dataset + # Only read first 2 nodes print("First two displacements:\n", u_ds[:2]) - + # Lazy element data elem_data = reader.read_element_data(0, lazy=True, file=f) stress_ds = elem_data["submesh_0"]["stress_vm"] - + print("First element stress:", stress_ds[0]) - + # Lazy Gauss-point data gp_data = reader.read_gausspoint_data(0, lazy=True, file=f) eps_ds = gp_data["submesh_0"]["strain_eq"] - + # Infer GP layout manually n_elems = f["mesh/submesh_0/elements"].shape[0] - + # GP0 for all elements gp0 = eps_ds[0:n_elems] print("GP0 strain:\n", gp0) @@ -997,11 +1009,11 @@ def read_fdh5(filename: str): # ------------------------------------------------------------------------- from xdmf_writer import XDMFExporter # your exporter class - + exporter = XDMFExporter(Path("example.fdh5")) exporter.export() - + import pyvista as pv + mesh = pv.read("example.xdmf") mesh.plot(show_edges=True) - diff --git a/fedoo/util/viewer.py b/fedoo/util/viewer.py index 6217954b..fc7adb92 100644 --- a/fedoo/util/viewer.py +++ b/fedoo/util/viewer.py @@ -2056,9 +2056,8 @@ def _picked(point): pid = self.data.mesh.nearest_node(point) else: pid = mesh.find_closest_point(point) - if ( - self.current_data_type == "GaussPoint" - and not _is_multimesh(self.data.mesh) + if self.current_data_type == "GaussPoint" and not _is_multimesh( + self.data.mesh ): pid = self.data.mesh.elements.ravel()[pid] except Exception: @@ -3762,9 +3761,8 @@ def _picked(point): pid = mainwin.data.mesh.nearest_node(point) else: pid = mesh.find_closest_point(point) - if ( - mainwin.current_data_type == "GaussPoint" - and not _is_multimesh(mainwin.data.mesh) + if mainwin.current_data_type == "GaussPoint" and not _is_multimesh( + mainwin.data.mesh ): pid = mainwin.data.mesh.elements.ravel()[pid] @@ -3791,9 +3789,7 @@ def _on_click(obj, event): if cell_id < 0: return display_cell_id = cell_id - actor_name, mesh = mainwin.active_dock.picked_actor_mesh( - picker.GetActor() - ) + actor_name, mesh = mainwin.active_dock.picked_actor_mesh(picker.GetActor()) if mesh is None: mesh = mainwin.active_dock.pv_mesh cell = mesh.extract_cells(cell_id) diff --git a/fedoo/util/xdmf_writer.py b/fedoo/util/xdmf_writer.py index 3c339292..90bdd890 100644 --- a/fedoo/util/xdmf_writer.py +++ b/fedoo/util/xdmf_writer.py @@ -1,4 +1,3 @@ - from __future__ import annotations from pathlib import Path @@ -78,9 +77,7 @@ def _list_iterations(self, f: h5py.File) -> Iterable[int]: if grp is None: return [] return sorted( - int(name.split("_")[1]) - for name in grp - if name.startswith("iter_") + int(name.split("_")[1]) for name in grp if name.startswith("iter_") ) def _build_iteration_grid(self, f: h5py.File, iteration: int) -> etree.Element: From 61a4d2641a2d099bc8c9a1c6a2d92167f6d9b5e8 Mon Sep 17 00:00:00 2001 From: pruliere Date: Tue, 7 Jul 2026 11:06:27 +0200 Subject: [PATCH 03/14] update dependencies: add h5py and update minimal simcoon version --- conda.recipe/meta.yaml | 3 ++- environment.yml | 3 ++- environment_doc.yml | 3 ++- pyproject.toml | 5 +++-- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml index 6096ea6e..8777e851 100644 --- a/conda.recipe/meta.yaml +++ b/conda.recipe/meta.yaml @@ -20,8 +20,9 @@ requirements: - python>=3.10 - numpy>=1.26 - scipy + - h5py - pyvista>=0.42 - - simcoon>=1.11.0 + - simcoon>=1.13.1 - {{ solver_backend }} test: diff --git a/environment.yml b/environment.yml index 16d9514d..ebeac1a7 100644 --- a/environment.yml +++ b/environment.yml @@ -8,13 +8,14 @@ dependencies: - python>=3.10 - numpy>=1.26 - scipy + - h5py - pypardiso # [x86 or x86_64] - scikit-umfpack # [(not x86) and (not x86_64)] - python-mumps # exercises _solver_mumps default path and the # _MumpsFactor reuse adapter, regardless of platform - vtk>=9.2 - pyvista>=0.39 - - simcoon>=1.11.0 + - simcoon>=1.13.1 - pytest - pytest-xdist - sphinx_rtd_theme diff --git a/environment_doc.yml b/environment_doc.yml index 3305d7d5..f211dcb7 100644 --- a/environment_doc.yml +++ b/environment_doc.yml @@ -8,12 +8,13 @@ dependencies: - python=3.10 - numpy>=1.24 - scipy + - h5py - pypardiso # [x86 or x86_64] # - scikit-umfpack # [(not x86) and (not x86_64)] - vtk>=9.2 - pypandoc - pyvista>=0.39 - - simcoon>=1.11.0 + - simcoon>=1.13.1 - pytest - sphinx - sphinx_rtd_theme diff --git a/pyproject.toml b/pyproject.toml index f6b05e31..828470a0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,8 @@ classifiers = [ dependencies = [ 'numpy>=2.00', 'scipy', - 'simcoon>=1.11.0', + 'h5py', + 'simcoon>=1.13.1', ] [project.optional-dependencies] @@ -42,7 +43,7 @@ plot = [ 'pyvista[io]', ] simcoon = [ - 'simcoon>=1.11.0' + 'simcoon>=1.13.1' ] test = [ 'pytest', From e7b8f939473f2177422b7a0b6fe917ee650583d5 Mon Sep 17 00:00:00 2001 From: pruliere Date: Tue, 7 Jul 2026 15:11:55 +0200 Subject: [PATCH 04/14] add missing file --- fedoo/core/multimeshdata.py | 218 ++++++++++++++++++++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 fedoo/core/multimeshdata.py diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py new file mode 100644 index 00000000..0649f84f --- /dev/null +++ b/fedoo/core/multimeshdata.py @@ -0,0 +1,218 @@ +"""Data containers for values attached to ``MultiMesh`` submeshes.""" + +from __future__ import annotations + +import numpy as np + +from fedoo.core.mesh import MultiMesh + + +class MultiMeshData: + """Array-like data container associated with a ``MultiMesh``. + + ``MultiMeshData`` stores one optional data block per ordered submesh of a + ``MultiMesh``. It is returned by :meth:`DataSet.get_data` for element and + Gauss point fields attached to a ``MultiMesh``. The object behaves like the + data block of its active submesh when converted to a NumPy array or indexed + directly, while still allowing explicit submesh selection. + + Submeshes may be selected by integer id, submesh name, element type, or a + list combining these selectors: + + >>> data = dataset.get_data("Stress", "XX") + >>> data.active # data for dataset.active_submesh + >>> data.submesh(1) # data for submesh id 1 + >>> data.submesh("tri3") # data for all tri3 submeshes + + Parameters + ---------- + mesh : MultiMesh + Mesh defining the ordered submeshes. + data : dict, sequence, MultiMeshData or array-like + Data to associate with submeshes. Dictionaries may use integer + submesh ids, submesh names, or unique element types as keys. A sequence + must have one entry per submesh. A plain array is attached only to the + active submesh. + active_submesh : int or str, default=0 + Submesh used when ``MultiMeshData`` is accessed as an array. + + Notes + ----- + Missing submesh data are stored as ``None``. ``to_global`` concatenates + selected submesh blocks in mesh order and fills missing blocks with + ``fill_value``. + """ + + def __init__(self, mesh: MultiMesh, data, active_submesh=0) -> None: + self.mesh = mesh + self.active_submesh = active_submesh + self._data = self._normalize_data(data) + + def _normalize_data(self, data) -> tuple: + n_submesh = len(self.mesh.submeshes) + if isinstance(data, MultiMeshData): + return data._data + if isinstance(data, dict): + return tuple(self._value_from_dict(data, i) for i in range(n_submesh)) + if isinstance(data, (list, tuple)) and len(data) == n_submesh: + return tuple(data) + return tuple( + data if i == self._resolve_single(self.active_submesh) else None + for i in range(n_submesh) + ) + + def _value_from_dict(self, data: dict, index: int): + mesh = self.mesh[index] + if index in data: + return data[index] + if mesh.name and mesh.name in data: + return data[mesh.name] + matches = [ + i + for i, submesh in enumerate(self.mesh.submeshes) + if submesh.elm_type == mesh.elm_type + ] + if len(matches) == 1 and mesh.elm_type in data: + return data[mesh.elm_type] + return None + + def _resolve_indices(self, selector=None) -> list[int]: + if selector is None: + return [i for i, data in enumerate(self._data) if data is not None] + if isinstance(selector, (list, tuple, set, np.ndarray)): + indices = [] + for item in selector: + indices.extend(self._resolve_indices(item)) + return list(dict.fromkeys(indices)) + if isinstance(selector, int): + return [selector] + if isinstance(selector, str): + name_matches = [ + i for i, mesh in enumerate(self.mesh.submeshes) if mesh.name == selector + ] + if name_matches: + return name_matches + type_matches = [ + i + for i, mesh in enumerate(self.mesh.submeshes) + if mesh.elm_type == selector + ] + if type_matches: + return type_matches + raise KeyError(selector) + + def _resolve_single(self, selector=None) -> int: + indices = self._resolve_indices(selector) + if len(indices) != 1: + raise KeyError( + f"Selector {selector!r} matches {len(indices)} submeshes, expected one." + ) + return indices[0] + + @property + def active(self): + return self.submesh(self.active_submesh) + + def submesh(self, selector): + """Return data associated with one submesh or a filtered view. + + If ``selector`` resolves to a single submesh, the raw data block is + returned. If it resolves to several submeshes, a new ``MultiMeshData`` + view backed by a filtered ``MultiMesh`` is returned. + """ + indices = self._resolve_indices(selector) + if len(indices) == 1: + return self._data[indices[0]] + mesh = MultiMesh.from_mesh_list( + [self.mesh[i] for i in indices], + name=self.mesh.name, + node_sets=self.mesh.node_sets, + register_name=False, + ) + return MultiMeshData( + mesh, + { + j: self._data[i] + for j, i in enumerate(indices) + if self._data[i] is not None + }, + active_submesh=0, + ) + + def map(self, func) -> "MultiMeshData": + return MultiMeshData( + self.mesh, + {i: func(data) for i, data in enumerate(self._data) if data is not None}, + active_submesh=self.active_submesh, + ) + + def keys(self): + return [i for i, data in enumerate(self._data) if data is not None] + + def items(self): + return [(i, data) for i, data in enumerate(self._data) if data is not None] + + def to_global(self, indices=None, fill_value=np.nan): + """Concatenate per-submesh data in the selected submesh order.""" + if indices is None: + indices = list(range(len(self.mesh.submeshes))) + else: + indices = self._resolve_indices(indices) + template = next( + (self._data[i] for i in indices if self._data[i] is not None), + None, + ) + arrays = [] + for i in indices: + data = self._data[i] + n_elements = self.mesh[i].n_elements + if data is None: + if template is not None and np.asarray(template).ndim > 1: + shape = np.asarray(template).shape[:-1] + (n_elements,) + else: + shape = (n_elements,) + arrays.append(np.full(shape, fill_value)) + else: + arrays.append(np.asarray(data)) + if not arrays: + return np.array([]) + if arrays[0].ndim == 1: + return np.concatenate(arrays) + return np.concatenate(arrays, axis=-1) + + def __array__(self, dtype=None): + data = np.asarray(self.active) + if dtype is not None: + data = data.astype(dtype) + return data + + def __getitem__(self, item): + return self.active[item] + + @property + def shape(self): + return np.asarray(self.active).shape + + @property + def ndim(self): + return np.asarray(self.active).ndim + + +def copy_data_value(value, deep: bool = False): + """Copy a dataset value, including per-submesh dictionaries.""" + if isinstance(value, MultiMeshData): + value = value._data + + if isinstance(value, dict): + return { + key: copy_data_value(item, deep=deep) + for key, item in value.items() + } + + if isinstance(value, (list, tuple)): + copied = [copy_data_value(item, deep=deep) for item in value] + return type(value)(copied) + + if deep and hasattr(value, "copy"): + return value.copy() + return value From 92b181bb8616bf9e1b5bec2d33b7b3290d1424df Mon Sep 17 00:00:00 2001 From: pruliere Date: Wed, 8 Jul 2026 15:25:11 +0200 Subject: [PATCH 05/14] - update viewer (bugs corrections, add a save menu) - add methods to convert local to global id in multi_mesh data - update dependencies (matlplotlib<3.5 no more supported) --- fedoo/core/dataset.py | 107 +++++++++++-- fedoo/core/multimeshdata.py | 79 +++++++++- fedoo/util/viewer.py | 298 ++++++++++++++++++++++++++++++------ pyproject.toml | 5 +- 4 files changed, 414 insertions(+), 75 deletions(-) diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 2a040b93..66120465 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -210,6 +210,89 @@ def _resolve_submesh_indices(self, selector=None) -> list[int]: def _active_submesh_index(self) -> int: return self._resolve_submesh_indices(self.active_submesh)[0] + def global_element_location(self, element_id: int) -> tuple[int, int]: + """Return ``(submesh_id, local_element_id)`` for a global element id. + + For regular ``Mesh`` datasets, the returned submesh id is always 0. + For ``MultiMesh`` datasets, global element ids follow the same + concatenated order used by ``to_pyvista`` and multimesh plotting: all + elements of submesh 0, then all elements of submesh 1, and so on. + """ + element_id = int(element_id) + if not self._is_multimesh(): + if element_id < 0 or element_id >= self.mesh.n_elements: + raise IndexError(element_id) + return 0, element_id + + offset = 0 + for submesh_id, submesh in enumerate(self.mesh.submeshes): + stop = offset + submesh.n_elements + if offset <= element_id < stop: + return submesh_id, element_id - offset + offset = stop + raise IndexError(element_id) + + def global_element_id(self, submesh_id: int, local_element_id: int) -> int: + """Return the global element id for a submesh-local element id.""" + submesh_id = int(submesh_id) + local_element_id = int(local_element_id) + + if not self._is_multimesh(): + if submesh_id != 0: + raise IndexError(submesh_id) + if local_element_id < 0 or local_element_id >= self.mesh.n_elements: + raise IndexError(local_element_id) + return local_element_id + + if submesh_id < 0 or submesh_id >= len(self.mesh.submeshes): + raise IndexError(submesh_id) + submesh = self.mesh[submesh_id] + if local_element_id < 0 or local_element_id >= submesh.n_elements: + raise IndexError(local_element_id) + return ( + sum(mesh.n_elements for mesh in self.mesh.submeshes[:submesh_id]) + + local_element_id + ) + + def split_global_element_indices(self, element_ids) -> dict[int, np.ndarray]: + """Split global element ids into submesh-local element ids. + + Parameters + ---------- + element_ids : array-like of int + Element ids in the concatenated global numbering used by + multimesh plotting and viewer selections. + + Returns + ------- + dict[int, numpy.ndarray] + Mapping ``submesh_id -> local_element_ids``. Empty submeshes are + omitted. + """ + element_ids = np.asarray(element_ids, dtype=int) + if element_ids.ndim == 0: + element_ids = element_ids.reshape(1) + + if not self._is_multimesh(): + if np.any((element_ids < 0) | (element_ids >= self.mesh.n_elements)): + raise IndexError(element_ids) + return {0: element_ids} + + split_indices = {} + offset = 0 + for submesh_id, submesh in enumerate(self.mesh.submeshes): + stop = offset + submesh.n_elements + local_ids = ( + element_ids[(element_ids >= offset) & (element_ids < stop)] - offset + ) + if len(local_ids): + split_indices[submesh_id] = local_ids + offset = stop + + if len(split_indices) == 0 and len(element_ids) > 0: + raise IndexError(element_ids) + return split_indices + def _selected_multimesh(self, selected_submeshes=None): if not self._is_multimesh(): return self.mesh, [0] @@ -992,16 +1075,6 @@ def _plot_multimesh( ): submesh_indices = [self._active_submesh_index()] - resolved_data_type = None - if field is not None: - _, resolved_data_type = self.get_data( - field, - component, - data_type, - True, - fill_unused_nodes=None, - ) - offsets = {} offset = 0 for i, submesh in enumerate(self.mesh.submeshes): @@ -1033,12 +1106,11 @@ def _plot_multimesh( continue subdataset = self._submesh_dataset(submesh_id) - if ( - field is not None - and resolved_data_type == "GaussPoint" - and field not in subdataset.gausspoint_data - ): - continue + if field is not None: + try: + subdataset.get_data(field, component, data_type) + except NameError: + continue sub_kargs = dict(kargs) if base_name is not None: sub_kargs["name"] = f"{base_name}_{submesh_id}" @@ -1095,6 +1167,9 @@ def _plot_multimesh( pl = pv.Plotter(window_size=window_size) backgroundplotter = False + if field is not None and not plotted: + raise NameError(f"Field data {field!r} not found on any selected submesh.") + if title is None: title = "" if field is None else f"{field}_{component}" pl.add_text(title, name="name", color="Black", font_size=title_size) diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py index 0649f84f..c4802920 100644 --- a/fedoo/core/multimeshdata.py +++ b/fedoo/core/multimeshdata.py @@ -152,6 +152,80 @@ def keys(self): def items(self): return [(i, data) for i, data in enumerate(self._data) if data is not None] + def global_element_location(self, element_id: int) -> tuple[int, int]: + """Return ``(submesh_id, local_element_id)`` for a global element id. + + Global element ids follow the same concatenated order as + :meth:`to_global`: all elements of submesh 0, then all elements of + submesh 1, and so on. + """ + element_id = int(element_id) + offset = 0 + for submesh_id, submesh in enumerate(self.mesh.submeshes): + stop = offset + submesh.n_elements + if offset <= element_id < stop: + return submesh_id, element_id - offset + offset = stop + raise IndexError(element_id) + + def global_element_id(self, submesh_id: int, local_element_id: int) -> int: + """Return the global element id for a submesh-local element id.""" + submesh_id = int(submesh_id) + local_element_id = int(local_element_id) + if submesh_id < 0 or submesh_id >= len(self.mesh.submeshes): + raise IndexError(submesh_id) + submesh = self.mesh[submesh_id] + if local_element_id < 0 or local_element_id >= submesh.n_elements: + raise IndexError(local_element_id) + return ( + sum(mesh.n_elements for mesh in self.mesh.submeshes[:submesh_id]) + + local_element_id + ) + + def global_element_value(self, element_id: int): + """Return the value attached to one global element id. + + The value is read from the corresponding submesh block. For vector or + tensor element fields stored with elements on the last axis, the + returned value is ``block[..., local_element_id]``. + """ + submesh_id, local_id = self.global_element_location(element_id) + block = self._data[submesh_id] + if block is None: + return None + return np.asarray(block)[..., local_id] + + def global_element_values(self, element_ids, fill_value=np.nan): + """Return values attached to several global element ids. + + Missing submesh blocks are filled with ``fill_value``. The returned + array preserves the order of ``element_ids``. + """ + element_ids = np.asarray(element_ids, dtype=int) + scalar_input = element_ids.ndim == 0 + if scalar_input: + return self.global_element_value(int(element_ids)) + + values = [] + template = next((data for data in self._data if data is not None), None) + for element_id in element_ids: + submesh_id, local_id = self.global_element_location(int(element_id)) + block = self._data[submesh_id] + if block is None: + if template is None or np.asarray(template).ndim == 1: + values.append(fill_value) + else: + values.append( + np.full(np.asarray(template).shape[:-1], fill_value) + ) + else: + values.append(np.asarray(block)[..., local_id]) + if not values: + return np.array([]) + if np.asarray(values[0]).ndim == 0: + return np.asarray(values) + return np.stack(values, axis=-1) + def to_global(self, indices=None, fill_value=np.nan): """Concatenate per-submesh data in the selected submesh order.""" if indices is None: @@ -204,10 +278,7 @@ def copy_data_value(value, deep: bool = False): value = value._data if isinstance(value, dict): - return { - key: copy_data_value(item, deep=deep) - for key, item in value.items() - } + return {key: copy_data_value(item, deep=deep) for key, item in value.items()} if isinstance(value, (list, tuple)): copied = [copy_data_value(item, deep=deep) for item in value] diff --git a/fedoo/util/viewer.py b/fedoo/util/viewer.py index fc7adb92..5f735e6b 100644 --- a/fedoo/util/viewer.py +++ b/fedoo/util/viewer.py @@ -42,39 +42,30 @@ def _global_data_array(data): return data -def _global_element_nodes(mesh, element_id): - """Return node ids for a global element id.""" - if not _is_multimesh(mesh): - return mesh.elements[element_id] +def _data_n_iter(data): + return getattr(data, "n_iter", 1) - submesh_id, local_id = _global_element_location(mesh, element_id) - return mesh[submesh_id].elements[local_id] - -def _global_element_location(mesh, element_id): - """Return ``(submesh_id, local_element_id)`` for a global element id.""" - if not _is_multimesh(mesh): - return 0, element_id - - offset = 0 - for submesh_id, submesh in enumerate(mesh.submeshes): - stop = offset + submesh.n_elements - if offset <= element_id < stop: - return submesh_id, element_id - offset - offset = stop - raise IndexError(element_id) +def _as_3d_tuple(values, fill=0.0): + if values is None: + return (fill, fill, fill) + values = tuple(values) + if len(values) >= 3: + return values[:3] + return values + (fill,) * (3 - len(values)) def _element_data_value(data, field, component, data_type, element_id): values = data.get_data(field, component, data_type) - values = _global_data_array(values) + if isinstance(values, MultiMeshData): + return values.global_element_value(element_id) return values[..., element_id] def _gausspoint_values_for_element(data, field, component, element_id): values = data.get_data(field, component, "GaussPoint") if isinstance(values, MultiMeshData): - submesh_id, local_id = _global_element_location(data.mesh, element_id) + submesh_id, local_id = data.global_element_location(element_id) block = values.submesh(submesh_id) if block is None: return np.array([]) @@ -89,7 +80,7 @@ def _gausspoint_values_for_element(data, field, component, element_id): def _gausspoint_global_indices_for_element(data, field, component, element_id): values = data.get_data(field, component, "GaussPoint") if isinstance(values, MultiMeshData): - submesh_id, local_id = _global_element_location(data.mesh, element_id) + submesh_id, local_id = data.global_element_location(element_id) offset = 0 for i in range(submesh_id): block = values.submesh(i) @@ -181,14 +172,16 @@ def _applyStyle(self): class PlotDock(QDockWidget): _n_created_dock = 0 # total dock created - def __init__(self, data, title, parent=None, opts=None): + def __init__(self, data, title, parent=None, opts=None, filename=None): PlotDock._n_created_dock += 1 self._dock_index = PlotDock._n_created_dock self.title = title + self.filename = filename title = f"{self._dock_index }: " + str(title) super().__init__(title, parent) self.data = data self.element_sets = None + self.unsaved_element_set_names = set() # visibility_mode to keep track on current_selection for element_sets dialog box # -> None for show all, "show" or "hide" self.visibility_mode = None @@ -458,6 +451,7 @@ def __init__(self, data=None): self.element_sets_dialog = None self._window_index = 1 self._picking_target = -1 # internal arg for picking tools. -1 = No target + self._syncing_line_widget = False self.window_layout = self._distribute_auto # if plane_widget or line wiget should be shown in the active dock @@ -671,9 +665,21 @@ def __init__(self, data=None): open_action.triggered.connect(self.open_file) file_menu.addAction(open_action) + # Save Data + save_action = QtWidgets.QAction("Save", self) + save_action.setShortcut("Ctrl+S") + save_action.triggered.connect(self.save_current_data) + file_menu.addAction(save_action) + + save_as_action = QtWidgets.QAction("Save as...", self) + save_as_action.setShortcut("Ctrl+Shift+S") + save_as_action.triggered.connect(self.save_current_data_as) + file_menu.addAction(save_as_action) + + file_menu.addSeparator() + # Save Image save_image_action = QtWidgets.QAction("Save Image...", self) - save_image_action.setShortcut("Ctrl+S") save_image_action.triggered.connect(self.save_image_dialog) file_menu.addAction(save_image_action) @@ -740,6 +746,7 @@ def __init__(self, data=None): clipAction.setShortcut("Ctrl+K") clipAction.triggered.connect(self.open_clip_dialog) tools_menu.addAction(clipAction) + self.clip_action = clipAction plotOverLineAction = QtWidgets.QAction("Plot over line...", self) plotOverLineAction.triggered.connect(self.open_plot_over_line_dialog) @@ -840,6 +847,8 @@ def __init__(self, data=None): clim_action, renderer_options_action, copy_action, + save_action, + save_as_action, clipAction, plotOverLineAction, element_sets_action, @@ -857,13 +866,15 @@ def __init__(self, data=None): # Dockable PyVista Widget # ------------------------- if data is not None: + filename = None if isinstance(data, str): + filename = data title = os.path.basename(data) data = fd.read_data(data) else: title = "Data" - self.add_dataset_dock(data, title=title) + self.add_dataset_dock(data, title=title, filename=filename) # self.plotter = QtInteractor(self) # dock_plotter = QDockWidget("3D View", self) # dock_plotter.setWidget(self.plotter.interactor) @@ -895,8 +906,8 @@ def update_plot(self, *args, **kargs): if self.active_dock: return self.active_dock.update_plot(*args, **kargs) - def add_dataset_dock(self, data, title, opts=None): - dock = PlotDock(data, title, self, opts) + def add_dataset_dock(self, data, title, opts=None, filename=None): + dock = PlotDock(data, title, self, opts, filename=filename) self.addDockWidget(Qt.RightDockWidgetArea, dock) dock.visibilityChanged.connect( @@ -1051,6 +1062,12 @@ def _distribute_auto(self): dock_index += 1 self.window_layout = self._distribute_auto + def _clip_is_supported(self, dock=None): + if dock is None: + dock = self.active_dock + mesh = getattr(getattr(dock, "data", None), "mesh", None) + return mesh is not None and getattr(mesh, "ndim", 3) == 3 + def _toggle_link_views(self, checked): # Synchronize both menu items self._link_views_action.setChecked(checked) @@ -1096,7 +1113,10 @@ def _toggle_sync_iter(self, checked): # Synchronize all docks to active dock's current iteration active_iter = self.active_dock.current_iter for dock in self.all_docks: - if dock is not self.active_dock and dock.data.n_iter > active_iter: + if ( + dock is not self.active_dock + and _data_n_iter(dock.data) > active_iter + ): dock.current_iter = active_iter dock.update_plot(iteration=active_iter) @@ -1159,13 +1179,17 @@ def _set_active(self, dock): self.active_dock = dock + if hasattr(self, "clip_action"): + clip_supported = self._clip_is_supported(dock) + self.clip_action.setEnabled(clip_supported) + if not clip_supported: + dock.opts["clip_args"] = None + self.disable_plane_widget() + for d in self.all_docks: d._titlebar._set_active(d is dock) - if hasattr(dock.data, "n_iter"): - max_iter = dock.data.n_iter - 1 - else: - max_iter = 0 + max_iter = _data_n_iter(dock.data) - 1 current_iter = dock.current_iter # may be modified by setRange self.setRange(0, max_iter) self.iter_slider.setValue(current_iter) @@ -1256,7 +1280,9 @@ def _on_slider_changed(self, val: int): # Sync iteration across all docks if enabled if self._sync_iter: for dock in self.all_docks: - if dock.data.n_iter > val: # only update if the dock has this iteration + if ( + _data_n_iter(dock.data) > val + ): # only update if the dock has this iteration dock.current_iter = val dock.update_plot(iteration=val) else: @@ -1275,7 +1301,9 @@ def _on_spin_changed(self, val: int): # Sync iteration across all docks if enabled if self._sync_iter: for dock in self.all_docks: - if dock.data.n_iter > val: # only update if the dock has this iteration + if ( + _data_n_iter(dock.data) > val + ): # only update if the dock has this iteration dock.current_iter = val dock.update_plot(iteration=val) else: @@ -1385,7 +1413,133 @@ def open_file(self): if fname: data = fd.read_data(fname) title = os.path.basename(fname) - self.add_dataset_dock(data, title=title) + self.add_dataset_dock(data, title=title, filename=fname) + + def _fdh5_save_filename(self, filename): + root, ext = os.path.splitext(filename) + if ext == "": + return root + ".fdh5" + if ext.lower() != ".fdh5": + return None + return filename + + def _default_element_sets_for_dock(self, dock): + element_sets = {"All": None} + unsaved = {"All"} + + mesh = dock.data.mesh + if _is_multimesh(mesh): + offset = 0 + for submesh_id, submesh in enumerate(mesh.submeshes): + label = submesh.name or submesh.elm_type + name = f"submesh_{submesh_id}: {label}" + if submesh.elm_type not in name: + name = f"{name} ({submesh.elm_type})" + base_name = name + counter = 1 + while name in element_sets or name in mesh.element_sets: + counter += 1 + name = f"{base_name} #{counter}" + element_sets[name] = list(range(offset, offset + submesh.n_elements)) + unsaved.add(name) + offset += submesh.n_elements + + element_sets.update(mesh.element_sets) + return element_sets, unsaved + + def _ensure_dock_element_sets(self, dock): + if dock.element_sets is None: + element_sets, unsaved = self._default_element_sets_for_dock(dock) + dock.element_sets = element_sets + dock.unsaved_element_set_names = unsaved + + def _element_sets_for_save(self, dock): + self._ensure_dock_element_sets(dock) + excluded = getattr(dock, "unsaved_element_set_names", set()) + return { + name: np.asarray(values, dtype=int) + for name, values in dock.element_sets.items() + if name not in excluded and values is not None + } + + def _apply_saved_element_sets_to_mesh(self, mesh, element_sets): + if not _is_multimesh(mesh): + mesh.element_sets = { + name: np.asarray(values, dtype=int) + for name, values in element_sets.items() + } + return + + submesh_sets = [{} for _ in mesh.submeshes] + dataset = fd.DataSet(mesh) + for name, values in element_sets.items(): + for submesh_id, local in dataset.split_global_element_indices( + values + ).items(): + submesh_sets[submesh_id][name] = local + + for submesh, sets in zip(mesh.submeshes, submesh_sets): + submesh.element_sets = sets + + def _data_for_fdh5_save(self, dock): + data = dock.data.deepcopy() + self._apply_saved_element_sets_to_mesh( + data.mesh, + self._element_sets_for_save(dock), + ) + return data + + def _save_data_to_fdh5(self, filename): + if self.active_dock is None: + return False + + filename = self._fdh5_save_filename(filename) + if filename is None: + QtWidgets.QMessageBox.warning( + self, + "Save data", + "Only FDH5 files (*.fdh5) are supported.", + ) + return False + + data = self._data_for_fdh5_save(self.active_dock) + if isinstance(data, fd.MultiFrameDataSet): + data.save_all(filename, file_format="fdh5") + else: + data.save(filename) + + self.active_dock.filename = filename + self.active_dock.title = os.path.basename(filename) + self.active_dock.setWindowTitle( + f"{self.active_dock._dock_index}: {self.active_dock.title}" + ) + QtWidgets.QMessageBox.information(self, "Save data", f"Saved:\n{filename}") + return True + + def save_current_data(self): + if self.active_dock is None: + return + filename = self.active_dock.filename + if not filename or os.path.splitext(filename)[1].lower() != ".fdh5": + self.save_current_data_as() + return + self._save_data_to_fdh5(filename) + + def save_current_data_as(self): + if self.active_dock is None: + return + start = self.active_dock.filename or self.active_dock.title + if not start or os.path.splitext(start)[1].lower() != ".fdh5": + start = os.path.splitext(start)[0] + ".fdh5" + + filename, _ = QtWidgets.QFileDialog.getSaveFileName( + self, + "Save data as", + start, + "Fedoo HDF5 files (*.fdh5)", + ) + if filename: + self._save_data_to_fdh5(filename) def open_clim_dialog(self): if self._clim_dialog is None: @@ -1761,7 +1915,7 @@ def on_use_current_toggled(checked: bool): ) pl.write_frame() - for iteration in range(1, self.data.n_iter): + for iteration in range(1, _data_n_iter(self.data)): rot_azimuth = float(rot_azimuth_spin.value()) rot_elevation = float(rot_elevation_spin.value()) if rot_azimuth != 0: @@ -1809,10 +1963,14 @@ def reset_camera(self): # Clip dialog : Open and close # ---------------------------- def open_clip_dialog(self): + if not self._clip_is_supported(): + self.statusBar().showMessage("Clip is only available for 3D meshes.", 3000) + return + # Create clip dialog only if not already exist if self._clip_dialog is None: origin = ( - tuple(self.data.mesh.bounding_box.center) + _as_3d_tuple(self.data.mesh.bounding_box.center) if self.data.mesh is not None else (0.0, 0.0, 0.0) ) @@ -1875,7 +2033,7 @@ def enable_plane_widget(self): bounds = tuple(np.array(self.data.mesh.as_3d().bounding_box).T.ravel()) if self.opts["clip_origin"] is None: - self.opts["clip_origin"] = tuple(self.data.mesh.bounding_box.center) + self.opts["clip_origin"] = _as_3d_tuple(self.data.mesh.bounding_box.center) self.opts["clip_normal"] = (1.0, 0.0, 0.0) origin = self.opts["clip_origin"] @@ -1960,14 +2118,31 @@ def _on_pol_dialog_closed(self, *args, **kwargs): self.plotter.clear_line_widgets() def _on_pol_dialog_changed(self, p1, p2): - # update widget (emit no signal) - self._line_widget.SetPoint1(p1) - self._line_widget.SetPoint2(p2) + self._set_line_widget_points(p1, p2) def on_line_changed(self, p1, p2): """Called interactively while moving the widget""" + if self._syncing_line_widget: + return + p1, p2 = self._set_line_widget_points(p1, p2) self._plot_over_line_dialog.update_line(p1, p2) + def _line_points_for_active_mesh(self, p1, p2): + p1, p2 = tuple(p1), tuple(p2) + if self.data.mesh.ndim == 2: + p1, p2 = (p1[0], p1[1], 0.0), (p2[0], p2[1], 0.0) + return p1, p2 + + def _set_line_widget_points(self, p1, p2): + p1, p2 = self._line_points_for_active_mesh(p1, p2) + self._syncing_line_widget = True + try: + self._line_widget.SetPoint1(p1) + self._line_widget.SetPoint2(p2) + finally: + self._syncing_line_widget = False + return p1, p2 + def _rebuild_line_widget(self): if self._plot_over_line_dialog: p1 = self._plot_over_line_dialog.p1 @@ -1998,6 +2173,7 @@ def _picked(point): p1 = tuple(point) else: p2 = tuple(point) + p1, p2 = self._line_points_for_active_mesh(p1, p2) # Sync UI and scene self._on_pol_dialog_changed(p1, p2) # update widget @@ -2548,7 +2724,7 @@ def closeEvent(self, event): def get_values(self): cmap_name = self.cmap_combo.currentText() n_colors = int(self.ncolors_spin.value()) - cmap = mpl.cm.get_cmap(cmap_name, n_colors) + cmap = mpl.colormaps[cmap_name].resampled(n_colors) if self.cmap_reverse.isChecked(): cmap = cmap.reversed() if self.rb_current.isChecked(): @@ -2583,7 +2759,7 @@ def __init__(self, height=24, parent=None): def update_preview(self, cmap_name: str, reverse: bool = False, n_colors=256): try: # cmap = mpl.colormaps.get(cmap_name) - cmap = mpl.cm.get_cmap(cmap_name, n_colors) + cmap = mpl.colormaps[cmap_name].resampled(n_colors) except Exception: self.setText(f"Colormap '{cmap_name}' introuvable") return @@ -2737,6 +2913,8 @@ def __init__( super().__init__(parent) self.setWindowTitle("Clip Plane") self.setModal(False) # non modale + default_origin = _as_3d_tuple(default_origin) + default_normal = _as_3d_tuple(default_normal) # --- Widgets --------------------------------------------------------- self.enableClipChk = QtWidgets.QCheckBox("Activate clipping") @@ -2840,6 +3018,8 @@ def _emit_clip_params(self): def set_values_from_widget(self, origin, normal, invert=None): """Update widget values without emiting signal.""" + origin = _as_3d_tuple(origin) + normal = _as_3d_tuple(normal) blockers = [ QSignalBlocker(self.originX), QSignalBlocker(self.originY), @@ -3196,20 +3376,21 @@ def update_plot(self): self.parent().update_plot() # -------------- Element sets UI helpers ------------- + def _default_element_sets(self): + return self.parent()._default_element_sets_for_dock(self.parent().active_dock) + def _refresh_sets_list(self): if self.element_sets is None: - element_sets = { - "All": None, # or set(range(self.mesh.n_cells)) - # add your pre-defined sets... - } - element_sets.update(self.parent().data.mesh.element_sets) + element_sets, unsaved = self._default_element_sets() self.parent().active_dock.element_sets = element_sets + self.parent().active_dock.unsaved_element_set_names = unsaved self.sets_list.clear() for name in sorted(self.element_sets.keys()): - if self.element_sets[name]: - n_elements = len(self.element_sets[name]) - else: + values = self.element_sets[name] + if values is None: n_elements = self.parent().data.mesh.n_elements + else: + n_elements = len(values) item = QtWidgets.QListWidgetItem(f"{name} ({n_elements} elements)") item.setData(Qt.UserRole, name) self.sets_list.addItem(item) @@ -3419,6 +3600,11 @@ def _rename_selected_set(self): self, "No set selected", "Please select a set to rename." ) return + if name in self.parent().active_dock.unsaved_element_set_names: + QtWidgets.QMessageBox.information( + self, "Built-in set", "Built-in element sets cannot be renamed." + ) + return new_name, ok = QtWidgets.QInputDialog.getText( self, "Rename Set", f"New name for '{name}':", text=name ) @@ -3441,6 +3627,11 @@ def _remove_selected_set(self): self, "No set selected", "Please select a set to rename." ) return + if name in self.parent().active_dock.unsaved_element_set_names: + QtWidgets.QMessageBox.information( + self, "Built-in set", "Built-in element sets cannot be removed." + ) + return del self.element_sets[name] self._refresh_sets_list() self.status_lbl.setText(f"Set '{name}' removed") @@ -3457,6 +3648,11 @@ def _save_selection_as_new_set(self): self, "Name exists", f"A set named '{new_name}' already exists." ) return + if new_name in self.parent().active_dock.unsaved_element_set_names: + QtWidgets.QMessageBox.warning( + self, "Reserved name", f"'{new_name}' is a built-in set name." + ) + return if self.current_selection: self.element_sets[new_name] = list(self.current_selection) else: diff --git a/pyproject.toml b/pyproject.toml index 828470a0..aa218177 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,12 +39,9 @@ solver = [ 'python-mumps ; platform_machine=="arm64" or platform_machine=="aarch64"', ] plot = [ - 'matplotlib', + 'matplotlib>=3.5', 'pyvista[io]', ] -simcoon = [ - 'simcoon>=1.13.1' -] test = [ 'pytest', 'pytest-cov' From e26d5c2ad0177b8dcb221ca0276e28a950217f49 Mon Sep 17 00:00:00 2001 From: pruliere Date: Wed, 8 Jul 2026 17:41:57 +0200 Subject: [PATCH 06/14] allow to steck mesh with different element types -> MultiMesh --- .../01-simple/spherical_shell_compression.py | 2 +- fedoo/__init__.py | 1 + fedoo/core/__init__.py | 2 +- fedoo/core/mesh.py | 67 +++++++++++++------ 4 files changed, 50 insertions(+), 22 deletions(-) diff --git a/examples/01-simple/spherical_shell_compression.py b/examples/01-simple/spherical_shell_compression.py index 8799d2ce..2bad562b 100644 --- a/examples/01-simple/spherical_shell_compression.py +++ b/examples/01-simple/spherical_shell_compression.py @@ -21,7 +21,7 @@ ############################################################################### # Create a simple sphere mesh using pyvista. -mesh = fd.Mesh.from_pyvista(pv.Sphere(radius)) +mesh = fd.Mesh.from_pyvista(pv.Sphere(radius=radius)) ############################################################################### # Define a linear isotropic material and an homogeneous shell section diff --git a/fedoo/__init__.py b/fedoo/__init__.py index 1177adcb..937fd736 100644 --- a/fedoo/__init__.py +++ b/fedoo/__init__.py @@ -18,6 +18,7 @@ DataSet, ListBC, Mesh, + MultiMesh, ModelingSpace, MultiFrameDataSet, MultiMeshData, diff --git a/fedoo/core/__init__.py b/fedoo/core/__init__.py index 5223bd78..3c1f4a32 100644 --- a/fedoo/core/__init__.py +++ b/fedoo/core/__init__.py @@ -1,6 +1,6 @@ from .modelingspace import ModelingSpace from .base import ConstitutiveLaw -from .mesh import Mesh +from .mesh import Mesh, MultiMesh from .weakform import WeakForm, WeakFormSum from .assembly import Assembly from .assembly_sum import AssemblySum diff --git a/fedoo/core/mesh.py b/fedoo/core/mesh.py index 06beb64e..165d6b03 100644 --- a/fedoo/core/mesh.py +++ b/fedoo/core/mesh.py @@ -395,56 +395,83 @@ def add_internal_nodes(self, nb_added: int) -> np.ndarray[int]: # warning , this method must be static @staticmethod - def stack(mesh1: "Mesh", mesh2: "Mesh", name: str = "") -> "Mesh": + def stack(mesh1: "Mesh", mesh2: "Mesh", name: str = "") -> "Mesh | MultiMesh": """Add two mesh together to make a new mesh. *Static method* - Make the spatial stack of two mesh objects which have the - same element shape. This function doesn't merge coindicent Nodes. + same element shape. If the two meshes use different element shapes, a + MultiMesh is returned. This function doesn't merge coincident nodes. For that purpose, use the Mesh methods 'find_coincident_nodes' and 'merge_nodes' on the resulting Mesh. Return --------- - Mesh object with is the spacial stack of mesh1 and mesh2 + Mesh or MultiMesh object which is the spatial stack of mesh1 and mesh2 """ if isinstance(mesh1, str): mesh1 = Mesh.get_all()[mesh1] if isinstance(mesh2, str): mesh2 = Mesh.get_all()[mesh2] - if mesh1.elm_type != mesh2.elm_type: - raise NameError("Can only stack meshes with the same element shape") - n_nodes = mesh1.n_nodes n_elements = mesh1.n_elements - new_crd = np.r_[mesh1.nodes, mesh2.nodes] - new_elm = np.r_[mesh1.elements, mesh2.elements + n_nodes] + ndim = max(mesh1.ndim, mesh2.ndim) + mesh1 = mesh1.as_ndim(ndim) + mesh2 = mesh2.as_ndim(ndim) - new_ndSets = dict(mesh1.node_sets) + nodes = np.r_[mesh1.nodes, mesh2.nodes] + node_sets = dict(mesh1.node_sets) for key in mesh2.node_sets: - if key in mesh1.node_sets: - new_ndSets[key] = np.r_[ + if key in node_sets: + node_sets[key] = np.r_[ mesh1.node_sets[key], np.array(mesh2.node_sets[key]) + n_nodes, ] else: - new_ndSets[key] = np.array(mesh2.node_sets[key]) + n_nodes + node_sets[key] = np.array(mesh2.node_sets[key]) + n_nodes + + if mesh1.elm_type != mesh2.elm_type: + return MultiMesh( + nodes, + { + mesh1.elm_type: ( + mesh1.elm_type, + mesh1.elements, + mesh1.element_sets, + ), + mesh2.elm_type: ( + mesh2.elm_type, + mesh2.elements + n_nodes, + mesh2.element_sets, + ), + }, + node_sets=node_sets, + ndim=ndim, + name=name, + ) - new_elSets = dict(mesh1.element_sets) + elements = np.r_[mesh1.elements, mesh2.elements + n_nodes] + + element_sets = dict(mesh1.element_sets) for key in mesh2.element_sets: - if key in mesh1.element_sets: - new_elSets[key] = np.r_[ + if key in element_sets: + element_sets[key] = np.r_[ mesh1.element_sets[key], np.array(mesh2.element_sets[key]) + n_elements, ] else: - new_elSets[key] = np.array(mesh2.element_sets[key]) + n_elements + element_sets[key] = np.array(mesh2.element_sets[key]) + n_elements - mesh3 = Mesh(new_crd, new_elm, mesh1.elm_type, name=name) - mesh3.node_sets = new_ndSets - mesh3.element_sets = new_elSets - return mesh3 + return Mesh( + nodes, + elements, + mesh1.elm_type, + node_sets=node_sets, + element_sets=element_sets, + ndim=ndim, + name=name, + ) def find_coincident_nodes(self, tol: float = 1e-8) -> np.ndarray[int]: """Find some nodes with the same position considering a tolerance given by the argument tol. From d3049fb1478730ad308a1eeaf61d19fb2b0ef9bf Mon Sep 17 00:00:00 2001 From: pruliere Date: Thu, 9 Jul 2026 14:30:17 +0200 Subject: [PATCH 07/14] correct few bugs and code cleaning --- fedoo/constitutivelaw/shell.py | 53 ++++++------- fedoo/constraint/distributed_load.py | 91 ++++++++++++++++------- fedoo/core/assembly.py | 53 ++++++++----- fedoo/core/assembly_sum.py | 14 ---- fedoo/core/dataset.py | 14 ++-- fedoo/core/multimeshdata.py | 23 ++++-- fedoo/core/output.py | 53 +++++++------ tests/test_distributed_load_neumann_bc.py | 65 +++++++++++++++- 8 files changed, 249 insertions(+), 117 deletions(-) diff --git a/fedoo/constitutivelaw/shell.py b/fedoo/constitutivelaw/shell.py index 542f70f8..694967f5 100644 --- a/fedoo/constitutivelaw/shell.py +++ b/fedoo/constitutivelaw/shell.py @@ -64,28 +64,28 @@ def update(self, assembly, pb): # rot = pb.get_rot() U = pb.get_dof_solution() if np.isscalar(U) and U == 0: - assembly.sv["GeneralizedStrain"] = 0 - assembly.sv["GeneralizedStress"] = 0 + assembly.sv["ShellStrain"] = 0 + assembly.sv["ShellStress"] = 0 else: - GeneralizedStrainOp = assembly.weakform.generalized_strain_operator() - GeneralizedStrain = [ + ShellStrainOp = assembly.weakform.generalized_strain_operator() + ShellStrain = [ op if np.isscalar(op) else assembly.get_gp_results(op, U) - for op in GeneralizedStrainOp + for op in ShellStrainOp ] H = self.get_shell_stiffness_matrix() # all terms are computed. Perhaps could be optimized by computed only the termes related to the associated weak form (eg shear for selective integration) - assembly.sv["GeneralizedStress"] = [ + assembly.sv["ShellStress"] = [ sum( [ - GeneralizedStrain[j] * assembly.convert_data(H[i][j]) + ShellStrain[j] * assembly.convert_data(H[i][j]) for j in range(8) ] ) for i in range(8) ] # H[i][j] are converted to gauss point excepted if scalar - assembly.sv["GeneralizedStrain"] = GeneralizedStrain + assembly.sv["ShellStrain"] = ShellStrain def get_strain(self, assembly, **kargs): """Return the last computed strain associated to the given assembly. @@ -110,15 +110,18 @@ def get_strain(self, assembly, **kargs): z = position * self.thickness / 2 Strain = StrainTensorList([0 for i in range(6)]) - GeneralizedStrain = assembly.sv["GeneralizedStrain"] + ShellStrain = assembly.sv["ShellStrain"] + if np.isscalar(ShellStrain) and ShellStrain == 0: + zeros = np.zeros(assembly.n_gauss_points) + return StrainTensorList([zeros.copy() for _ in range(6)]) Strain[0] = ( - GeneralizedStrain[0] + z * GeneralizedStrain[4] + ShellStrain[0] + z * ShellStrain[4] ) # epsXX -> membrane and bending Strain[1] = ( - GeneralizedStrain[1] - z * GeneralizedStrain[3] + ShellStrain[1] - z * ShellStrain[3] ) # epsYY -> membrane and bending - Strain[3] = GeneralizedStrain[2] # 2epsXY - Strain[4:6] = GeneralizedStrain[6:8] # 2epsXZ and 2epsYZ -> shear + Strain[3] = ShellStrain[2] # 2epsXY + Strain[4:6] = ShellStrain[6:8] # 2epsXZ and 2epsYZ -> shear return Strain @@ -204,16 +207,16 @@ def GetStressDistribution(self, assembly, pg, resolution=100): z = np.arange(-h / 2, h / 2, h / resolution) Strain = StrainTensorList([0 for i in range(6)]) - GeneralizedStrain = assembly.sv["GeneralizedStrain"] + ShellStrain = assembly.sv["ShellStrain"] Strain[0] = ( - GeneralizedStrain[0][pg] + z * GeneralizedStrain[4][pg] + ShellStrain[0][pg] + z * ShellStrain[4][pg] ) # epsXX -> membrane and bending Strain[1] = ( - GeneralizedStrain[1][pg] - z * GeneralizedStrain[3][pg] + ShellStrain[1][pg] - z * ShellStrain[3][pg] ) # epsYY -> membrane and bending - Strain[3] = GeneralizedStrain[2][pg] * np.ones_like(z) # 2epsXY - Strain[4] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear - Strain[5] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear + Strain[3] = ShellStrain[2][pg] * np.ones_like(z) # 2epsXY + Strain[4] = ShellStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear + Strain[5] = ShellStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear Hplane = self.material.get_elastic_matrix( "2Dstress" @@ -389,17 +392,17 @@ def GetStressDistribution(self, assembly, pg, resolution=100): z = np.linspace(-h / 2, h / 2, resolution) Strain = StrainTensorList([0 for i in range(6)]) - GeneralizedStrain = assembly.sv["GeneralizedStrain"] + ShellStrain = assembly.sv["ShellStrain"] Strain[0] = ( - GeneralizedStrain[0][pg] + z * GeneralizedStrain[4][pg] + ShellStrain[0][pg] + z * ShellStrain[4][pg] ) # epsXX -> membrane and bending Strain[1] = ( - GeneralizedStrain[1][pg] - z * GeneralizedStrain[3][pg] + ShellStrain[1][pg] - z * ShellStrain[3][pg] ) # epsYY -> membrane and bending - Strain[3] = GeneralizedStrain[2][pg] * np.ones_like(z) # 2epsXY - Strain[4] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear - Strain[5] = GeneralizedStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear + Strain[3] = ShellStrain[2][pg] * np.ones_like(z) # 2epsXY + Strain[4] = ShellStrain[6][pg] * np.ones_like(z) # 2epsXZ -> shear + Strain[5] = ShellStrain[6][pg] * np.ones_like(z) # 2epsYZ -> shear layer_z = [ list((pos - self.__layer) <= 0).index(True) - 1 for pos in z diff --git a/fedoo/constraint/distributed_load.py b/fedoo/constraint/distributed_load.py index b3fdda45..d887ee94 100644 --- a/fedoo/constraint/distributed_load.py +++ b/fedoo/constraint/distributed_load.py @@ -17,10 +17,9 @@ class _AssemblyNeumannBC(BCBase): """Neumann boundary condition generated from a load assembly. - The wrapped assembly computes the equivalent nodal load vector, but the - vector is injected through ``Problem.B`` like a standard Neumann boundary - condition. This keeps distributed loads out of the problem's internal - residual assembly during nonlinear solves. + The wrapped assembly computes equivalent nodal load reference vectors. + Fixed-geometry loads cache those vectors, then apply the time factor without + reassembling. Follower loads refresh the reference vectors when generated. """ def __init__(self, assembly: Assembly, name: str = "", time_func=None): @@ -28,15 +27,22 @@ def __init__(self, assembly: Assembly, name: str = "", time_func=None): self.bc_type = "Neumann" self.assembly = assembly self._start_value_default = 0 + self.start_value = self._start_value_default + self.value = 0 if time_func is None: def time_func(t_fact): return t_fact self.time_func = time_func - self._assembly_t_fact = None self._assembly_vector = None + def str_condensed(self): + """Return a condensed one line str describing the object.""" + if self.name == "": + return "Neumann -> assembly load" + return "Neumann (name = '{}') -> assembly load".format(self.name) + def initialize(self, problem: ProblemBase): if self.assembly.mesh.n_nodes != problem.mesh.n_nodes: raise ValueError( @@ -45,8 +51,58 @@ def initialize(self, problem: ProblemBase): self.assembly.initialize(problem) self._update_during_inc = bool(self.assembly._nlgeom) - self._assembly_t_fact = None self._assembly_vector = None + self.value = 0 + self.start_value = self._start_value_default + + def _get_factor(self, t_fact=1, t_fact_old=None): + return self.time_func(t_fact) + + def get_value(self, t_fact=1, t_fact_old=None): + factor = self._get_factor(t_fact, t_fact_old) + if factor == 0: + return self.start_value + elif self.start_value is None: + return factor * self.value + else: + return factor * (self.value - self.start_value) + self.start_value + + def get_true_value(self, t_fact=1, t_fact_old=None): + return self.get_value(t_fact, t_fact_old) + + def _format_load_vector(self, problem: ProblemBase, value): + if np.isscalar(value) and value == 0: + return 0 + if problem.n_global_dof and len(value) < problem.n_dof: + value = np.pad(value, (0, problem.n_dof - len(value))) + return value.copy() if hasattr(value, "copy") else value + + def _assemble_load_vector(self, problem: ProblemBase, load_factor): + self.assembly.set_start(problem, t_fact=load_factor) + if self._update_during_inc: + self.assembly.update(problem, compute="vector") + else: + self.assembly.assemble_global_mat(compute="vector") + return self._format_load_vector( + problem, + self.assembly.current.get_global_vector(), + ) + + def _has_initial_load(self): + return any( + getattr(self.assembly, attr, None) is not None + for attr in ("initial_pressure", "initial_force") + ) + + def _refresh_reference_values(self, problem: ProblemBase): + if self._has_initial_load(): + self.start_value = self._assemble_load_vector(problem, 0) + else: + self.start_value = 0 + self.value = self._assemble_load_vector(problem, 1) + self._assembly_vector = ( + self.value.copy() if hasattr(self.value, "copy") else self.value + ) def _format_current_value(self, problem: ProblemBase, value): if np.isscalar(value) and value == 0: @@ -54,33 +110,16 @@ def _format_current_value(self, problem: ProblemBase, value): self._current_value = np.array([]) return - if problem.n_global_dof and len(value) < problem.n_dof: - value = np.pad(value, (0, problem.n_dof - len(value))) self._dof_index = np.arange(problem.n_dof, dtype=int) self._current_value = value def generate(self, problem: ProblemBase, t_fact=1, t_fact_old=None): - assembly_t_fact = self.time_func(t_fact) - factor_changed = not np.array_equal(self._assembly_t_fact, assembly_t_fact) - - if factor_changed: - self.assembly.set_start(problem, t_fact=assembly_t_fact) - self._assembly_t_fact = assembly_t_fact + if self._update_during_inc or self._assembly_vector is None: + self._refresh_reference_values(problem) - if self._update_during_inc: - self.assembly.update(problem, compute="vector") - value = self.assembly.current.get_global_vector() - elif factor_changed or self._assembly_vector is None: - self.assembly.assemble_global_mat(compute="vector") - value = self.assembly.current.get_global_vector() - self._assembly_vector = value.copy() if hasattr(value, "copy") else value - else: - value = self._assembly_vector - - self._format_current_value(problem, value) + self._format_current_value(problem, self.get_value(t_fact, t_fact_old)) return [self] - class Pressure(Assembly): """Pressure load. diff --git a/fedoo/core/assembly.py b/fedoo/core/assembly.py index b9cc9369..bc05559a 100644 --- a/fedoo/core/assembly.py +++ b/fedoo/core/assembly.py @@ -1177,7 +1177,7 @@ def _get_assembled_operator(self, operator, n_elm_gp=None): return matrix - def operator_apply(self, wf, nodal_vector): + def operator_apply(self, wf, nodal_vector, use_local_dof=False): """Apply a discrete nodal vector to the trial space of a bilinear operator. This function performs a contraction of a bilinear form a(v, u*) with @@ -1194,14 +1194,39 @@ def operator_apply(self, wf, nodal_vector): nodal_vector : array_like The discrete values (Degrees of Freedom) used to instantiate - the trial field. + the trial field. By default, the vector is interpreted in the + global nodal basis. + use_local_dof : bool, default: False + If True, nodal_vector is interpreted as local element dofs, ordered + by variable blocks with ``n_elements * n_elm_nodes`` values each. Return: DiffOp : An operator containing only the virtual (test) operator, ready for vector assembly. """ - n_nodes = self.mesh.n_nodes new_wf = type(wf)([], [], []) + mat_change_of_basis = self.get_change_of_basis_mat() + if use_local_dof: + local_nodal_vector = nodal_vector + elif np.isscalar(mat_change_of_basis) and mat_change_of_basis == 1: + local_nodal_vector = nodal_vector + else: + n_nodal_dof = self.space.nvar * self.mesh.n_nodes + n_global_dof = mat_change_of_basis.shape[1] + if len(nodal_vector) == n_global_dof: + local_nodal_vector = mat_change_of_basis @ nodal_vector + elif len(nodal_vector) == n_nodal_dof and n_nodal_dof <= n_global_dof: + padded_vector = np.zeros(n_global_dof) + padded_vector[:n_nodal_dof] = nodal_vector + local_nodal_vector = mat_change_of_basis @ padded_vector + else: + raise ValueError( + "nodal_vector length is inconsistent with this assembly: " + f"got {len(nodal_vector)}, expected {n_nodal_dof} or " + f"{n_global_dof}. Pass use_local_dof=True for local " + "element dofs." + ) + associatedVariables = ( self._get_associated_variables() ) # for element requiring many variable such as beam with disp and rot dof @@ -1218,16 +1243,16 @@ def operator_apply(self, wf, nodal_vector): var.extend(associatedVariables[var[0]][0]) coef.extend(associatedVariables[var[0]][1]) - gp_values = ( - self._get_elementary_operator(wf.op[ii])[0] - @ nodal_vector[var[0] * n_nodes : (var[0] + 1) * n_nodes] - ) + elementary_operator = self._get_elementary_operator(wf.op[ii]) + n_operator_cols = elementary_operator[0].shape[1] + gp_values = elementary_operator[0] @ local_nodal_vector[ + var[0] * n_operator_cols : (var[0] + 1) * n_operator_cols + ] if len(var) > 1: for jj, v in enumerate(var[1:]): - gp_values += ( - self._get_elementary_operator(wf.op[ii])[jj + 1] - @ nodal_vector[v * n_nodes : (v + 1) * n_nodes] - ) + gp_values += elementary_operator[jj + 1] @ local_nodal_vector[ + v * n_operator_cols : (v + 1) * n_operator_cols + ] # or: # gp_values = self.get_gp_results(type(wf)([wf.op[ii]], [1], [1]), nodal_vector) @@ -1677,13 +1702,7 @@ def create(weakform, mesh="", elm_type="", name="", **kargs): ] list_assembly.append(Assembly(wf, mesh, elm_type, "", **kargs)) - if list_weakform[0] in l_wf: - assembly_output = list_assembly[ - -1 - ] # by default, the assembly used for output is the one associated to the 1st weakform - # list_assembly = [Assembly(wf, mesh, elm_type, "", **kargs) for wf in weakform.list_weakform] - kargs["assembly_output"] = kargs.get("assembly_output", assembly_output) return AssemblySum(list_assembly, name, **kargs) else: diff --git a/fedoo/core/assembly_sum.py b/fedoo/core/assembly_sum.py index 9f941a6f..eef4db5f 100644 --- a/fedoo/core/assembly_sum.py +++ b/fedoo/core/assembly_sum.py @@ -21,8 +21,6 @@ class AssemblySum(AssemblyBase): list of Assembly objects to sum name: str name of the Assembly - assembly_output: Assembly (optional keyword arg) - Assembly object used to extract output values (using Problem.get_results or Problem.save_results) """ def __init__(self, list_assembly, name="", **kargs): @@ -52,10 +50,6 @@ def __init__(self, list_assembly, name="", **kargs): self.current = _AssemblySumCurrent(list_assembly) - # for post-treatment only - self.__assembly_output = kargs.get("assembly_output", None) - if self.__assembly_output is not None: - self.sv = self.__assembly_output.sv # alias def __add__(self, another_assembly): if isinstance(another_assembly, AssemblySum): @@ -184,20 +178,12 @@ def _load_sv(self): def list_assembly(self): return self._list_assembly - @property - def assembly_output(self): - return self.__assembly_output - class _AssemblySumCurrent( AssemblySum ): # same as AssemblySum using the current assemblies and without output function def __init__(self, list_assembly, **kargs): self._list_assembly = list_assembly - self.__assembly_output = None - # self.__assembly_output = kargs.get('assembly_output', None) - # if self.__assembly_output is not None: self.sv = self.__assembly_output.sv #alias - self._reload = kargs.pop("reload", "all") AssemblyBase.__init__(self, name="") diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 66120465..0f21a786 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -7,6 +7,7 @@ from zipfile import ZipFile, Path as ZipPath from fedoo.core.mesh import Mesh, MultiMesh from fedoo.core.multimeshdata import MultiMeshData, copy_data_value +from fedoo.lib_elements.element_list import get_default_n_gp from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList try: @@ -953,19 +954,21 @@ def _multimesh_block_for_submesh( if block is not None or not fill_missing: return block - template = next( - (value for _, value in multimesh_data.items() if value is not None), + template_entry = next( + (item for item in multimesh_data.items() if item[1] is not None), None, ) - if template is None: + if template_entry is None: return None + template_id, template = template_entry template = np.asarray(template) if template.ndim > 1: shape = template.shape[:-1] + (n_items,) else: shape = (n_items,) - return np.full(shape, np.nan) + return np.zeros(shape, dtype=template.dtype) + def _submesh_dataset(self, submesh_id: int) -> "DataSet": """Build a lightweight DataSet view attached to one submesh.""" @@ -992,8 +995,7 @@ def _submesh_dataset(self, submesh_id: int) -> "DataSet": block := self._multimesh_block_for_submesh( value, submesh_id, - submesh.n_elements, - fill_missing=False, + submesh.n_elements * get_default_n_gp(submesh.elm_type, submesh), ) ) is not None diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py index c4802920..f502cce8 100644 --- a/fedoo/core/multimeshdata.py +++ b/fedoo/core/multimeshdata.py @@ -5,6 +5,7 @@ import numpy as np from fedoo.core.mesh import MultiMesh +from fedoo.lib_elements.element_list import get_default_n_gp class MultiMeshData: @@ -195,7 +196,7 @@ def global_element_value(self, element_id: int): return None return np.asarray(block)[..., local_id] - def global_element_values(self, element_ids, fill_value=np.nan): + def global_element_values(self, element_ids, fill_value=0.0): """Return values attached to several global element ids. Missing submesh blocks are filled with ``fill_value``. The returned @@ -226,14 +227,14 @@ def global_element_values(self, element_ids, fill_value=np.nan): return np.asarray(values) return np.stack(values, axis=-1) - def to_global(self, indices=None, fill_value=np.nan): + def to_global(self, indices=None, fill_value=0.0): """Concatenate per-submesh data in the selected submesh order.""" if indices is None: indices = list(range(len(self.mesh.submeshes))) else: indices = self._resolve_indices(indices) - template = next( - (self._data[i] for i in indices if self._data[i] is not None), + template_entry = next( + ((i, self._data[i]) for i in indices if self._data[i] is not None), None, ) arrays = [] @@ -241,8 +242,18 @@ def to_global(self, indices=None, fill_value=np.nan): data = self._data[i] n_elements = self.mesh[i].n_elements if data is None: - if template is not None and np.asarray(template).ndim > 1: - shape = np.asarray(template).shape[:-1] + (n_elements,) + if template_entry is not None and np.asarray(template_entry[1]).ndim > 1: + template_id, template = template_entry + template = np.asarray(template) + n_template_elements = self.mesh[template_id].n_elements + factor = template.shape[-1] // n_template_elements + if factor > 1: + n_missing = n_elements * get_default_n_gp( + self.mesh[i].elm_type, self.mesh[i] + ) + else: + n_missing = n_elements + shape = template.shape[:-1] + (n_missing,) else: shape = (n_elements,) arrays.append(np.full(shape, fill_value)) diff --git a/fedoo/core/output.py b/fedoo/core/output.py index 782224e5..ef23e8b3 100644 --- a/fedoo/core/output.py +++ b/fedoo/core/output.py @@ -114,12 +114,12 @@ def _output_mesh_for_assembly(assemb, element_set=None): """Return the mesh associated with an output request. For a regular assembly, this is its mesh, optionally restricted to - ``element_set``. For an ``AssemblySum`` without ``assembly_output``, this - builds a ``MultiMesh`` from the unique elementary assembly meshes, in + ``element_set``. For an ``AssemblySum``, this builds a ``MultiMesh`` from + the unique elementary assembly meshes, in ``list_assembly`` order. If all elementary assemblies share one mesh, the mesh itself is returned. """ - if hasattr(assemb, "list_assembly") and assemb.assembly_output is None: + if hasattr(assemb, "list_assembly"): entries = _unique_assembly_mesh_entries(assemb, element_set) meshes = [entry["mesh"] for entry in entries] if not meshes: @@ -132,14 +132,29 @@ def _output_mesh_for_assembly(assemb, element_set=None): register_name=False, ) - if hasattr(assemb, "list_assembly"): - assemb = assemb.assembly_output - if element_set is None: return assemb.mesh return assemb.mesh.extract_elements(element_set) +def _find_constitutivelaw_with_method(weakform, method_name): + """Return a constitutive law from a possibly nested weakform.""" + law = getattr(weakform, "constitutivelaw", None) + if law is not None and hasattr(law, method_name): + return law + + for child in getattr(weakform, "list_weakform", []): + law = _find_constitutivelaw_with_method(child, method_name) + if law is not None: + return law + + wrapped = getattr(weakform, "weakform", None) + if wrapped is not None: + return _find_constitutivelaw_with_method(wrapped, method_name) + + return None + + def _dataset_field(data_set, field): """Return a field and its data type from a DataSet, or (None, None).""" if field in data_set.node_data: @@ -286,7 +301,7 @@ def _get_results( if isinstance(assemb, str): assemb = AssemblyBase.get_all()[assemb] - if hasattr(assemb, "list_assembly") and assemb.assembly_output is None: + if hasattr(assemb, "list_assembly"): return _get_assemblysum_results( pb, assemb, @@ -312,9 +327,6 @@ def _get_results( data_sav = {} # dict to keep data in memory that may be used more that one time - if hasattr(assemb, "list_assembly"): # AssemblySum object - assemb = assemb.assembly_output - sv = assemb.sv # state variables associated to the assembly if include_mesh: @@ -356,19 +368,16 @@ def _get_results( data = sv[res] else: # attent to compute - try: - if res == "Strain": - data = assemb.weakform.constitutivelaw.get_strain( - assemb, position=position - ) - elif res == "Stress": - data = assemb.weakform.constitutivelaw.get_stress( - assemb, position=position - ) - else: - assert 0 - except: + method_name = "get_strain" if res == "Strain" else "get_stress" + law = _find_constitutivelaw_with_method(assemb.weakform, method_name) + if law is None: raise NameError('Field "{}" not available'.format(res)) + try: + data = getattr(law, method_name)(assemb, position=position) + except Exception as exc: + raise NameError( + 'Field "{}" not available'.format(res) + ) from exc # keep data in memory in case it may be used later for vm, pc or pdir stress computation data_sav[res] = data diff --git a/tests/test_distributed_load_neumann_bc.py b/tests/test_distributed_load_neumann_bc.py index 887c0d0d..799e3f0d 100644 --- a/tests/test_distributed_load_neumann_bc.py +++ b/tests/test_distributed_load_neumann_bc.py @@ -101,11 +101,74 @@ def counted_assemble_global_mat(compute="all"): pb.apply_boundary_conditions(t_fact=1.0, t_fact_old=0.5) full_load = pb.get_B().copy() - assert call_count == 2 + assert call_count == 1 np.testing.assert_allclose(quarter_load_again, quarter_load) np.testing.assert_allclose(quarter_load, 0.25 * full_load) +def test_assembly_neumann_load_caches_initial_and_final_vectors(): + fd.Assembly.delete_memory() + fd.ModelingSpace("3D") + + mesh = fd.mesh.box_mesh(2, 2, 2, name="initial_distributed_load_box") + material = fd.constitutivelaw.ElasticIsotrop(1000.0, 0.3, name="InitialLoadMat") + wf = fd.weakform.StressEquilibrium(material) + solid = fd.Assembly.create(wf, mesh, name="initial_load_solid") + pb = fd.problem.NonLinear(solid, nlgeom=False) + load = fd.constraint.DistributedForce( + mesh, + [0.0, 0.0, -10.0], + initial_force=[0.0, 0.0, -2.0], + nlgeom=False, + ) + bc = pb.bc.add(load) + assemble_global_mat = load.assemble_global_mat + call_count = 0 + + def counted_assemble_global_mat(compute="all"): + nonlocal call_count + call_count += 1 + return assemble_global_mat(compute) + + load.assemble_global_mat = counted_assemble_global_mat + + pb.apply_boundary_conditions(t_fact=0.25, t_fact_old=0.0) + quarter_load = pb.get_B().copy() + start_load = bc.start_value.copy() + final_load = bc.value.copy() + + pb.apply_boundary_conditions(t_fact=0.50, t_fact_old=0.25) + half_load = pb.get_B().copy() + + assert call_count == 2 + np.testing.assert_allclose( + quarter_load, start_load + 0.25 * (final_load - start_load) + ) + np.testing.assert_allclose( + half_load, start_load + 0.50 * (final_load - start_load) + ) + + +def test_follower_assembly_neumann_load_without_initial_value_skips_start_vector(): + pb, bc = _build_problem( + load_nlgeom=True, problem_nlgeom=True, load_kind="distributed" + ) + load = bc.assembly + set_start = load.set_start + load_factors = [] + + def counted_set_start(problem, t_fact=None): + load_factors.append(t_fact) + return set_start(problem, t_fact) + + load.set_start = counted_set_start + + pb.apply_boundary_conditions(t_fact=0.25, t_fact_old=0.0) + + assert load_factors == [1] + assert bc.start_value == 0 + + def test_raw_assembly_without_neumann_conversion_is_rejected(): fd.Assembly.delete_memory() fd.ModelingSpace("3D") From 50aa5ff0aff226c834840970031f532e95dd3eaa Mon Sep 17 00:00:00 2001 From: pruliere Date: Thu, 9 Jul 2026 14:30:56 +0200 Subject: [PATCH 08/14] ruff format --- fedoo/constitutivelaw/shell.py | 15 +++------------ fedoo/constraint/distributed_load.py | 1 + fedoo/core/assembly.py | 18 ++++++++++++------ fedoo/core/assembly_sum.py | 1 - fedoo/core/dataset.py | 1 - fedoo/core/multimeshdata.py | 9 +++++---- fedoo/core/output.py | 8 ++++---- tests/test_distributed_load_neumann_bc.py | 4 +--- 8 files changed, 26 insertions(+), 31 deletions(-) diff --git a/fedoo/constitutivelaw/shell.py b/fedoo/constitutivelaw/shell.py index 694967f5..4c72f3f2 100644 --- a/fedoo/constitutivelaw/shell.py +++ b/fedoo/constitutivelaw/shell.py @@ -77,12 +77,7 @@ def update(self, assembly, pb): # all terms are computed. Perhaps could be optimized by computed only the termes related to the associated weak form (eg shear for selective integration) assembly.sv["ShellStress"] = [ - sum( - [ - ShellStrain[j] * assembly.convert_data(H[i][j]) - for j in range(8) - ] - ) + sum([ShellStrain[j] * assembly.convert_data(H[i][j]) for j in range(8)]) for i in range(8) ] # H[i][j] are converted to gauss point excepted if scalar assembly.sv["ShellStrain"] = ShellStrain @@ -114,12 +109,8 @@ def get_strain(self, assembly, **kargs): if np.isscalar(ShellStrain) and ShellStrain == 0: zeros = np.zeros(assembly.n_gauss_points) return StrainTensorList([zeros.copy() for _ in range(6)]) - Strain[0] = ( - ShellStrain[0] + z * ShellStrain[4] - ) # epsXX -> membrane and bending - Strain[1] = ( - ShellStrain[1] - z * ShellStrain[3] - ) # epsYY -> membrane and bending + Strain[0] = ShellStrain[0] + z * ShellStrain[4] # epsXX -> membrane and bending + Strain[1] = ShellStrain[1] - z * ShellStrain[3] # epsYY -> membrane and bending Strain[3] = ShellStrain[2] # 2epsXY Strain[4:6] = ShellStrain[6:8] # 2epsXZ and 2epsYZ -> shear diff --git a/fedoo/constraint/distributed_load.py b/fedoo/constraint/distributed_load.py index d887ee94..9d20f030 100644 --- a/fedoo/constraint/distributed_load.py +++ b/fedoo/constraint/distributed_load.py @@ -120,6 +120,7 @@ def generate(self, problem: ProblemBase, t_fact=1, t_fact_old=None): self._format_current_value(problem, self.get_value(t_fact, t_fact_old)) return [self] + class Pressure(Assembly): """Pressure load. diff --git a/fedoo/core/assembly.py b/fedoo/core/assembly.py index bc05559a..203794bf 100644 --- a/fedoo/core/assembly.py +++ b/fedoo/core/assembly.py @@ -1245,14 +1245,20 @@ def operator_apply(self, wf, nodal_vector, use_local_dof=False): elementary_operator = self._get_elementary_operator(wf.op[ii]) n_operator_cols = elementary_operator[0].shape[1] - gp_values = elementary_operator[0] @ local_nodal_vector[ - var[0] * n_operator_cols : (var[0] + 1) * n_operator_cols - ] + gp_values = ( + elementary_operator[0] + @ local_nodal_vector[ + var[0] * n_operator_cols : (var[0] + 1) * n_operator_cols + ] + ) if len(var) > 1: for jj, v in enumerate(var[1:]): - gp_values += elementary_operator[jj + 1] @ local_nodal_vector[ - v * n_operator_cols : (v + 1) * n_operator_cols - ] + gp_values += ( + elementary_operator[jj + 1] + @ local_nodal_vector[ + v * n_operator_cols : (v + 1) * n_operator_cols + ] + ) # or: # gp_values = self.get_gp_results(type(wf)([wf.op[ii]], [1], [1]), nodal_vector) diff --git a/fedoo/core/assembly_sum.py b/fedoo/core/assembly_sum.py index eef4db5f..b1e40ded 100644 --- a/fedoo/core/assembly_sum.py +++ b/fedoo/core/assembly_sum.py @@ -50,7 +50,6 @@ def __init__(self, list_assembly, name="", **kargs): self.current = _AssemblySumCurrent(list_assembly) - def __add__(self, another_assembly): if isinstance(another_assembly, AssemblySum): return AssemblySum(self.list_assembly + another_assembly.list_assembly) diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 0f21a786..fb9dd8de 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -969,7 +969,6 @@ def _multimesh_block_for_submesh( shape = (n_items,) return np.zeros(shape, dtype=template.dtype) - def _submesh_dataset(self, submesh_id: int) -> "DataSet": """Build a lightweight DataSet view attached to one submesh.""" submesh = self.mesh[submesh_id] diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py index f502cce8..bd61799d 100644 --- a/fedoo/core/multimeshdata.py +++ b/fedoo/core/multimeshdata.py @@ -216,9 +216,7 @@ def global_element_values(self, element_ids, fill_value=0.0): if template is None or np.asarray(template).ndim == 1: values.append(fill_value) else: - values.append( - np.full(np.asarray(template).shape[:-1], fill_value) - ) + values.append(np.full(np.asarray(template).shape[:-1], fill_value)) else: values.append(np.asarray(block)[..., local_id]) if not values: @@ -242,7 +240,10 @@ def to_global(self, indices=None, fill_value=0.0): data = self._data[i] n_elements = self.mesh[i].n_elements if data is None: - if template_entry is not None and np.asarray(template_entry[1]).ndim > 1: + if ( + template_entry is not None + and np.asarray(template_entry[1]).ndim > 1 + ): template_id, template = template_entry template = np.asarray(template) n_template_elements = self.mesh[template_id].n_elements diff --git a/fedoo/core/output.py b/fedoo/core/output.py index ef23e8b3..3020f0a8 100644 --- a/fedoo/core/output.py +++ b/fedoo/core/output.py @@ -369,15 +369,15 @@ def _get_results( else: # attent to compute method_name = "get_strain" if res == "Strain" else "get_stress" - law = _find_constitutivelaw_with_method(assemb.weakform, method_name) + law = _find_constitutivelaw_with_method( + assemb.weakform, method_name + ) if law is None: raise NameError('Field "{}" not available'.format(res)) try: data = getattr(law, method_name)(assemb, position=position) except Exception as exc: - raise NameError( - 'Field "{}" not available'.format(res) - ) from exc + raise NameError('Field "{}" not available'.format(res)) from exc # keep data in memory in case it may be used later for vm, pc or pdir stress computation data_sav[res] = data diff --git a/tests/test_distributed_load_neumann_bc.py b/tests/test_distributed_load_neumann_bc.py index 799e3f0d..706d0e3a 100644 --- a/tests/test_distributed_load_neumann_bc.py +++ b/tests/test_distributed_load_neumann_bc.py @@ -144,9 +144,7 @@ def counted_assemble_global_mat(compute="all"): np.testing.assert_allclose( quarter_load, start_load + 0.25 * (final_load - start_load) ) - np.testing.assert_allclose( - half_load, start_load + 0.50 * (final_load - start_load) - ) + np.testing.assert_allclose(half_load, start_load + 0.50 * (final_load - start_load)) def test_follower_assembly_neumann_load_without_initial_value_skips_start_vector(): From 538b11f9ccd6aa06aa4e7e9a73a49616c8cd6df1 Mon Sep 17 00:00:00 2001 From: pruliere Date: Thu, 9 Jul 2026 17:36:46 +0200 Subject: [PATCH 09/14] remove confusion between local and global dof in assembly saved data meshChange -> mesh_change --- docs/simple_examples.rst | 2 +- examples/FE2/FE2.py | 2 +- .../2DplasticBending_from_scratch.py | 4 +- examples/plasticity/plastic_bending_3D.py | 4 +- examples/plasticity/shear_test.py | 4 +- .../thermal/thermo_meca_weak_coupling_3D.py | 2 +- fedoo/core/assembly.py | 77 +++++++++++++------ fedoo/core/mesh.py | 17 +++- fedoo/util/simple_plot.py | 2 +- tests/test_2DDynamicPlasticBending.py | 4 +- tests/test_2DDynamicPlasticBending_v2.py | 4 +- tests/test_BeamElement.py | 6 +- tests/test_distributed_load_neumann_bc.py | 31 ++++++++ tests/test_factor_reuse.py | 2 +- tests/test_platewithhol.py | 2 +- 15 files changed, 118 insertions(+), 45 deletions(-) diff --git a/docs/simple_examples.rst b/docs/simple_examples.rst index 3e303ab9..3b1eb89f 100644 --- a/docs/simple_examples.rst +++ b/docs/simple_examples.rst @@ -24,7 +24,7 @@ ___________________________ fd.weakform.StressEquilibrium("ElasticLaw", name = "WeakForm") #Create a global assembly - fd.Assembly.create("WeakForm", "Domain", name="Assembly", MeshChange = True) + fd.Assembly.create("WeakForm", "Domain", name="Assembly", mesh_change = True) #Define a new static problem pb = fd.problem.Linear("Assembly") diff --git a/examples/FE2/FE2.py b/examples/FE2/FE2.py index eab64e26..0d82fbe7 100644 --- a/examples/FE2/FE2.py +++ b/examples/FE2/FE2.py @@ -46,7 +46,7 @@ # Create a global assembly macro_assembly = fd.Assembly.create( - "WeakForm", "macro", name="Assembly", MeshChange=True + "WeakForm", "macro", name="Assembly", mesh_change=True ) # Define a new static problem diff --git a/examples/plasticity/2DplasticBending_from_scratch.py b/examples/plasticity/2DplasticBending_from_scratch.py index 70f72a11..89e03d4e 100644 --- a/examples/plasticity/2DplasticBending_from_scratch.py +++ b/examples/plasticity/2DplasticBending_from_scratch.py @@ -52,8 +52,8 @@ top_center = np.hstack((nodes_top1, nodes_top2)) fd.Assembly.create( - wf, "Domain", "quad4", name="Assembling", MeshChange=False -) # uses MeshChange=True when the mesh change during the time + wf, "Domain", "quad4", name="Assembling", mesh_change=False +) # uses mesh_change=True when the mesh change during the time pb = fd.problem.NonLinear("Assembling") # incremental non linear problems diff --git a/examples/plasticity/plastic_bending_3D.py b/examples/plasticity/plastic_bending_3D.py index a652017f..4c01db63 100644 --- a/examples/plasticity/plastic_bending_3D.py +++ b/examples/plasticity/plastic_bending_3D.py @@ -80,8 +80,8 @@ nodes_top2 = mesh.find_nodes("XY", (3 * L / 4, h)) nodes_topCenter = np.hstack((nodes_top1, nodes_top2)) -# Assembly.create("constitutivelaw", meshname, 'hex8', name="Assembling", MeshChange = False, n_elm_gp = 27) #uses MeshChange=True when the mesh change during the time -# assemb = fd.Assembly.create("constitutivelaw", meshname, 'hex8', name="Assembling", MeshChange = False, n_elm_gp = 8) #uses MeshChange=True when the mesh change during the time +# Assembly.create("constitutivelaw", meshname, 'hex8', name="Assembling", mesh_change = False, n_elm_gp = 27) #uses mesh_change=True when the mesh change during the time +# assemb = fd.Assembly.create("constitutivelaw", meshname, 'hex8', name="Assembling", mesh_change = False, n_elm_gp = 8) #uses mesh_change=True when the mesh change during the time assemb = fd.Assembly.create(wf, meshname, name="Assembling") diff --git a/examples/plasticity/shear_test.py b/examples/plasticity/shear_test.py index 8f091668..3957519b 100644 --- a/examples/plasticity/shear_test.py +++ b/examples/plasticity/shear_test.py @@ -66,10 +66,10 @@ wf.fbar = True wf.corate = "log" -# fd.Assembly.create("ConstitutiveLaw", meshname, 'hex8', name="Assembling", MeshChange = False, n_elm_gp = 27) #uses MeshChange=True when the mesh change during the time +# fd.Assembly.create("ConstitutiveLaw", meshname, 'hex8', name="Assembling", mesh_change = False, n_elm_gp = 27) #uses mesh_change=True when the mesh change during the time fd.Assembly.create( wf, meshname, name="Assembling" -) # uses MeshChange=True when the mesh change during the time +) # uses mesh_change=True when the mesh change during the time pb = fd.problem.NonLinear("Assembling") pb.set_nr_criterion("Force", err0=None, tol=5e-3, max_subiter=10) diff --git a/examples/thermal/thermo_meca_weak_coupling_3D.py b/examples/thermal/thermo_meca_weak_coupling_3D.py index f65a5709..bcd1aab1 100644 --- a/examples/thermal/thermo_meca_weak_coupling_3D.py +++ b/examples/thermal/thermo_meca_weak_coupling_3D.py @@ -82,7 +82,7 @@ fd.Assembly.create( wf_mech, meshname, name="Assembling_M" -) # uses MeshChange=True when the mesh change during the time +) # uses mesh_change=True when the mesh change during the time pb_m = fd.problem.NonLinear("Assembling_M") # pb_m.set_solver('cg', precond = True) diff --git a/fedoo/core/assembly.py b/fedoo/core/assembly.py index 203794bf..d4d17efb 100644 --- a/fedoo/core/assembly.py +++ b/fedoo/core/assembly.py @@ -92,7 +92,7 @@ def __init__(self, weakform, mesh="", elm_type="", name="", **kargs): # used for update lagrangian method. self.current = self - self.meshChange = kargs.pop("MeshChange", False) + self.mesh_change = kargs.pop("mesh_change", False) if elm_type == "" or elm_type is None: elm_type = weakform.assembly_options.get( "elm_type", mesh.elm_type, mesh.elm_type @@ -155,10 +155,11 @@ def assemble_global_mat(self, compute="all"): n_elm_gp = self.n_elm_gp if ( - self.meshChange == True + self.mesh_change == True ): # only node position change is considered here. For change of sparsity, use also, self.mesh.init_interpolation - if self.mesh in Assembly._saved_change_of_basis_mat: - del Assembly._saved_change_of_basis_mat[self.mesh] + for k in tuple(Assembly._saved_change_of_basis_mat): + if k[0] == self.mesh: + Assembly._saved_change_of_basis_mat.pop(k) self.compute_elementary_operators() nvar = self.space.nvar @@ -557,8 +558,13 @@ def get_change_of_basis_mat(self): mesh = self.mesh - if mesh in Assembly._saved_change_of_basis_mat: - return Assembly._saved_change_of_basis_mat[mesh] + change_of_basis_key = ( + mesh, + self._use_local_csys, + mesh._local_frame_cache_key(self._element_local_frame), + ) + if change_of_basis_key in Assembly._saved_change_of_basis_mat: + return Assembly._saved_change_of_basis_mat[change_of_basis_key] mat_change_of_basis = 1 compute_mat_change_of_basis = False @@ -675,7 +681,7 @@ def get_change_of_basis_mat(self): mat_change_of_basis = mat_change_of_basis.tocsr() - Assembly._saved_change_of_basis_mat[mesh] = mat_change_of_basis + Assembly._saved_change_of_basis_mat[change_of_basis_key] = mat_change_of_basis return mat_change_of_basis def initialize(self, pb): @@ -773,7 +779,7 @@ def delete_memory(): but it is not done by default. In this case, deleting the memory should resolve the problem. - Note: it the MeshChange argument is set to True when creating the Assembly object, the + Note: it the mesh_change argument is set to True when creating the Assembly object, the memory will be recomputed by default, which may cause a decrease in assembling performances """ Assembly._saved_elementary_operators = {} @@ -808,19 +814,27 @@ def compute_elementary_operators( # we compute the operators directly from the element library elmRef = get_element(elm_type)(n_elm_gp) OP = elmRef.computeOperator(nodes, elements) - mesh._saved_gaussian_quadrature_mat[n_elm_gp] = sparse.identity( + mesh._saved_gaussian_quadrature_mat[ + mesh._gaussian_quadrature_cache_key(n_elm_gp) + ] = sparse.identity( OP[0, 0][0].shape[0], "d", format="csr" ) # No gaussian quadrature in this case : nodal identity matrix mesh._saved_gausspoint2node_mat[n_elm_gp] = ( 1 # no need to translate between pg and nodes because no pg ) mesh._saved_node2gausspoint_mat[n_elm_gp] = 1 - Assembly._saved_change_of_basis_mat[mesh] = ( - 1 # No change of basis: mat_change_of_basis = 1 #this line could be deleted because the coordinate should in principle defined as 'global' - ) - Assembly._saved_elementary_operators[(mesh, elm_type, n_elm_gp)] = ( - OP # elmRef.computeOperator(nodes,elements) - ) + Assembly._saved_change_of_basis_mat[ + (mesh, self._use_local_csys, mesh._local_frame_cache_key(None)) + ] = 1 # No change of basis: mat_change_of_basis = 1 + Assembly._saved_elementary_operators[ + ( + mesh, + elm_type, + n_elm_gp, + self._use_local_csys, + mesh._local_frame_cache_key(self._element_local_frame), + ) + ] = OP # elmRef.computeOperator(nodes,elements) return # ------------------------------------------------------------------- @@ -935,7 +949,15 @@ def compute_elementary_operators( i + n_diff_interpolations ] # as index and indptr should be the same, perhaps it will be more memory efficient to only store the data field - Assembly._saved_elementary_operators[(mesh, elm_type.name, n_elm_gp)] = data + Assembly._saved_elementary_operators[ + ( + mesh, + elm_type.name, + n_elm_gp, + self._use_local_csys, + mesh._local_frame_cache_key(self._element_local_frame), + ) + ] = data def _get_elementary_operator(self, deriv, n_elm_gp=None): # Gives a list of sparse matrix that convert node values for one variable to the pg values of a simple derivative op (for instance d/dz) @@ -950,10 +972,17 @@ def _get_elementary_operator(self, deriv, n_elm_gp=None): if hasattr(get_element(elm_type), "get_elm_type"): elm_type = get_element(elm_type).get_elm_type(deriv.u_name).name - if not ((mesh, elm_type, n_elm_gp) in Assembly._saved_elementary_operators): + key = ( + mesh, + elm_type, + n_elm_gp, + self._use_local_csys, + mesh._local_frame_cache_key(self._element_local_frame), + ) + if key not in Assembly._saved_elementary_operators: self.compute_elementary_operators(n_elm_gp) - data = Assembly._saved_elementary_operators[(mesh, elm_type, n_elm_gp)] + data = Assembly._saved_elementary_operators[key] if deriv.ordre == 0: # in this case deriv.x should be 0 execpt if several interpolations @@ -982,9 +1011,12 @@ def _get_elementary_operator(self, deriv, n_elm_gp=None): def _get_gaussian_quadrature_mat( self, ): # calcul la discrétision relative à un seul opérateur dérivé - if not (self.n_elm_gp in self.mesh._saved_gaussian_quadrature_mat): + cache_key = self.mesh._gaussian_quadrature_cache_key( + self.n_elm_gp, self._element_local_frame + ) + if cache_key not in self.mesh._saved_gaussian_quadrature_mat: self.compute_elementary_operators() - return self.mesh._saved_gaussian_quadrature_mat[self.n_elm_gp] + return self.mesh._saved_gaussian_quadrature_mat[cache_key] def _get_associated_variables( self, @@ -1330,8 +1362,9 @@ def set_disp(self, disp): self.current.mesh._saved_gaussian_quadrature_mat = {} self.current.mesh._saved_gausspoint2node_l2 = {} - if self.current.mesh in self._saved_change_of_basis_mat: - del self._saved_change_of_basis_mat[self.current.mesh] + for k in tuple(self._saved_change_of_basis_mat): + if k[0] == self.current.mesh: + self._saved_change_of_basis_mat.pop(k) for k in tuple(self._saved_elementary_operators): if k[0] == self.current.mesh: diff --git a/fedoo/core/mesh.py b/fedoo/core/mesh.py index 165d6b03..216e47b5 100644 --- a/fedoo/core/mesh.py +++ b/fedoo/core/mesh.py @@ -1230,6 +1230,13 @@ def init_interpolation(self, n_elm_gp: int | None = None) -> None: self._sparse_structure[n_elm_gp] = (row.reshape(-1), col.reshape(-1)) self._elements_geom = elm_geom # dont depend on n_elm_gp + @staticmethod + def _local_frame_cache_key(local_frame): + return None if local_frame is None else id(local_frame) + + def _gaussian_quadrature_cache_key(self, n_elm_gp, local_frame=None): + return (n_elm_gp, self._local_frame_cache_key(local_frame)) + def _compute_gaussian_quadrature_mat( self, n_elm_gp: int | None = None, @@ -1237,7 +1244,8 @@ def _compute_gaussian_quadrature_mat( ) -> None: if n_elm_gp is None: n_elm_gp = get_default_n_gp(self.elm_type) - if n_elm_gp not in self._saved_gaussian_quadrature_mat: + cache_key = self._gaussian_quadrature_cache_key(n_elm_gp, local_frame) + if cache_key not in self._saved_gaussian_quadrature_mat: self.init_interpolation(n_elm_gp) elm_interpol = self._elm_interpolation[n_elm_gp] @@ -1252,7 +1260,7 @@ def _compute_gaussian_quadrature_mat( # Compute the diag matrix used for the gaussian quadrature # ------------------------------------------------------------------- gaussianQuadrature = (elm_interpol.detJ * elm_interpol.w_pg).T.reshape(-1) - self._saved_gaussian_quadrature_mat[n_elm_gp] = sparse.diags( + self._saved_gaussian_quadrature_mat[cache_key] = sparse.diags( gaussianQuadrature, 0, format="csr" ) # matrix to get the gaussian quadrature (integration over each element) @@ -1404,9 +1412,10 @@ def _get_node2gausspoint_mat(self, n_elm_gp=None): def _get_gaussian_quadrature_mat(self, n_elm_gp=None): if n_elm_gp is None: n_elm_gp = get_default_n_gp(self.elm_type) - if not (n_elm_gp in self._saved_gaussian_quadrature_mat): + cache_key = self._gaussian_quadrature_cache_key(n_elm_gp) + if cache_key not in self._saved_gaussian_quadrature_mat: self._compute_gaussian_quadrature_mat(n_elm_gp) - return self._saved_gaussian_quadrature_mat[n_elm_gp] + return self._saved_gaussian_quadrature_mat[cache_key] def determine_data_type(self, data: np.ndarray, n_elm_gp: int | None = None) -> str: if n_elm_gp is None: diff --git a/fedoo/util/simple_plot.py b/fedoo/util/simple_plot.py index cc29b6d0..e1eda92b 100644 --- a/fedoo/util/simple_plot.py +++ b/fedoo/util/simple_plot.py @@ -203,7 +203,7 @@ def field_plot_2d( U = ((disp.reshape(2, -1).T[elm.ravel()]).T).ravel() # reload the assembly with the new mesh - assemb_visu = Assembly(wf, "visu", type_el, name="visu", MeshChange=True) + assemb_visu = Assembly(wf, "visu", type_el, name="visu", mesh_change=True) assemb_visu.compute_elementary_operators() else: raise NameError("type_plot should be either 'real' or 'smooth'") diff --git a/tests/test_2DDynamicPlasticBending.py b/tests/test_2DDynamicPlasticBending.py index 734c87b6..fe6184e7 100644 --- a/tests/test_2DDynamicPlasticBending.py +++ b/tests/test_2DDynamicPlasticBending.py @@ -56,8 +56,8 @@ def test_2DDynamicPlasticBending(): nodes_top2 = np.where((crd[:, 0] == 3 * L / 4) * (crd[:, 1] == h))[0] fd.Assembly.create( - "weakform", "Domain", "quad4", name="Assembling", MeshChange=False - ) # uses MeshChange=True when the mesh change during the time + "weakform", "Domain", "quad4", name="Assembling", mesh_change=False + ) # uses mesh_change=True when the mesh change during the time pb = fd.problem.NonLinearNewmark("Assembling", 0.25, 0.5) diff --git a/tests/test_2DDynamicPlasticBending_v2.py b/tests/test_2DDynamicPlasticBending_v2.py index ed5f8367..2db4b5bf 100644 --- a/tests/test_2DDynamicPlasticBending_v2.py +++ b/tests/test_2DDynamicPlasticBending_v2.py @@ -58,8 +58,8 @@ def test_2DDynamicPlasticBending_v2(): nodes_top2 = np.where((crd[:, 0] == 3 * L / 4) * (crd[:, 1] == h))[0] fd.Assembly.create( - wf, "Domain", "quad4", name="Assembling", MeshChange=False - ) # uses MeshChange=True when the mesh change during the time + wf, "Domain", "quad4", name="Assembling", mesh_change=False + ) # uses mesh_change=True when the mesh change during the time pb = fd.problem.NonLinear("Assembling") diff --git a/tests/test_BeamElement.py b/tests/test_BeamElement.py index 2f52addd..2489530b 100644 --- a/tests/test_BeamElement.py +++ b/tests/test_BeamElement.py @@ -43,13 +43,13 @@ def test_beam_element(): "ElasticLaw", Section, Jx, Iyy, Izz, name="WFbeam" ) # by default k=0 i.e. no shear effect fd.Assembly.create( - "WFbeam", "beam", "bernoullibeam", name="beam", MeshChange=True + "WFbeam", "beam", "bernoullibeam", name="beam", mesh_change=True ) elif computeShear == 1: fd.weakform.BeamEquilibrium( "ElasticLaw", Section, Jx, Iyy, Izz, k=k, name="WFbeam" ) - fd.Assembly.create("WFbeam", "beam", "beam", name="beam", MeshChange=True) + fd.Assembly.create("WFbeam", "beam", "beam", name="beam", mesh_change=True) else: # computeShear = 2 fd.Mesh["beam"].add_internal_nodes( 1 @@ -58,7 +58,7 @@ def test_beam_element(): "ElasticLaw", Section, Jx, Iyy, Izz, k=k, name="WFbeam" ) fd.Assembly.create( - "WFbeam", "beam", "beamfcq", name="beam", MeshChange=True + "WFbeam", "beam", "beamfcq", name="beam", mesh_change=True ) pb = fd.problem.Linear("beam") diff --git a/tests/test_distributed_load_neumann_bc.py b/tests/test_distributed_load_neumann_bc.py index 706d0e3a..d81c723f 100644 --- a/tests/test_distributed_load_neumann_bc.py +++ b/tests/test_distributed_load_neumann_bc.py @@ -147,6 +147,37 @@ def counted_assemble_global_mat(compute="all"): np.testing.assert_allclose(half_load, start_load + 0.50 * (final_load - start_load)) +def test_distributed_load_after_shell_operator_cache_is_not_dropped(): + fd.Assembly.delete_memory() + fd.ModelingSpace("3D") + + nodes = np.array( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0], + ] + ) + elements = np.array([[0, 1, 2], [0, 2, 3]]) + mesh = fd.Mesh(nodes, elements, "tri3", ndim=3) + material = fd.constitutivelaw.ElasticIsotrop(1000.0, 0.3) + shell = fd.constitutivelaw.ShellHomogeneous(material, 0.1) + shell_wf = fd.weakform.PlateEquilibrium(shell, nlgeom=True) + shell_assembly = fd.Assembly.create(shell_wf, mesh, name="shell_cache_test") + pb = fd.problem.NonLinear(shell_assembly, nlgeom=True) + load = fd.constraint.DistributedForce(mesh, [0.0, 0.0, -1.0], nlgeom=False) + bc = pb.bc.add(load) + + pb.initialize() + shell_assembly.assemble_global_mat(compute="matrix") + pb.apply_boundary_conditions(t_fact=1.0, t_fact_old=0.0) + + z_dofs = slice(2 * mesh.n_nodes, 3 * mesh.n_nodes) + assert np.linalg.norm(bc.value[z_dofs]) > 0.0 + np.testing.assert_allclose(pb.get_B()[z_dofs], bc.value[z_dofs]) + + def test_follower_assembly_neumann_load_without_initial_value_skips_start_vector(): pb, bc = _build_problem( load_nlgeom=True, problem_nlgeom=True, load_kind="distributed" diff --git a/tests/test_factor_reuse.py b/tests/test_factor_reuse.py index beee2293..cefdae0c 100644 --- a/tests/test_factor_reuse.py +++ b/tests/test_factor_reuse.py @@ -23,7 +23,7 @@ def _build_plate_problem(): ) fd.constitutivelaw.ElasticIsotrop(2e5, 0.3, name="ElasticLaw") fd.weakform.StressEquilibrium("ElasticLaw", name="WeakForm") - fd.Assembly.create("WeakForm", "Domain", name="Assembly", MeshChange=True) + fd.Assembly.create("WeakForm", "Domain", name="Assembly", mesh_change=True) pb = fd.problem.Linear("Assembly") mesh = fd.Mesh["Domain"] diff --git a/tests/test_platewithhol.py b/tests/test_platewithhol.py index 4011c198..5a2c3130 100644 --- a/tests/test_platewithhol.py +++ b/tests/test_platewithhol.py @@ -19,7 +19,7 @@ def test_plate_with_hole(): fd.weakform.StressEquilibrium("ElasticLaw", name="WeakForm") # Create a global assembly - fd.Assembly.create("WeakForm", "Domain", name="Assembly", MeshChange=True) + fd.Assembly.create("WeakForm", "Domain", name="Assembly", mesh_change=True) # Define a new static problem pb = fd.problem.Linear("Assembly") From 497a892e739bbfe2210f250ed899640041608c28 Mon Sep 17 00:00:00 2001 From: pruliere Date: Thu, 9 Jul 2026 17:38:36 +0200 Subject: [PATCH 10/14] remove deprecated simple_plot.py --- fedoo/util/simple_plot.py | 276 -------------------------------------- 1 file changed, 276 deletions(-) delete mode 100644 fedoo/util/simple_plot.py diff --git a/fedoo/util/simple_plot.py b/fedoo/util/simple_plot.py deleted file mode 100644 index e1eda92b..00000000 --- a/fedoo/util/simple_plot.py +++ /dev/null @@ -1,276 +0,0 @@ -try: - import matplotlib.pyplot as plt - import matplotlib.tri as mtri - from matplotlib.patches import Polygon - from matplotlib.collections import PatchCollection - - USE_MPL = True -except ImportError: - USE_MPL = False - - -import numpy as np -from fedoo.core.mesh import Mesh -from fedoo.core.assembly import Assembly -from fedoo.core.base import ConstitutiveLaw - - -def mesh_plot_2d( - mesh, - disp=None, - data=None, - data_min=None, - data_max=None, - scale_factor=1, - plot_edge=True, - nb_level=6, - cm="hsv", -): - """ - Simple function for ploting 2D mesh with node data and node displacement - - Parameters - ---------- - mesh : TYPE - DESCRIPTION. - disp : TYPE - DESCRIPTION. - data : TYPE - DESCRIPTION. - scale_factor : TYPE, optional - DESCRIPTION. The default is 1. - plot_edge : TYPE, optional - DESCRIPTION. The default is False. - nb_level : TYPE, optional - DESCRIPTION. The default is 6. - - Returns - ------- - None. - - """ - if not (USE_MPL): - raise NameError( - "Matplotlib need to be installed to use the function mesh_plot_2D." - ) - - # Create triangulation. - if cm == "binary": - color = plt.cm.binary - else: - color = plt.cm.hsv - - if isinstance(mesh, str): - mesh = Mesh.get_all()[mesh] - - crd = mesh.nodes - elm = mesh.elements - type_el = mesh.elm_type - - if disp is not None: - # Get the displacement vector on nodes for export to vtk - U = np.reshape(disp, (2, -1)).T - N = mesh.n_nodes - U = np.c_[U, np.zeros(N)] - x = crd[:, 0] + U[:, 0] * scale_factor - y = crd[:, 1] + U[:, 1] * scale_factor - else: - x = crd[:, 0] - y = crd[:, 1] - - if type_el in ["tri3", "tri6"]: - triang = mtri.Triangulation(x, y, elm[:, 0:3]) - else: - triang = mtri.Triangulation(x, y, np.vstack((elm[:, 0:3], elm[:, [0, 2, 3]]))) - - fig = plt.figure() - ax = fig.gca() - - if data is not None: - if data_min is None and data_max is None: - plt.tricontourf(triang, data, nb_level, cmap=color) - else: - if data_min is None: - data_min = data.min() - if data_max is None: - data_max = data.max() - plt.tricontourf( - triang, - data, - vmax=data_max, - levels=np.linspace(data_min, data_max, nb_level), - extend="both", - cmap=color, - ) - - if plot_edge == True: - if type_el in ["tri3", "tri6"]: - plt.triplot(triang, lw=0.5, color="k") - else: # quad - crd_scaled = np.c_[x, y] # crd[:,0:2] + U[:,0:2]*50 - Nnde_elm = len(elm[0]) - patches = [] - - for i in range(len(elm)): - polygon = Polygon(crd_scaled[elm[i, 0:4], 0:2], closed=True) - patches.append(polygon) - - p = PatchCollection(patches, edgecolors="k", fc="None", lw=0.5, alpha=1) - ax.add_collection(p) - ax.set_xlim(crd_scaled[:, 0].min(), crd_scaled[:, 0].max()) - ax.set_ylim(crd_scaled[:, 1].min(), crd_scaled[:, 1].max()) - - ax.set_aspect("equal") - # ax.set_title('Stress') - # ax = fig.gca() - - # fig, ax = plt.subplots() - if data is not None: - plt.colorbar(orientation="horizontal") - # plt.clim(data_min, data_max) - plt.xlabel("x") - plt.ylabel("y") - # plt.ion() - - -# def fieldPlot2d(mesh, Matname, disp, dataname =None, component=0, data_min=None,data_max=None, scale_factor = 1, plot_edge = True, nb_level = 6, type_plot = "real", cm = 'hsv'): -def field_plot_2d( - assemb, - disp, - dataname=None, - component=0, - data_min=None, - data_max=None, - scale_factor=1, - plot_edge=True, - nb_level=6, - type_plot="real", - cm="hsv", -): - if not (USE_MPL): - raise NameError( - "Matplotlib need to be installed to use the function field_plot_2d." - ) - - if isinstance(assemb, str): - assemb = Assembly.get_all()[assemb] - mesh = assemb.mesh - wf = assemb.weakform - - # type_plot is "real" or "smooth" - if dataname is None: # no data, just plot mesh - mesh_plot_2d(mesh, disp, None, None, None, scale_factor, plot_edge, nb_level) - return - elif ( - dataname.lower() == "disp" - ): # if data is disp, no difference between "real" and "smooth" plot - try: - mesh_plot_2d( - mesh, - disp, - disp.reshape(2, -1)[component], - data_min, - data_max, - scale_factor, - plot_edge, - nb_level, - ) - except: - raise NameError( - "Dataname: " - + str(dataname) - + " and component: " - + str(component) - + " doesn't exist" - ) - - plt.gca().set_title(dataname + "_" + str(component)) - return - - crd = mesh.nodes - elm = mesh.elements - type_el = mesh.elm_type - - if type_plot.lower() == "smooth": - # mesh2 = Mesh(crd, elm, type_el, name ='visu') - # crd2=crd ; elm2=elm - U = disp.ravel() - assemb_visu = assemb - elif type_plot.lower() == "real": - crd2 = crd[elm.ravel()] - elm2 = np.arange(elm.shape[0] * elm.shape[1]).reshape(-1, elm.shape[1]) - mesh2 = Mesh(crd2, elm2, type_el, name="visu") - U = ((disp.reshape(2, -1).T[elm.ravel()]).T).ravel() - - # reload the assembly with the new mesh - assemb_visu = Assembly(wf, "visu", type_el, name="visu", mesh_change=True) - assemb_visu.compute_elementary_operators() - else: - raise NameError("type_plot should be either 'real' or 'smooth'") - - # compute tensorstrain and tensorstress - # TensorStrain = assemb_visu.get_strain(U, "Nodal", nlgeom = False) - # TensorStress = ConstitutiveLaw.get_all()[Matname].GetStressFromStrain(TensorStrain) - - TensorStrain = assemb_visu.get_strain(U, "GaussPoint", nlgeom=False) - TensorStress = wf.constitutivelaw.get_stress_from_strain(assemb_visu, TensorStrain) - TensorStrain = TensorStrain.convert(assemb_visu, convert_to="Node") - TensorStress = TensorStress.convert(assemb_visu, convert_to="Node") - - try: - if dataname.lower() == "stress": - if isinstance(component, str) and component.lower() == "vm": - data = TensorStress.vonMises() - else: - data = TensorStress[component] - elif dataname.lower() == "strain": - data = TensorStrain[component] - except: - raise NameError( - "Dataname: " - + str(dataname) - + " and component: " - + str(component) - + " doesn't exist" - ) - - if type_plot.lower() == "smooth": - mesh_plot_2d( - mesh, U, data, data_min, data_max, scale_factor, plot_edge, nb_level, cm - ) - else: # type_plot.lower() == "real": - mesh_plot_2d( - "visu", U, data, data_min, data_max, scale_factor, plot_edge, nb_level, cm - ) - - # # Create triangulation. - # color = plt.cm.hsv - - # U = np.reshape(U,(2,-1)).T - # N = mesh2.n_nodes - # U = np.c_[U,np.zeros(N)] - - # x = crd2[:,0] + U[:,0]*scale_factor - # y = crd2[:,1] + U[:,1]*scale_factor - # triang = mtri.Triangulation(x, y, elm2[:,0:3]) - - # #plot the new mesh - # # Set up the figure - # fig = plt.figure() - - # # color = plt.cm.get_cmap('hsv',256) - # color = plt.cm.hsv - # nb_level = 10 - # plt.tricontourf(triang, data, nb_level, cmap=color) - # if plot_edge == True: - # plt.triplot(triang, lw=0.2, color='k') - # ax = fig.gca() - # ax.set_aspect('equal') - # ax.set_title('Triangular grid') - # # ax = fig.gca() - - # # fig, ax = plt.subplots() - # plt.colorbar(orientation = 'horizontal') - # plt.xlabel('x') - # plt.ylabel('y') - plt.gca().set_title(dataname + "_" + str(component)) From 6700835a5d9c19f2134855ecd7b9bd4fe18a50c6 Mon Sep 17 00:00:00 2001 From: Yves Chemisky Date: Fri, 10 Jul 2026 00:20:05 +0200 Subject: [PATCH 11/14] Improve FDH5/XDMF docs and caching Add documentation, minor API clarifications and behavior fixes across IO and core modules. Key changes: - environment: add meshio dependency. - constraint/distributed_load: cache load magnitude and skip re-assembly when unchanged to avoid unnecessary regeneration. - core: fix typo in Assembly docstring; Mesh now documents register_name; MultiMesh.find_elements documents unsupported name-saving; DataSet.get_data uses fill_unused_nodes=None when converting to preserve source node data. - multimeshdata: add concise docstrings for properties and helpers (active, map, keys, items, shape, ndim) to clarify semantics (operate on active submesh blocks). - util/fdh5: introduce CompressionConfig dataclass; add comprehensive docstrings and parameter validation to FDH5Writer/FDH5Reader methods; preserve submesh ordering and names when reconstructing MultiMesh to avoid silent collisions; remove large example and commented dead code. - util/xdmf_writer: add parameter docs and minor behavior tweak to only export gauss-point fields with a single GP as cell-centered attributes; various code cleanups and formatting. These changes improve usability, documentation, and correctness of IO routines and avoid unnecessary recomputation during assembly. --- environment_doc.yml | 1 + fedoo/constraint/distributed_load.py | 27 +- fedoo/core/assembly.py | 4 +- fedoo/core/dataset.py | 7 +- fedoo/core/mesh.py | 11 +- fedoo/core/multimeshdata.py | 7 + fedoo/util/fdh5.py | 517 +++++++++++++-------------- fedoo/util/xdmf_writer.py | 36 +- 8 files changed, 318 insertions(+), 292 deletions(-) diff --git a/environment_doc.yml b/environment_doc.yml index f211dcb7..96a9478a 100644 --- a/environment_doc.yml +++ b/environment_doc.yml @@ -14,6 +14,7 @@ dependencies: - vtk>=9.2 - pypandoc - pyvista>=0.39 + - meshio - simcoon>=1.13.1 - pytest - sphinx diff --git a/fedoo/constraint/distributed_load.py b/fedoo/constraint/distributed_load.py index 9d20f030..3bacbb4d 100644 --- a/fedoo/constraint/distributed_load.py +++ b/fedoo/constraint/distributed_load.py @@ -36,6 +36,7 @@ def time_func(t_fact): self.time_func = time_func self._assembly_vector = None + self._cached_magnitude = None def str_condensed(self): """Return a condensed one line str describing the object.""" @@ -52,6 +53,7 @@ def initialize(self, problem: ProblemBase): self.assembly.initialize(problem) self._update_during_inc = bool(self.assembly._nlgeom) self._assembly_vector = None + self._cached_magnitude = None self.value = 0 self.start_value = self._start_value_default @@ -94,6 +96,23 @@ def _has_initial_load(self): for attr in ("initial_pressure", "initial_force") ) + def _load_magnitude(self): + # Cached fixed-geometry vectors are only valid while the load magnitude + # is unchanged; generate() compares this snapshot to detect a change. + return tuple( + getattr(self.assembly, attr, None) + for attr in ("pressure", "force", "initial_pressure", "initial_force") + ) + + @staticmethod + def _magnitude_changed(current, cached): + if cached is None: + return True + try: + return not all(np.array_equal(a, b) for a, b in zip(current, cached)) + except (ValueError, TypeError): + return True + def _refresh_reference_values(self, problem: ProblemBase): if self._has_initial_load(): self.start_value = self._assemble_load_vector(problem, 0) @@ -114,8 +133,14 @@ def _format_current_value(self, problem: ProblemBase, value): self._current_value = value def generate(self, problem: ProblemBase, t_fact=1, t_fact_old=None): - if self._update_during_inc or self._assembly_vector is None: + magnitude = self._load_magnitude() + if ( + self._update_during_inc + or self._assembly_vector is None + or self._magnitude_changed(magnitude, self._cached_magnitude) + ): self._refresh_reference_values(problem) + self._cached_magnitude = magnitude self._format_current_value(problem, self.get_value(t_fact, t_fact_old)) return [self] diff --git a/fedoo/core/assembly.py b/fedoo/core/assembly.py index d4d17efb..725ab22b 100644 --- a/fedoo/core/assembly.py +++ b/fedoo/core/assembly.py @@ -779,7 +779,7 @@ def delete_memory(): but it is not done by default. In this case, deleting the memory should resolve the problem. - Note: it the mesh_change argument is set to True when creating the Assembly object, the + Note: if the mesh_change argument is set to True when creating the Assembly object, the memory will be recomputed by default, which may cause a decrease in assembling performances """ Assembly._saved_elementary_operators = {} @@ -834,7 +834,7 @@ def compute_elementary_operators( self._use_local_csys, mesh._local_frame_cache_key(self._element_local_frame), ) - ] = OP # elmRef.computeOperator(nodes,elements) + ] = OP return # ------------------------------------------------------------------- diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index fb9dd8de..205e68e7 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -1300,11 +1300,16 @@ def get_data( if field in self.dict_data[data_type]: data = self.dict_data[data_type][field] else: # if field is not present whith the given data_type search if it exist elsewhere and convert it + # Fetch the source field with all nodes intact: the fill applies + # only to a Node request restricted to the active submesh. When + # converting to Element/GaussPoint data every submesh is used, + # so masking the source nodes here would zero the other + # submeshes' converted values. data, current_data_type = self.get_data( field, component, return_data_type=True, - fill_unused_nodes=fill_unused_nodes, + fill_unused_nodes=None, ) if current_data_type != "Scalar": if self._is_multimesh(): diff --git a/fedoo/core/mesh.py b/fedoo/core/mesh.py index 216e47b5..04f60247 100644 --- a/fedoo/core/mesh.py +++ b/fedoo/core/mesh.py @@ -86,6 +86,10 @@ class Mesh(MeshBase): using ndim = nodes.shape[1] name: str, optional The name of the mesh + register_name: bool, default=True + If True and ``name`` is non-empty, the mesh is registered in the global + mesh dictionary under ``name``. Set to False for transient meshes (e.g. + submeshes) that should not be globally reachable. Example -------- @@ -2030,7 +2034,12 @@ def find_elements( all_nodes: bool = True, name: str | None = None, ): - """Return global element ids matching a selection criterion.""" + """Return global element ids matching a selection criterion. + + Unlike :meth:`Mesh.find_elements`, passing ``name`` to save the result + as an element set is not supported and raises ``NotImplementedError``; + add the set to the relevant submesh instead. + """ element_ids = [] offset = 0 for mesh in self._mesh_list: diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py index bd61799d..efe1f18c 100644 --- a/fedoo/core/multimeshdata.py +++ b/fedoo/core/multimeshdata.py @@ -112,6 +112,7 @@ def _resolve_single(self, selector=None) -> int: @property def active(self): + """Data block of the currently active submesh (see ``active_submesh``).""" return self.submesh(self.active_submesh) def submesh(self, selector): @@ -141,6 +142,8 @@ def submesh(self, selector): ) def map(self, func) -> "MultiMeshData": + """Return a new ``MultiMeshData`` with ``func`` applied to each non-empty + submesh block.""" return MultiMeshData( self.mesh, {i: func(data) for i, data in enumerate(self._data) if data is not None}, @@ -148,9 +151,11 @@ def map(self, func) -> "MultiMeshData": ) def keys(self): + """Submesh ids that hold a (non-empty) data block.""" return [i for i, data in enumerate(self._data) if data is not None] def items(self): + """``(submesh_id, data_block)`` pairs for submeshes that hold data.""" return [(i, data) for i, data in enumerate(self._data) if data is not None] def global_element_location(self, element_id: int) -> tuple[int, int]: @@ -277,10 +282,12 @@ def __getitem__(self, item): @property def shape(self): + """Shape of the active submesh block only (not the whole MultiMesh).""" return np.asarray(self.active).shape @property def ndim(self): + """Number of dimensions of the active submesh block only.""" return np.asarray(self.active).ndim diff --git a/fedoo/util/fdh5.py b/fedoo/util/fdh5.py index 812a3db4..4e078cb6 100644 --- a/fedoo/util/fdh5.py +++ b/fedoo/util/fdh5.py @@ -70,32 +70,31 @@ def iteration_name(iteration: int) -> str: @dataclass(frozen=True) class CompressionConfig: + """HDF5 dataset storage options used by :class:`FDH5Writer`. + + Parameters + ---------- + compression : {"gzip", "lzf"} or None, optional + Compression filter applied to non-scalar datasets; ``None`` disables + compression. Default is ``"gzip"``. + compression_opts : int or None, optional + Compression level for the gzip filter (0-9). Default is 4. + chunks : bool or tuple of int or None, optional + Chunk shape passed to h5py; ``True`` enables auto-chunking and ``None`` + disables chunking. Default is ``True``. + """ + compression: Optional[Literal["gzip", "lzf"]] = "gzip" compression_opts: Optional[int] = 4 chunks: Optional[Union[bool, tuple[int, ...]]] = True class FDH5Writer: - """ - Writer for Finite Element HDF5 files (FDH5). - - This writer implements the following structure: - - mesh/ - nodes - node_sets/ - submesh_X/ - elements - element_sets/ - metadata/ - - results/ - iter_0/ - node_data/ - element_data/submesh_X/ - gausspoint_data/submesh_X/ - scalars/ - metadata/ + """Writer for the Fedoo HDF5 (FDH5) result format. + + Creates or appends to an FDH5 file following the layout described in the + module docstring: a ``mesh/`` group holding the nodes and submeshes, and a + ``results/`` group holding one subgroup per result iteration. """ FILE_VERSION = "1.0" @@ -108,6 +107,20 @@ def __init__( validate: bool = True, create_parents: bool = True, ) -> None: + """Open (or create) an FDH5 file for writing. + + Parameters + ---------- + file_path : str or Path + Path of the FDH5 file to create or append to. + compression : CompressionConfig, optional + HDF5 storage options applied to the written datasets. + validate : bool, optional + If True (default), check array shapes and dtypes before writing. + create_parents : bool, optional + If True (default), create the parent directories of ``file_path`` + when they do not already exist. + """ self.path = Path(file_path) if create_parents: self.path.parent.mkdir(parents=True, exist_ok=True) @@ -121,9 +134,6 @@ def __init__( f.attrs["version"] = self.FILE_VERSION f.attrs["created_utc"] = datetime.now(timezone.utc).isoformat() - # ------------------------------------------------------------------ - # Internal helpers - # ------------------------------------------------------------------ def _create_dataset( self, parent: h5py.Group, @@ -156,9 +166,6 @@ def _next_submesh_id(self, f: h5py.File) -> str: nums = [int(k.split("_")[1]) for k in existing] return f"submesh_{max(nums) + 1}" - # ------------------------------------------------------------------ - # Mesh writing - # ------------------------------------------------------------------ def write_mesh( self, nodes: NDArray, @@ -166,6 +173,18 @@ def write_mesh( node_sets: Optional[Mapping[str, NDArray]] = None, overwrite: bool = False, ) -> None: + """Write the global node coordinates and optional node sets. + + Parameters + ---------- + nodes : numpy.ndarray + Node coordinates with shape ``(n_nodes, dim)``. + node_sets : mapping of str to numpy.ndarray, optional + Named sets of node indices stored under ``mesh/node_sets``. + overwrite : bool, optional + If True, replace existing node data and node sets. If False + (default), writing over existing data raises ``ValueError``. + """ nodes = np.asarray(nodes) if self.validate: @@ -202,6 +221,31 @@ def add_submesh( submesh_id: Optional[str] = None, overwrite: bool = False, ) -> str: + """Add a submesh (a single element type) to the mesh group. + + Parameters + ---------- + element_type : str + Fedoo element type stored in the submesh metadata (e.g. ``"tri3"``). + elements : numpy.ndarray + Connectivity table ``(n_elements, nodes_per_element)`` of 0-based + node indices. + element_sets : mapping of str to numpy.ndarray, optional + Named sets of element indices stored under the submesh. + name : str, optional + Human-readable submesh name stored in the metadata. + submesh_id : str, optional + Target submesh group id. If None (default), the next free + ``submesh_`` id is used. + overwrite : bool, optional + If True, replace an existing submesh with the same id. If False + (default), a collision raises ``ValueError``. + + Returns + ------- + str + The submesh id under which the data was written. + """ elements = np.asarray(elements) if self.validate: @@ -223,7 +267,6 @@ def add_submesh( sm = mesh.create_group(sid) self._create_dataset(sm, "elements", elements) - # Metadata meta = sm.create_group("metadata") meta.attrs["element_type"] = element_type meta.attrs["n_elements"] = elements.shape[0] @@ -231,7 +274,6 @@ def add_submesh( if name: meta.attrs["name"] = name - # Element sets if element_sets: es = sm.create_group("element_sets") for set_name, idx in element_sets.items(): @@ -239,9 +281,6 @@ def add_submesh( return sid - # ------------------------------------------------------------------ - # Iteration writing - # ------------------------------------------------------------------ def write_iteration( self, iteration: int, @@ -254,6 +293,35 @@ def write_iteration( dt: Optional[float] = None, overwrite: bool = False, ) -> str: + """Write all result data for a single iteration. + + Parameters + ---------- + iteration : int + Iteration index; stored as the ``iter_`` results subgroup. + node_data : mapping of str to numpy.ndarray, optional + Nodal fields, each with shape ``(n_nodes, ...)``. + element_data : mapping, optional + Element fields keyed by submesh id, then by field name. + gausspoint_data : mapping, optional + Gauss-point fields keyed by submesh id, then by field name. Each + field uses a GP-major flattened layout of shape + ``(n_elements * n_gauss_points, n_components)``. + scalars : mapping of str to (int, float or numpy.ndarray), optional + Scalar quantities stored under the iteration ``scalars`` group. + time : float, optional + Physical time of the iteration, stored in the metadata. + dt : float, optional + Time increment of the iteration, stored in the metadata. + overwrite : bool, optional + If True, replace an existing iteration. If False (default), a + collision raises ``ValueError``. + + Returns + ------- + str + The HDF5 path of the written iteration group. + """ iter_name = iteration_name(iteration) with h5py.File(self.path, "a") as f: @@ -267,7 +335,6 @@ def write_iteration( it = results.create_group(iter_name) - # Metadata md = it.create_group("metadata") md.attrs["iteration"] = iteration if time is not None: @@ -276,13 +343,11 @@ def write_iteration( md.attrs["dt"] = float(dt) md.attrs["created_utc"] = datetime.now(timezone.utc).isoformat() - # Node data if node_data: nd = it.create_group("node_data") for name, arr in node_data.items(): self._create_dataset(nd, name, np.asarray(arr)) - # Element data if element_data: ed = it.create_group("element_data") for sid, fields in element_data.items(): @@ -290,7 +355,6 @@ def write_iteration( for fname, arr in fields.items(): self._create_dataset(smg, fname, np.asarray(arr)) - # Gauss-point data (GP-major flattened) if gausspoint_data: gd = it.create_group("gausspoint_data") for sid, fields in gausspoint_data.items(): @@ -308,7 +372,6 @@ def write_iteration( else: ds.attrs["n_components"] = arr.shape[1] - # Scalars if scalars: sc = it.create_group("scalars") for name, val in scalars.items(): @@ -333,9 +396,6 @@ def __init__(self, file_path: PathLike) -> None: if not self.path.exists(): raise FileNotFoundError(self.path) - # ------------------------------------------------------------------ - # File handling - # ------------------------------------------------------------------ @contextlib.contextmanager def open(self) -> Iterator[h5py.File]: """ @@ -368,15 +428,28 @@ def _require_lazy_file(self, file: h5py.File | None) -> h5py.File: ) return file - # ------------------------------------------------------------------ - # Mesh - # ------------------------------------------------------------------ def read_nodes( self, *, lazy: bool = False, file: h5py.File | None = None, ) -> Union[NDArray, h5py.Dataset]: + """Read the global node coordinates. + + Parameters + ---------- + lazy : bool, optional + If True, return the open ``h5py.Dataset`` instead of a NumPy array; + this requires ``file`` and that the file remains open while it is + used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + numpy.ndarray or h5py.Dataset + Node coordinates ``(n_nodes, dim)``, eager or lazy. + """ if not lazy: with self._open() as f: return f["mesh/nodes"][...] @@ -390,6 +463,21 @@ def read_node_sets( lazy: bool = False, file: h5py.File | None = None, ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + """Read the named node sets. + + Parameters + ---------- + lazy : bool, optional + If True, return open ``h5py.Dataset`` objects; the file must stay + open while they are used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + dict of str to (numpy.ndarray or h5py.Dataset) + Node index arrays keyed by set name; empty if no node sets exist. + """ if not lazy: out: Dict[str, NDArray] = {} with self._open() as f: @@ -405,6 +493,18 @@ def read_node_sets( return {} if grp is None else {name: ds for name, ds in grp.items()} def list_submeshes(self, *, file: h5py.File | None = None) -> List[str]: + """Return the sorted submesh ids stored in the mesh group. + + Parameters + ---------- + file : h5py.File, optional + Open file handle to reuse; a temporary one is opened if omitted. + + Returns + ------- + list of str + Submesh ids (``submesh_``) in ascending order. + """ if file is None: with self._open() as f: return self.list_submeshes(file=f) @@ -484,10 +584,19 @@ def read_mesh( mesh["submeshes"] = subs return mesh - # ------------------------------------------------------------------ - # Iterations - # ------------------------------------------------------------------ def list_iterations(self, *, file: h5py.File | None = None) -> List[int]: + """Return the sorted result iteration indices stored in the file. + + Parameters + ---------- + file : h5py.File, optional + Open file handle to reuse; a temporary one is opened if omitted. + + Returns + ------- + list of int + Iteration indices in ascending order. + """ if file is None: with self._open() as f: return self.list_iterations(file=f) @@ -500,6 +609,18 @@ def list_iterations(self, *, file: h5py.File | None = None) -> List[int]: ) def read_iteration_metadata(self, iteration: int) -> Dict[str, Any]: + """Read the metadata attributes of one result iteration. + + Parameters + ---------- + iteration : int + Iteration index to read. + + Returns + ------- + dict of str to Any + Metadata attributes (e.g. ``iteration``, ``time``, ``dt``). + """ with self._open() as f: md = f[f"results/{iteration_name(iteration)}/metadata"] return {k: self._decode(v) for k, v in md.attrs.items()} @@ -532,6 +653,23 @@ def read_scalars( lazy: bool = False, file: h5py.File | None = None, ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + """Read the scalar quantities stored for one iteration. + + Parameters + ---------- + iteration : int + Iteration index to read. + lazy : bool, optional + If True, return open ``h5py.Dataset`` objects; the file must stay + open while they are used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + dict of str to (numpy.ndarray or h5py.Dataset) + Scalar values keyed by name; empty if none are stored. + """ if not lazy: out: Dict[str, NDArray] = {} with self._open() as f: @@ -553,6 +691,23 @@ def read_node_data( lazy: bool = False, file: h5py.File | None = None, ) -> Dict[str, Union[NDArray, h5py.Dataset]]: + """Read the nodal fields stored for one iteration. + + Parameters + ---------- + iteration : int + Iteration index to read. + lazy : bool, optional + If True, return open ``h5py.Dataset`` objects; the file must stay + open while they are used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + dict of str to (numpy.ndarray or h5py.Dataset) + Nodal field arrays keyed by field name; empty if none are stored. + """ if not lazy: out: Dict[str, NDArray] = {} with self._open() as f: @@ -574,6 +729,24 @@ def read_element_data( lazy: bool = False, file: h5py.File | None = None, ) -> Dict[str, Dict[str, Union[NDArray, h5py.Dataset]]]: + """Read the element fields stored for one iteration. + + Parameters + ---------- + iteration : int + Iteration index to read. + lazy : bool, optional + If True, return open ``h5py.Dataset`` objects; the file must stay + open while they are used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + dict of str to (dict of str to (numpy.ndarray or h5py.Dataset)) + Element fields keyed by submesh id, then by field name; empty if + none are stored. + """ if not lazy: out: Dict[str, Dict[str, NDArray]] = {} with self._open() as f: @@ -600,6 +773,27 @@ def read_gausspoint_data( lazy: bool = False, file: h5py.File | None = None, ) -> Dict[str, Dict[str, Union[NDArray, h5py.Dataset]]]: + """Read the Gauss-point fields stored for one iteration. + + Each field uses a GP-major flattened layout of shape + ``(n_elements * n_gauss_points, n_components)``. + + Parameters + ---------- + iteration : int + Iteration index to read. + lazy : bool, optional + If True, return open ``h5py.Dataset`` objects; the file must stay + open while they are used. + file : h5py.File, optional + Open file handle, required when ``lazy`` is True. + + Returns + ------- + dict of str to (dict of str to (numpy.ndarray or h5py.Dataset)) + Gauss-point fields keyed by submesh id, then by field name; empty + if none are stored. + """ if not lazy: out: Dict[str, Dict[str, NDArray]] = {} with self._open() as f: @@ -648,18 +842,24 @@ def mesh_to_fedoo(mesh_data: dict): register_name=False, ) - elements_dict = {} - for i, submesh in enumerate(ordered_submeshes): - name = submesh.get("name") or submesh_id(i) - elements_dict[name] = ( - submesh["element_type"], + # Reconstruct submeshes by order, not by name: two submeshes may share the + # same name (e.g. duplicate element types with empty names, stored with + # ``name = elm_type``). Keying a dict by name would silently drop the + # collision, so build ordered Mesh objects and preserve their stored names. + submesh_list = [ + Mesh( + nodes, submesh["elements"], - submesh.get("element_sets", {}), + submesh["element_type"], + element_sets=submesh.get("element_sets", {}), + name=submesh.get("name", ""), + register_name=False, ) + for submesh in ordered_submeshes + ] - return MultiMesh( - nodes, - elements_dict, + return MultiMesh.from_mesh_list( + submesh_list, node_sets=node_sets, register_name=False, ) @@ -808,212 +1008,3 @@ def read_fdh5(filename: str): dataset = MultiFrameDataSet(mesh) dataset.list_data = [("fdh5", str(path), iteration) for iteration in iterations] return dataset - - -# class FDH5File: -# # Class that allow both to write and read data -# # Don't know if this may be usefull -# def __init__(self, path, mode="r"): -# self.path = path -# self.mode = mode - -# if mode == "r": -# self.reader = FDH5Reader(path) -# self.writer = None -# elif mode in ("a", "w"): -# self.writer = FDH5Writer(path) -# self.reader = FDH5Reader(path) -# else: -# raise ValueError("mode must be 'r', 'a', or 'w'") - -# def write_iteration(self, *args, **kwargs): -# if self.writer is None: -# raise RuntimeError("File opened read-only") -# return self.writer.write_iteration(*args, **kwargs) - -# def read_node_data(self, *args, **kwargs): -# return self.reader.read_node_data(*args, **kwargs) - -# def open(self): -# return self.reader.open() - -if __name__ == "__main__": - # example of use - import numpy as np - from pathlib import Path - - # Import your classes - # from fdh5 import FDH5Writer, FDH5Reader - - # -------------------------------------------------- - # File path - # -------------------------------------------------- - path = Path("example.fdh5") - - writer = FDH5Writer(path, validate=True) - - # -------------------------------------------------- - # Mesh definition - # -------------------------------------------------- - - # 4 nodes, 2D - nodes = np.array( - [ - [0.0, 0.0], - [1.0, 0.0], - [1.0, 1.0], - [0.0, 1.0], - ], - dtype=float, - ) - - # Write mesh + node sets - writer.write_mesh( - nodes, - node_sets={ - "boundary": np.array([0, 1, 2, 3], dtype=int), - }, - overwrite=True, - ) - - # Two TRI3 elements - elements = np.array( - [ - [0, 1, 2], - [0, 2, 3], - ], - dtype=int, - ) - - submesh_id = writer.add_submesh( - element_type="tri3", - elements=elements, - name="square_triangles", - ) - - # -------------------------------------------------- - # Iteration 0 - # -------------------------------------------------- - - # Node field: displacement (n_nodes, 2) - node_data = { - "displacement": np.array( - [ - [0.0, 0.0], - [0.1, 0.0], - [0.1, 0.1], - [0.0, 0.1], - ], - dtype=float, - ) - } - - # Element field: von Mises stress (n_elements,) - element_data = { - submesh_id: { - "stress_vm": np.array([100.0, 120.0], dtype=float), - } - } - - # Gauss‑point field (GP‑major, flattened) - # Here: 2 elements, 2 Gauss points each - # Order: - # GP0: elem0, elem1 - # GP1: elem0, elem1 - gausspoint_data = { - submesh_id: { - "strain_eq": np.array( - [ - 0.01, - 0.02, # GP0 - 0.015, - 0.025, # GP1 - ], - dtype=float, - ).reshape(-1, 1) # (n_elem * n_gp, n_comp) - } - } - - # Scalars - scalars = { - "time": 0.0, - "total_energy": 42.0, - } - - # Write iteration - writer.write_iteration( - iteration=0, - node_data=node_data, - element_data=element_data, - gausspoint_data=gausspoint_data, - scalars=scalars, - time=0.0, - dt=0.1, - ) - - print("✅ File written:", path) - - # ------------------------------------------------------------------------- - # Read written file - # ------------------------------------------------------------------------- - reader = FDH5Reader(path) - - # Read mesh - mesh = reader.read_mesh() - print("Nodes:\n", mesh["nodes"]) - print("Submeshes:", list(mesh["submeshes"].keys())) - - # Read iteration 0 - it0 = reader.read_iteration(0) - - u = it0["node_data"]["displacement"] - stress = it0["element_data"]["submesh_0"]["stress_vm"] - strain_gp = it0["gausspoint_data"]["submesh_0"]["strain_eq"] - - print("\nDisplacement:\n", u) - print("\nElement stress:\n", stress) - print("\nGauss-point strain (flattened):\n", strain_gp) - - # ------------------------------------------------------------------------- - # Read written file lazily - # ------------------------------------------------------------------------- - reader = FDH5Reader(path) - - with reader.open() as f: - # Lazy node data - node_data = reader.read_node_data(0, lazy=True, file=f) - u_ds = node_data["displacement"] # h5py.Dataset - - # Only read first 2 nodes - print("First two displacements:\n", u_ds[:2]) - - # Lazy element data - elem_data = reader.read_element_data(0, lazy=True, file=f) - stress_ds = elem_data["submesh_0"]["stress_vm"] - - print("First element stress:", stress_ds[0]) - - # Lazy Gauss-point data - gp_data = reader.read_gausspoint_data(0, lazy=True, file=f) - eps_ds = gp_data["submesh_0"]["strain_eq"] - - # Infer GP layout manually - n_elems = f["mesh/submesh_0/elements"].shape[0] - - # GP0 for all elements - gp0 = eps_ds[0:n_elems] - print("GP0 strain:\n", gp0) - - # ------------------------------------------------------------------------- - # Read written file lazily - # ------------------------------------------------------------------------- - - from xdmf_writer import XDMFExporter # your exporter class - - exporter = XDMFExporter(Path("example.fdh5")) - exporter.export() - - import pyvista as pv - - mesh = pv.read("example.xdmf") - mesh.plot(show_edges=True) diff --git a/fedoo/util/xdmf_writer.py b/fedoo/util/xdmf_writer.py index 90bdd890..5d241f64 100644 --- a/fedoo/util/xdmf_writer.py +++ b/fedoo/util/xdmf_writer.py @@ -34,12 +34,18 @@ def __init__( h5_path: Path, xdmf_path: Optional[Path] = None, ) -> None: + """ + Parameters + ---------- + h5_path : Path + Path to the source FDH5 file to describe. + xdmf_path : Path, optional + Output path for the ``.xdmf`` file. Defaults to ``h5_path`` with a + ``.xdmf`` suffix. + """ self.h5_path = Path(h5_path) self.xdmf_path = xdmf_path or self.h5_path.with_suffix(".xdmf") - # ------------------------------------------------------------------ - # Public API - # ------------------------------------------------------------------ def export(self) -> Path: """ Generate the XDMF file for all iterations (time series). @@ -69,9 +75,6 @@ def export(self) -> Path: return self.xdmf_path - # ------------------------------------------------------------------ - # Iterations - # ------------------------------------------------------------------ def _list_iterations(self, f: h5py.File) -> Iterable[int]: grp = f.get("results") if grp is None: @@ -99,9 +102,6 @@ def _build_iteration_grid(self, f: h5py.File, iteration: int) -> etree.Element: return it_grid - # ------------------------------------------------------------------ - # Submeshes - # ------------------------------------------------------------------ def _list_submeshes(self, f: h5py.File) -> Iterable[str]: mesh = f.get("mesh") if mesh is None: @@ -131,9 +131,7 @@ def _build_submesh_grid( grid = etree.Element("Grid", Name=submesh_id, GridType="Uniform") - # ------------------------------ # Topology - # ------------------------------ topo = etree.SubElement( grid, "Topology", @@ -149,9 +147,7 @@ def _build_submesh_grid( Format="HDF", ).text = f"{self.h5_path.name}:/mesh/{submesh_id}/elements" - # ------------------------------ # Geometry (global nodes) - # ------------------------------ nodes = f["mesh/nodes"] n_nodes, dim = nodes.shape geom_type = "XYZ" if dim == 3 else "XY" @@ -171,9 +167,7 @@ def _build_submesh_grid( Format="HDF", ).text = f"{self.h5_path.name}:/mesh/nodes" - # ------------------------------ # Node data (global) - # ------------------------------ node_grp = it_grp.get("node_data") if node_grp is not None: for name, ds in node_grp.items(): @@ -185,9 +179,7 @@ def _build_submesh_grid( shape=ds.shape, ) - # ------------------------------ # Element data (per submesh) - # ------------------------------ elem_grp = it_grp.get(f"element_data/{submesh_id}") if elem_grp is not None: for name, ds in elem_grp.items(): @@ -202,13 +194,12 @@ def _build_submesh_grid( shape=ds.shape, ) - # ------------------------------ - # Gauss-point data (flattened GP-major) - # Exported as Cell-centered attributes - # ------------------------------ gp_grp = it_grp.get(f"gausspoint_data/{submesh_id}") if gp_grp is not None: for name, ds in gp_grp.items(): + # Only single-Gauss-point fields map onto a cell-centered + # attribute; multi-GP fields have no per-cell value and are + # skipped (XDMF has no native multi-value-per-cell attribute). n_gauss_points = int(ds.attrs.get("n_gauss_points", 1)) if n_gauss_points != 1: continue @@ -225,9 +216,6 @@ def _build_submesh_grid( return grid - # ------------------------------------------------------------------ - # Attributes - # ------------------------------------------------------------------ def _add_attribute( self, grid: etree.Element, From a5db6a142815b4fa49e21ab47d54b9fab2a830a5 Mon Sep 17 00:00:00 2001 From: pruliere Date: Fri, 10 Jul 2026 09:59:45 +0200 Subject: [PATCH 12/14] update doc files and improve the viewer plotting for multi-mesh data --- fedoo/core/dataset.py | 12 +++++ fedoo/core/multimeshdata.py | 8 ++- fedoo/post_processing/__init__.py | 69 ++++++++++++++++++++++--- fedoo/util/viewer.py | 86 +++++++++++++++++++------------ 4 files changed, 135 insertions(+), 40 deletions(-) diff --git a/fedoo/core/dataset.py b/fedoo/core/dataset.py index 205e68e7..4cbf4a2b 100644 --- a/fedoo/core/dataset.py +++ b/fedoo/core/dataset.py @@ -1282,6 +1282,18 @@ def get_data( With a ``MultiMesh``, element and Gauss point fields may be stored as a dictionary keyed by submesh id, submesh name, or unique element type. The active submesh is controlled by ``dataset.active_submesh``. + + Element ids in a ``MultiMesh`` use a global, concatenated numbering: + all elements of submesh 0, then all elements of submesh 1, and so on. + Retrieve values using those ids from the returned ``MultiMeshData``:: + + stress = dataset.get_data("Stress", "vm", "Element") + value = stress.global_element_value(global_element_id) + values = stress.global_element_values(global_element_ids) + + Use ``stress.to_global()`` when a single NumPy array in this global + order is needed. ``np.asarray(stress)`` instead returns the data of + the active submesh. """ if data_type is None: # search if field exist somewhere diff --git a/fedoo/core/multimeshdata.py b/fedoo/core/multimeshdata.py index efe1f18c..432ae897 100644 --- a/fedoo/core/multimeshdata.py +++ b/fedoo/core/multimeshdata.py @@ -24,6 +24,8 @@ class MultiMeshData: >>> data.active # data for dataset.active_submesh >>> data.submesh(1) # data for submesh id 1 >>> data.submesh("tri3") # data for all tri3 submeshes + >>> data.global_element_value(42) # value for global element 42 + >>> data.global_element_values([3, 18, 42]) Parameters ---------- @@ -41,7 +43,11 @@ class MultiMeshData: ----- Missing submesh data are stored as ``None``. ``to_global`` concatenates selected submesh blocks in mesh order and fills missing blocks with - ``fill_value``. + ``fill_value``. Global element ids use that same concatenated order. + + NumPy-style access, including ``np.asarray(data)`` and ``data[...]``, + intentionally applies to the active submesh only. Use ``to_global()`` or + ``global_element_value(s)`` when working with global element ids. """ def __init__(self, mesh: MultiMesh, data, active_submesh=0) -> None: diff --git a/fedoo/post_processing/__init__.py b/fedoo/post_processing/__init__.py index 5d7a1a8a..f805238e 100644 --- a/fedoo/post_processing/__init__.py +++ b/fedoo/post_processing/__init__.py @@ -12,7 +12,9 @@ :py:meth:`fedoo.Problem.get_results` method of the problem class. The get_results method returns a :py:class:`fedoo.DataSet` object which comes -with several methods for plotting, saving and loading mesh dependent results. +with several methods for plotting, saving and loading mesh-dependent results. +A dataset can be associated with either a :py:class:`fedoo.Mesh` or a +:py:class:`fedoo.MultiMesh`. To avoid a redundent call of the get_results function, especially for time @@ -26,8 +28,8 @@ solved using :py:meth:`fedoo.Problem.nlsolve`, results are automatically saved at some iterations dependending on the choosen parameters. -The :py:class:`MultiFrameDataSet` store the path of the saved files for each -iteration in the MultiFrameDataSet.list_data attribute. The method +The :py:class:`MultiFrameDataSet` stores references to the saved iterations in +the MultiFrameDataSet.list_data attribute. The method :py:meth:`MultiFrameDataSet.load` is called to read the data of a given iteration. @@ -50,6 +52,51 @@ MultiFrameDataSet +MultiMesh results +----------------- + +A :py:class:`DataSet` or :py:class:`MultiFrameDataSet` associated with a +:py:class:`MultiMesh` can contain submeshes with different element types. +Node fields use the common global node numbering. Element and Gauss-point +fields returned by :py:meth:`DataSet.get_data` are stored in a +:py:class:`MultiMeshData` object, with one data block per submesh. + +Global element ids use the concatenated submesh order: all elements of +submesh 0, then all elements of submesh 1, and so on. For example: + +.. code-block:: python + + stress = results.get_data("Stress", "vm", "Element") + + # Access data by submesh or by global element id + shell_stress = stress.submesh("tri3") + value = stress.global_element_value(42) + values = stress.global_element_values([3, 18, 42]) + + # Build one NumPy array using the global element order + global_stress = stress.to_global() + +Normal NumPy-style access to a :py:class:`MultiMeshData` object refers to its +active submesh. Use :py:meth:`MultiMeshData.to_global` when a globally ordered +array is required. A global selection spanning several submeshes can be +plotted directly: + +.. code-block:: python + + results.plot( + "Stress", + component="vm", + data_type="Element", + element_set=[3, 18, 42], + global_element_set=True, + ) + +.. autosummary:: + :toctree: generated/ + :template: custom-class-template.rst + + MultiMeshData + Save data to disk ----------------- @@ -57,8 +104,16 @@ :py:meth:`fedoo.Problem.get_results` method, the data can easily be saved on disk using for instance the :py:meth:`fedoo.DataSet.save` method. +FDH5 is the default and recommended format. If ``DataSet.save``, +``MultiFrameDataSet.save_all``, or ``read_data`` is called with a filename +without an extension, the ``.fdh5`` extension is assumed. + The available file types are: - * 'fdz': A zipped archive containing the mesh using the 'vtk' format + * 'fdh5': The native HDF5-based Fedoo format. It stores the mesh, node and + element sets, node/element/Gauss-point fields, and multiple iterations + in one file. It also preserves the complete MultiMesh structure and its + per-submesh fields. This is the default format. + * 'fdz': A legacy zipped archive containing the mesh using the 'vtk' format named '_mesh_.vtk', and data from several iterations named 'iter_x.npz' where x is the iteration number (x=0 for the 1st iteration). * 'vtk': The vtk format contains the mesh and the data in a single files. @@ -86,6 +141,7 @@ To read data saved on disk, use the function :py:func:`read_data`. The data are imported as :py:class:`DataSet` or :py:class:`MultiFrameDataSet` objects depending on the imported file(s). +FDH5 is used by default when the filename has no extension. .. _very_simple_problem: @@ -169,11 +225,11 @@ .. code-block:: python import fedoo as fd - result = fd.read_data('myfile.fdz') # load a DataSet from file + result = fd.read_data('myfile.fdh5') # load a DataSet from file fd.viewer() # start the viewer with no file opened fd.viewer(result) # start the viewer and open the result DataSet - fd.viewer('myfile.fdz') # start the viewer with the data from a file + fd.viewer('myfile.fdh5') # start the viewer with the data from a file The viewer includes the following tools and features: @@ -233,6 +289,7 @@ DataSet.to_excel DataSet.to_vtk DataSet.to_msh + DataSet.to_fdh5 MultiFrameDataSet.save_all diff --git a/fedoo/util/viewer.py b/fedoo/util/viewer.py index 5f735e6b..40ee5b2d 100644 --- a/fedoo/util/viewer.py +++ b/fedoo/util/viewer.py @@ -1,6 +1,7 @@ import fedoo as fd import numpy as np import sys +from contextlib import contextmanager from qtpy import QtWidgets, QtGui from qtpy.QtWidgets import ( QDockWidget, @@ -46,6 +47,23 @@ def _data_n_iter(data): return getattr(data, "n_iter", 1) +@contextmanager +def _defer_rendering(plotter, enabled): + """Render a rebuilt scene only after all its actors are ready.""" + if not enabled or not hasattr(plotter, "suppress_rendering"): + yield + return + + was_suppressed = plotter.suppress_rendering + plotter.suppress_rendering = True + try: + yield + finally: + plotter.suppress_rendering = was_suppressed + if not was_suppressed: + plotter.render() + + def _as_3d_tuple(values, fill=0.0): if values is None: return (fill, fill, fill) @@ -379,45 +397,47 @@ def update_plot(self, val=None, iteration=None, lock_view=True, plotter=None): # "n_colors": self.opts['n_colors'], "n_labels": self.opts["n_labels"], } - # plotter.clear() # not compatible with pbr ??? - plotter.renderer.clear_actors() + is_multimesh = _is_multimesh(self.data.mesh) multimesh_kargs = {} - if _is_multimesh(self.data.mesh): + if is_multimesh: multimesh_kargs["global_element_set"] = True multimesh_kargs["name"] = "data" - self.data.plot( - field=self.current_field, - component=self.current_comp, - data_type=self.current_data_type, - clim=clim, - plotter=plotter, - show=False, - show_edges=self.opts["show_edges"], - show_scalar_bar=self.opts["show_scalar_bar"], - scalar_bar_args=sargs, - node_labels=self.opts["node_labels"], - element_labels=self.opts["element_labels"], - scale=self.opts["scale"], - opacity=self.opts["opacity"], - pbr=self.opts["pbr"], - metallic=self.opts["metallic"], - roughness=self.opts["roughness"], - diffuse=self.opts["diffuse"], - clip_args=self.opts["clip_args"], - lock_view=lock_view, - title=self.opts["title_plot"], - element_set=self.opts["element_set"], - # element_set_invert = self.opts["element_set_invert"], - cmap=self.opts["cmap"], - **multimesh_kargs, - ) + with _defer_rendering(plotter, is_multimesh): + # plotter.clear() # not compatible with pbr ??? + plotter.renderer.clear_actors() + self.data.plot( + field=self.current_field, + component=self.current_comp, + data_type=self.current_data_type, + clim=clim, + plotter=plotter, + show=False, + show_edges=self.opts["show_edges"], + show_scalar_bar=self.opts["show_scalar_bar"], + scalar_bar_args=sargs, + node_labels=self.opts["node_labels"], + element_labels=self.opts["element_labels"], + scale=self.opts["scale"], + opacity=self.opts["opacity"], + pbr=self.opts["pbr"], + metallic=self.opts["metallic"], + roughness=self.opts["roughness"], + diffuse=self.opts["diffuse"], + clip_args=self.opts["clip_args"], + lock_view=lock_view, + title=self.opts["title_plot"], + element_set=self.opts["element_set"], + # element_set_invert = self.opts["element_set_invert"], + cmap=self.opts["cmap"], + **multimesh_kargs, + ) - if self.parent()._plane_widget_enabled: - self.parent().enable_plane_widget() + if self.parent()._plane_widget_enabled: + self.parent().enable_plane_widget() - if self.parent()._line_widget_enabled: - self.parent()._rebuild_line_widget() + if self.parent()._line_widget_enabled: + self.parent()._rebuild_line_widget() def closeEvent(self, event): self.plotter.close() From e36417de260669a180ba82108e8b75b9ae3b48f3 Mon Sep 17 00:00:00 2001 From: pruliere Date: Fri, 10 Jul 2026 12:28:21 +0200 Subject: [PATCH 13/14] Correct an update error for Shell Add a small launch test feature at the end of all tests --- fedoo/constitutivelaw/shell.py | 39 +++++++++++----------- fedoo/weakform/plate.py | 21 ------------ tests/test_BeamElement.py | 4 +-- tests/test_PlateElement.py | 1 + tests/test_axi_phase3.py | 4 +++ tests/test_axi_simcoon_structural.py | 4 +++ tests/test_axi_to_3d_rotation.py | 4 +++ tests/test_cantilever_beam_3D_model.py | 6 ++++ tests/test_distributed_load_neumann_bc.py | 4 +++ tests/test_factor_reuse.py | 4 +++ tests/test_octet.py | 3 +- tests/test_periodicity.py | 4 +++ tests/test_platewithhol.py | 6 ++++ tests/test_poromech_mandel.py | 5 +-- tests/test_poromech_mandel_analytical.py | 5 +-- tests/test_poromech_taylor_hood.py | 5 +-- tests/test_poromech_terzaghi.py | 6 ++-- tests/test_poromech_terzaghi_analytical.py | 5 +-- tests/test_recovery.py | 4 +++ tests/test_simcoon_temperature.py | 6 ++++ tests/test_stress_equilibrium_mixed.py | 6 ++-- tests/test_thermal3D.py | 6 ++++ 22 files changed, 95 insertions(+), 57 deletions(-) diff --git a/fedoo/constitutivelaw/shell.py b/fedoo/constitutivelaw/shell.py index 4c72f3f2..04c35b34 100644 --- a/fedoo/constitutivelaw/shell.py +++ b/fedoo/constitutivelaw/shell.py @@ -60,27 +60,28 @@ def _material_density(material, owner_name): return density def update(self, assembly, pb): - # disp = pb.get_disp() - # rot = pb.get_rot() - U = pb.get_dof_solution() - if np.isscalar(U) and U == 0: - assembly.sv["ShellStrain"] = 0 + """Update shell stress from the strain computed by the weak form.""" + shell_strain = assembly.sv["ShellStrain"] + if np.isscalar(shell_strain) and shell_strain == 0: assembly.sv["ShellStress"] = 0 - else: - ShellStrainOp = assembly.weakform.generalized_strain_operator() - ShellStrain = [ - op if np.isscalar(op) else assembly.get_gp_results(op, U) - for op in ShellStrainOp - ] - - H = self.get_shell_stiffness_matrix() + return - # all terms are computed. Perhaps could be optimized by computed only the termes related to the associated weak form (eg shear for selective integration) - assembly.sv["ShellStress"] = [ - sum([ShellStrain[j] * assembly.convert_data(H[i][j]) for j in range(8)]) - for i in range(8) - ] # H[i][j] are converted to gauss point excepted if scalar - assembly.sv["ShellStrain"] = ShellStrain + H = self.get_shell_stiffness_matrix() + assembly.sv["_ShellStiffnessMatrix"] = H + stress = [ + sum( + [ + ( + shell_strain[j] * assembly.convert_data(H[i][j]) + if not np.array_equal(shell_strain[j], 0) + else 0 + ) + for j in range(8) + ] + ) + for i in range(8) + ] + assembly.sv["ShellStress"] = type(shell_strain)(stress) def get_strain(self, assembly, **kargs): """Return the last computed strain associated to the given assembly. diff --git a/fedoo/weakform/plate.py b/fedoo/weakform/plate.py index 86beb1f1..de67d2d9 100644 --- a/fedoo/weakform/plate.py +++ b/fedoo/weakform/plate.py @@ -302,16 +302,12 @@ def update(self, assembly, pb): if np.array_equal(pb.get_dof_solution(), 0): assembly.sv["ShellStrain"] = 0 - assembly.sv["ShellStress"] = 0 assembly.sv["_DrillConstraint"] = 0 return op_plate_strain = self.generalized_strain_operator() op_drill_constraint = self.drill_constraint_operator() - H = self.constitutivelaw.get_shell_stiffness_matrix() - assembly.sv["_ShellStiffnessMatrix"] = H - if self.nlgeom: use_local_dof = True dof = self._compute_local_dof(assembly, pb) @@ -340,23 +336,6 @@ def update(self, assembly, pb): op_drill_constraint, dof, use_local_dof=use_local_dof ) - # Evaluate Stresses (Shared linear/nonlinear) - assembly.sv["ShellStress"] = _ShellComponentList( - [ - sum( - [ - ( - assembly.sv["ShellStrain"][j] * H[i][j] - if not np.array_equal(assembly.sv["ShellStrain"][j], 0) - else 0 - ) - for j in range(8) - ] - ) - for i in range(8) - ] - ) - def to_start(self, assembly, pb): if self.nlgeom == "UL": assembly.set_disp(pb.get_disp()) diff --git a/tests/test_BeamElement.py b/tests/test_BeamElement.py index 2489530b..09ad4b78 100644 --- a/tests/test_BeamElement.py +++ b/tests/test_BeamElement.py @@ -3,9 +3,7 @@ # import numpy as np - import fedoo as fd -import pytest def test_beam_element(): @@ -109,4 +107,6 @@ def test_beam_element(): if __name__ == "__main__": + import pytest + pytest.main([__file__]) diff --git a/tests/test_PlateElement.py b/tests/test_PlateElement.py index 29451833..b54ba770 100644 --- a/tests/test_PlateElement.py +++ b/tests/test_PlateElement.py @@ -4,6 +4,7 @@ import numpy as np import pytest + import fedoo as fd # Define the combinations to test: (geom_type, element_formulation, test_nlgeom, ref_sol) diff --git a/tests/test_axi_phase3.py b/tests/test_axi_phase3.py index 4da2070a..f3dfdb88 100644 --- a/tests/test_axi_phase3.py +++ b/tests/test_axi_phase3.py @@ -76,3 +76,7 @@ def test_rotary_inertia_axi_raises(): fd.ModelingSpace("2Daxi") with pytest.raises(NotImplementedError, match="2Daxi"): fd.weakform.RotaryInertia(1.0) + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_axi_simcoon_structural.py b/tests/test_axi_simcoon_structural.py index 133a8933..c8fad3d9 100644 --- a/tests/test_axi_simcoon_structural.py +++ b/tests/test_axi_simcoon_structural.py @@ -245,3 +245,7 @@ def test_set_start_rotation_preserves_axi_zeros_for_pure_meridional_F(): assert DR0[0, 2] == 0.0 and DR0[1, 2] == 0.0 assert DR0[2, 0] == 0.0 and DR0[2, 1] == 0.0 assert DR0[2, 2] == pytest.approx(1.0, abs=1e-14) + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_axi_to_3d_rotation.py b/tests/test_axi_to_3d_rotation.py index e0207900..f0310c4d 100644 --- a/tests/test_axi_to_3d_rotation.py +++ b/tests/test_axi_to_3d_rotation.py @@ -205,3 +205,7 @@ def test_invalid_kind_raises(): theta = np.array([0.0]) with pytest.raises(ValueError, match="kind"): _revolve_axi_field(field, theta, n_nodes, kind="bogus") + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_cantilever_beam_3D_model.py b/tests/test_cantilever_beam_3D_model.py index 440407dc..d78f6986 100644 --- a/tests/test_cantilever_beam_3D_model.py +++ b/tests/test_cantilever_beam_3D_model.py @@ -70,3 +70,9 @@ def test_cantilever_beam_3D_model(): ) assert np.abs(TensorStress[5][-1] + 0.9007983467254552) < 1e-10 + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_distributed_load_neumann_bc.py b/tests/test_distributed_load_neumann_bc.py index d81c723f..28eb673d 100644 --- a/tests/test_distributed_load_neumann_bc.py +++ b/tests/test_distributed_load_neumann_bc.py @@ -247,3 +247,7 @@ def test_assembly_neumann_load_is_ignored_by_standard_start_value_update(): pb.set_D(0) pb.init_bc_start_value() + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_factor_reuse.py b/tests/test_factor_reuse.py index cefdae0c..38f21a5d 100644 --- a/tests/test_factor_reuse.py +++ b/tests/test_factor_reuse.py @@ -127,3 +127,7 @@ def test_solver_mumps_default_path(): x = _solver_mumps(A, b) assert np.allclose(A @ x, b, atol=1e-10) + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_octet.py b/tests/test_octet.py index b31364d3..d011ca51 100644 --- a/tests/test_octet.py +++ b/tests/test_octet.py @@ -2,7 +2,6 @@ import numpy as np import fedoo as fd -import pytest def test_octet(): @@ -82,4 +81,6 @@ def test_octet(): if __name__ == "__main__": + import pytest + pytest.main([__file__]) diff --git a/tests/test_periodicity.py b/tests/test_periodicity.py index dcfde00e..9fed0f29 100644 --- a/tests/test_periodicity.py +++ b/tests/test_periodicity.py @@ -97,3 +97,7 @@ def test_periodic_bc_works_on_clean_box_at_default_tol(): bc = fd.constraint.PeriodicBC("small_strain") # default tol=1e-8 pb.bc.add(bc) pb.apply_boundary_conditions() # initializes the BC list + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_platewithhol.py b/tests/test_platewithhol.py index 5a2c3130..a71d0805 100644 --- a/tests/test_platewithhol.py +++ b/tests/test_platewithhol.py @@ -55,3 +55,9 @@ def test_plate_with_hole(): # res.plot('Stress_vm') # fd.util.field_plot_2d("Assembly", disp = pb.get_disp(), dataname = 'stress', component='vm', data_min=None, data_max = None, scale_factor = 6, plot_edge = True, nb_level = 10, type_plot = "smooth") # fd.util.field_plot_2d("Assembly", disp = pb.get_disp(), dataname = 'disp', component=0, scale_factor = 6, plot_edge = True, nb_level = 6, type_plot = "smooth") + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_poromech_mandel.py b/tests/test_poromech_mandel.py index 69c45fe3..cdea8793 100644 --- a/tests/test_poromech_mandel.py +++ b/tests/test_poromech_mandel.py @@ -151,5 +151,6 @@ def test_mandel_cryer_overshoot(): if __name__ == "__main__": - test_mandel_cryer_overshoot() - print("test_poromech_mandel: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_poromech_mandel_analytical.py b/tests/test_poromech_mandel_analytical.py index 53b7732d..890e2d48 100644 --- a/tests/test_poromech_mandel_analytical.py +++ b/tests/test_poromech_mandel_analytical.py @@ -165,5 +165,6 @@ def p_centre(t): if __name__ == "__main__": - test_mandel_vs_abousleiman() - print("test_poromech_mandel_analytical: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_poromech_taylor_hood.py b/tests/test_poromech_taylor_hood.py index e6764a40..808bc43e 100644 --- a/tests/test_poromech_taylor_hood.py +++ b/tests/test_poromech_taylor_hood.py @@ -99,5 +99,6 @@ def test_taylor_hood_poro_undrained_magnitude(): if __name__ == "__main__": - test_taylor_hood_poro_undrained_magnitude() - print("test_poromech_taylor_hood: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_poromech_terzaghi.py b/tests/test_poromech_terzaghi.py index 733c22f0..2aade22d 100644 --- a/tests/test_poromech_terzaghi.py +++ b/tests/test_poromech_terzaghi.py @@ -16,7 +16,6 @@ """ import numpy as np -import pytest import fedoo as fd @@ -179,5 +178,6 @@ def test_terzaghi_consolidation_smoke(): if __name__ == "__main__": - test_terzaghi_consolidation_smoke() - print("test_poromech_terzaghi: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_poromech_terzaghi_analytical.py b/tests/test_poromech_terzaghi_analytical.py index 436638c9..a8c8bc35 100644 --- a/tests/test_poromech_terzaghi_analytical.py +++ b/tests/test_poromech_terzaghi_analytical.py @@ -142,5 +142,6 @@ def test_terzaghi_consolidation_vs_analytical(): if __name__ == "__main__": - test_terzaghi_consolidation_vs_analytical() - print("test_poromech_terzaghi_analytical: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_recovery.py b/tests/test_recovery.py index c4289827..21e26d35 100644 --- a/tests/test_recovery.py +++ b/tests/test_recovery.py @@ -266,3 +266,7 @@ def test_recover_hessian_validates_shape(): mesh = fd.mesh.rectangle_mesh(nx=5, ny=5, elm_type="quad4") with pytest.raises(ValueError, match="shape"): fd.recover_hessian(mesh, np.zeros(mesh.n_nodes - 1)) + + +if __name__ == "__main__": + pytest.main([__file__]) diff --git a/tests/test_simcoon_temperature.py b/tests/test_simcoon_temperature.py index b7d2acfb..b2ed4876 100644 --- a/tests/test_simcoon_temperature.py +++ b/tests/test_simcoon_temperature.py @@ -62,3 +62,9 @@ def test_scalar_zero_temp_is_none(): space, mesh, law, assembly, pb = _build(register_temp=False, tag="zero") assembly.sv["Temp"] = 0 assert law.get_temp_gp(assembly, pb) is None + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) diff --git a/tests/test_stress_equilibrium_mixed.py b/tests/test_stress_equilibrium_mixed.py index e396e44a..0b22dea7 100644 --- a/tests/test_stress_equilibrium_mixed.py +++ b/tests/test_stress_equilibrium_mixed.py @@ -102,6 +102,6 @@ def test_mixed_bulk_none_matches_explicit_bulk(): if __name__ == "__main__": - test_mixed_nearly_incompressible_default_bulk() - test_mixed_bulk_none_matches_explicit_bulk() - print("test_stress_equilibrium_mixed: OK") + import pytest + + pytest.main([__file__]) diff --git a/tests/test_thermal3D.py b/tests/test_thermal3D.py index 7d8fa760..5f3ad483 100644 --- a/tests/test_thermal3D.py +++ b/tests/test_thermal3D.py @@ -90,3 +90,9 @@ def test_thermal3D(): assert np.abs(legacy_temp[8712] - 2.610859332847924) < 1e-8 np.testing.assert_allclose(new_temp, legacy_temp, rtol=1e-10, atol=1e-12) + + +if __name__ == "__main__": + import pytest + + pytest.main([__file__]) From 5b50769fff9edb5ffe075ef99721b3083465ef03 Mon Sep 17 00:00:00 2001 From: pruliere Date: Fri, 10 Jul 2026 16:14:49 +0200 Subject: [PATCH 14/14] - add new dynamic example for multimesh (shell+solid) + dynamic +ipc_contact - small bug corrections --- .../dynamique/multi_mesh_dynamic_contact.py | 208 ++++++++++++++++++ fedoo/core/problem.py | 1 + fedoo/time/common.py | 15 ++ fedoo/weakform/damping_stabilization.py | 16 +- 4 files changed, 238 insertions(+), 2 deletions(-) create mode 100644 examples/dynamique/multi_mesh_dynamic_contact.py diff --git a/examples/dynamique/multi_mesh_dynamic_contact.py b/examples/dynamique/multi_mesh_dynamic_contact.py new file mode 100644 index 00000000..da0da457 --- /dev/null +++ b/examples/dynamique/multi_mesh_dynamic_contact.py @@ -0,0 +1,208 @@ +""" +Dynamic MultiMesh example: shell ping-pong ball dropped on a solid plate. + +The example combines: + - a shell sphere meshed with ``tri3`` elements, + - a solid plate meshed with ``hex8`` elements, + - a ``MultiMesh`` obtained with ``mesh_sphere + mesh_plate``, + - gravity as a distributed load on the shell sphere, + - IPC contact on the sphere and the top surface of the plate, + - FDH5 time output from the structural ``AssemblySum``. + +The discretisation is intentionally modest so the file stays usable as an +example. Increase the sphere resolution, plate mesh density, and time horizon +for production-quality results. +""" + +import os + +import fedoo as fd +import numpy as np +import pyvista as pv + + +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +RESULTS_DIR = os.path.join(SCRIPT_DIR, "results") +os.makedirs(RESULTS_DIR, exist_ok=True) + + +############################################################################### +# Parameters +# + +# Units: N, mm, tonne, s, MPa = N/mm2 +E = 2.0e3 +NU = 0.37 +DENSITY = 10.0e-9 # tonne / mm3 +GRAVITY = 9.81e3 # mm / s2 + +RADIUS = 20.0 +THICKNESS = 0.45 +GAP = 30 # 0.5 + +PLATE_X = 100.0 +PLATE_Y = 100.0 +PLATE_H = 2.0 + +SPHERE_CENTER = np.array([0.5 * PLATE_X, 0.5 * PLATE_Y, PLATE_H + RADIUS + GAP]) + + +############################################################################### +# Mesh +# + +sphere = pv.Sphere( + radius=RADIUS, + center=SPHERE_CENTER, + theta_resolution=24, + phi_resolution=24, +) +sphere_mesh = fd.Mesh.from_pyvista(sphere) +sphere_mesh.element_sets["sphere"] = np.arange(sphere_mesh.n_elements) + +plate_mesh = fd.mesh.box_mesh( + 20, + 20, + 3, + x_max=PLATE_X, + y_max=PLATE_Y, + z_max=PLATE_H, +) +plate_mesh.element_sets["plate"] = np.arange(plate_mesh.n_elements) + +mesh = sphere_mesh + plate_mesh +print(f"MultiMesh: {mesh.n_nodes} nodes, {mesh.n_elements} elements, {mesh.elm_types}") + + +############################################################################### +# Models and structural assemblies +# + +fd.ModelingSpace("3D") + +material = fd.constitutivelaw.ElasticIsotrop(E, NU) +material.set_density(DENSITY) + +shell_section = fd.constitutivelaw.ShellHomogeneous(material, THICKNESS) +shell_wf = fd.weakform.PlateEquilibrium(shell_section, nlgeom=True) +# # Add numerical damping +# shell_wf += fd.weakform.ArtificialDamping( +# c_stab=0.001, # target about 1% dissipated energy +# energy_fraction=True, +# variables=["Disp", "Rot"], +# ) +sphere_assembly = fd.Assembly.create(shell_wf, mesh[0], name="sphere_shell") + +solid_wf = fd.weakform.StressEquilibrium(material, nlgeom=True) +plate_assembly = fd.Assembly.create(solid_wf, mesh[1], name="plate_solid") + +structural_assembly = sphere_assembly + plate_assembly + + +############################################################################### +# Gravity and IPC contact +# + +gravity_force = [0.0, 0.0, -DENSITY * THICKNESS * GRAVITY] +gravity = fd.constraint.DistributedForce( + mesh[0], + gravity_force, + nlgeom=True, + time_func=lambda t: 1.0, + name="sphere_gravity", +) + +# IPC contact in 3D expects a tri3 surface mesh with node ids matching the full +# mechanical mesh. The shell sphere already is a surface mesh; the plate uses +# only its top face, converted from quad4 to tri3. +plate_top_nodes = mesh[1].find_nodes("Z", PLATE_H) +plate_surface = fd.mesh.extract_surface( + mesh[1], + node_set=plate_top_nodes, + quad2tri=True, +) +ipc_surface = fd.Mesh( + mesh.nodes, + np.vstack((mesh[0].elements, plate_surface.elements)), + "tri3", + ndim=3, +) + +ipc_contact = fd.constraint.IPCContact( + mesh, + surface_mesh=ipc_surface, + dhat=0.5, + dhat_is_relative=False, + friction_coefficient=0.0, + use_ccd=True, +) + +############################################################################### +# Problem, boundary conditions, and output +# + +pb = fd.problem.NonLinearNewmark( + structural_assembly + ipc_contact, + nlgeom=True, +) +# # or with GeneralizedAlpha with high frequency damping: +# pb = fd.problem.NonLinear(structural_assembly + ipc_contact, nlgeom=True) +# pb.set_time_integrator( +# fd.time.SECOND_ORDER, +# fd.time.GeneralizedAlpha(alpha_m=0.0, alpha_f=1.0 / 3.0), +# ) + +plate_nodes = np.unique(mesh[1].elements) +plate_xyz = mesh.nodes[plate_nodes] +plate_support = plate_nodes[ + np.isclose(plate_xyz[:, 2], 0.0) + | np.isclose(plate_xyz[:, 0], 0.0) + | np.isclose(plate_xyz[:, 0], PLATE_X) + | np.isclose(plate_xyz[:, 1], 0.0) + | np.isclose(plate_xyz[:, 1], PLATE_Y) +] + +pb.bc.add(gravity) +pb.bc.add("Dirichlet", plate_support, "Disp", 0) +pb.bc.add("Dirichlet", plate_nodes, "Rot", 0) + +results = pb.add_output( + os.path.join(RESULTS_DIR, "ping_pong_multimesh"), + structural_assembly, + ["Disp", "Stress", "Strain", "Acceleration"], + file_format="fdh5", +) + + +############################################################################### +# Solve and post-process +# + +sphere_nodes = np.unique(mesh[0].elements) +history = {"time": [], "sphere_z": [], "contact_force_norm": []} + + +def track(pb): + disp = pb.get_disp() + history["time"].append(pb.time) + history["sphere_z"].append( + np.mean(mesh.nodes[sphere_nodes, 2] + disp[2, sphere_nodes]) + ) + history["contact_force_norm"].append(np.linalg.norm(ipc_contact.global_vector)) + + +pb.set_nr_criterion("Force", tol=5e-2, max_subiter=20) +pb.nlsolve( + dt=1.0e-3, + tmax=2e-1, + dt_min=2.0e-5, + update_dt=True, + print_info=1, + # interval_output=2.0e-3, # same as dt by default + callback=track, + exec_callback_at_each_iter=False, +) + +# results.load(-1) +# results.plot("Stress", "vm", "Node", scale=1.0) +# fd.viewer(results) diff --git a/fedoo/core/problem.py b/fedoo/core/problem.py index 7e58a4b2..4d6981e2 100644 --- a/fedoo/core/problem.py +++ b/fedoo/core/problem.py @@ -58,6 +58,7 @@ def __init__(self, A=None, B=0, D=0, mesh=None, name="MainProblem", space=None): self._dof_slave = np.array([]) # dof from dirichlet bc or mpc self._dof_free = np.array([]) + self._MFext = None # prepering output demand to export results self._problem_output = _ProblemOutput() diff --git a/fedoo/time/common.py b/fedoo/time/common.py index b14db9c0..947ec163 100644 --- a/fedoo/time/common.py +++ b/fedoo/time/common.py @@ -1,5 +1,7 @@ from dataclasses import dataclass +import numpy as np + @dataclass(frozen=True) class RayleighDamping: @@ -23,6 +25,19 @@ def newmark_acceleration_velocity(beta, gamma, dt, delta_disp, v_n, a_n): stiffness, dissipation) so the recurrence constants live in exactly one place. """ + if delta_disp is None or (np.isscalar(delta_disp) and delta_disp == 0): + delta_disp = np.zeros_like(v_n) + else: + delta_disp = np.asarray(delta_disp) + if delta_disp.shape != np.shape(v_n): + if delta_disp.size != np.size(v_n): + raise ValueError( + "delta_disp and velocity must contain the same number of DOFs." + ) + # Assemblies can map the increment to the global flat DOF layout while + # storing dynamic state variables in their node-shaped representation. + delta_disp = delta_disp.reshape(np.shape(v_n)) + a0 = 1.0 / (beta * dt**2) acc = a0 * (delta_disp - dt * v_n) + (1.0 - 0.5 / beta) * a_n vel = v_n + dt * ((1.0 - gamma) * a_n + gamma * acc) diff --git a/fedoo/weakform/damping_stabilization.py b/fedoo/weakform/damping_stabilization.py index 2cbc1ece..49aa847c 100644 --- a/fedoo/weakform/damping_stabilization.py +++ b/fedoo/weakform/damping_stabilization.py @@ -111,6 +111,8 @@ def initialize(self, assembly, pb): # Start with a very small global fraction self._c_stab = 1e-3 * self.target_ratio self._c_stab_initialized = False + else: + self._c_stab = self.c_stab def set_start(self, assembly, pb): """Update historical variables and adapt c_stab based on energy ratio.""" @@ -118,7 +120,13 @@ def set_start(self, assembly, pb): dt = pb.dtime # 1. Skip if it's the very first initialization or a zero-time step - if dt == 0: + if ( + dt == 0 + or not self._c_stab_initialized + or (np.isscalar(pb.get_A()) and pb.get_A() == 0) + ): + self._c_stab_initialized = True + assembly.sv["_DeltaDisp"] = 0 return # 2. Get the converged displacement increment from the PREVIOUS step @@ -137,7 +145,11 @@ def set_start(self, assembly, pb): # F_damp = c_stab * M* * (delta_u / dt) # M = assembly.get_global_matrix() # delta_E_damp = self._c_stab/dt * (delta_u @ M @ delta_u) - delta_E_damp = delta_u @ assembly.get_global_vector() + damping_force = assembly.get_global_vector() + if np.isscalar(damping_force): + assembly.sv["_DeltaDisp"] = 0 + return + delta_E_damp = delta_u @ damping_force # 5. Adaptive Adjustment # We only adjust if there is significant external work being done