diff --git a/examples/06-plotting/00-basic_plotting.py b/examples/06-plotting/00-basic_plotting.py index 83fae5c9a0d..47dbe74ec6f 100644 --- a/examples/06-plotting/00-basic_plotting.py +++ b/examples/06-plotting/00-basic_plotting.py @@ -52,9 +52,8 @@ # - The "off_screen" keyword only works when "notebook=False" to prevent the GUI from appearing. -# Plot a field on its supporting mesh (field location must be Elemental or Nodal) +# Plot a field on its supporting mesh stress = model.results.stress() -stress.inputs.requested_location.connect(dpf.locations.nodal) fc = stress.outputs.fields_container() field = fc[0] field.plot(notebook=False, shell_layers=None, show_axes=True, title="Field", text="Field plot") diff --git a/src/ansys/dpf/core/elements.py b/src/ansys/dpf/core/elements.py index 2cff1b09988..074f4e051a7 100644 --- a/src/ansys/dpf/core/elements.py +++ b/src/ansys/dpf/core/elements.py @@ -677,7 +677,7 @@ def map_scoping(self, external_scope): Examples -------- - Return the indices that map a field to an elements collection. + Return the indices that map a field to an Elements collection. >>> import ansys.dpf.core as dpf >>> from ansys.dpf.core import examples @@ -689,7 +689,7 @@ def map_scoping(self, external_scope): """ if external_scope.location in ["Nodal", "NodalElemental"]: - raise ValueError('Input scope location must be "Nodal"') + raise ValueError('Input scope location must be "Elemental"') arr = np.array(list(map(self.mapping_id_to_index.get, external_scope.ids))) mask = arr != None ind = arr[mask].astype(np.int32) diff --git a/src/ansys/dpf/core/meshed_region.py b/src/ansys/dpf/core/meshed_region.py index 80ba95ad61e..a6b1c54c70c 100644 --- a/src/ansys/dpf/core/meshed_region.py +++ b/src/ansys/dpf/core/meshed_region.py @@ -33,6 +33,8 @@ import traceback import warnings +import numpy as np + from ansys.dpf.core import field, property_field, scoping, server as server_module from ansys.dpf.core.cache import class_handling_cache from ansys.dpf.core.check_version import meets_version, server_meet_version, version_requires @@ -713,3 +715,48 @@ def is_empty(self) -> bool: no_elements = self.elements.n_elements == 0 no_nodes = self.nodes.n_nodes == 0 return no_nodes and no_faces and no_elements + + def location_data_len(self, location: locations) -> int: + """Return the data length for a given mesh location. + + Accepted mesh locations are nodal, elemental, faces, and elemental_nodal. + + Parameters + ---------- + location: + The mesh location to compute data length for. + Can be nodal, elemental, faces, or elemental_nodal. + + Returns + ------- + data_size: + If location is nodal, return the number of nodes. + If location is elemental, return the number of elements. + If location is faces, return the number of faces. + If location is elemental nodal, return the sum of the number of nodes per element. + """ + if location == locations.nodal: + return len(self.nodes) + elif location == locations.elemental: + return len(self.elements) + elif location == locations.faces: + return len(self.faces) + elif location == locations.elemental_nodal: + return np.sum(self.get_elemental_nodal_size_list()) + elif location == locations.overall: + return len(self.elements) + else: + raise TypeError(f"Location {location} is not recognized.") + + def get_elemental_nodal_size_list(self) -> np.ndarray: + """Return the array of number of nodes per element in the mesh.""" + # Get the field of element types + element_types_field = self.elements.element_types_field + # get the number of nodes for each possible element type + size_map = dict( + [(e_type.value, element_types.descriptor(e_type).n_nodes) for e_type in element_types] + ) + keys = list(size_map.keys()) + sort_idx = np.argsort(keys) + idx = np.searchsorted(keys, element_types_field.data, sorter=sort_idx) + return np.asarray(list(size_map.values()))[sort_idx][idx] diff --git a/src/ansys/dpf/core/plotter.py b/src/ansys/dpf/core/plotter.py index 35051f25911..f0c5b331a69 100644 --- a/src/ansys/dpf/core/plotter.py +++ b/src/ansys/dpf/core/plotter.py @@ -147,6 +147,7 @@ def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, as_linear=Tr grid = meshed_region._as_vtk( meshed_region.nodes.coordinates_field, as_linear=as_linear ) + meshed_region._full_grid = grid meshed_region.as_linear = as_linear else: grid = meshed_region.grid @@ -317,8 +318,14 @@ def add_field( show_min = False elif location == locations.overall: mesh_location = meshed_region.elements + elif location == locations.elemental_nodal: + mesh_location = meshed_region.elements + # If ElementalNodal, first extend results to mid-nodes + field = dpf.core.operators.averaging.extend_to_mid_nodes(field=field).eval() else: - raise ValueError("Only elemental, nodal or faces location are supported for plotting.") + raise ValueError( + "Only elemental, elemental nodal, nodal, faces, or overall location are supported for plotting." + ) # Treat multilayered shells if not isinstance(shell_layer, eshell_layers): @@ -333,13 +340,27 @@ def add_field( ) field = change_shell_layer_op.get_output(0, core.types.field) + location_data_len = meshed_region.location_data_len(location) component_count = field.component_count if component_count > 1: - overall_data = np.full((len(mesh_location), component_count), np.nan) + overall_data = np.full((location_data_len, component_count), np.nan) else: - overall_data = np.full(len(mesh_location), np.nan) + overall_data = np.full(location_data_len, np.nan) if location != locations.overall: ind, mask = mesh_location.map_scoping(field.scoping) + + # Rework ind and mask to take into account n_nodes per element if ElementalNodal + if location == locations.elemental_nodal: + n_nodes_list = meshed_region.get_elemental_nodal_size_list().astype(np.int32) + first_index = np.insert(np.cumsum(n_nodes_list)[:-1], 0, 0).astype(np.int32) + mask_2 = np.asarray( + [mask_i for i, mask_i in enumerate(mask) for _ in range(n_nodes_list[ind[i]])] + ) + ind_2 = np.asarray( + [first_index[ind_i] + j for ind_i in ind for j in range(n_nodes_list[ind_i])] + ) + mask = mask_2 + ind = ind_2 overall_data[ind] = field.data[mask] else: overall_data[:] = field.data[0] @@ -348,12 +369,22 @@ def add_field( # Have to remove any active scalar field from the pre-existing grid object, # otherwise we get two scalar bars when calling several plot_contour on the same mesh # but not for the same field. The PyVista UnstructuredGrid keeps memory of it. - if not deform_by: - grid = meshed_region.grid - else: + if location == locations.elemental_nodal: + as_linear = False + if deform_by: grid = meshed_region._as_vtk( - meshed_region.deform_by(deform_by, scale_factor), as_linear + meshed_region.deform_by(deform_by, scale_factor), as_linear=as_linear ) + else: + if as_linear != meshed_region.as_linear: + grid = meshed_region._as_vtk( + meshed_region.nodes.coordinates_field, as_linear=as_linear + ) + meshed_region.as_linear = as_linear + else: + grid = meshed_region.grid + if location == locations.elemental_nodal: + grid = grid.shrink(1.0) grid.set_active_scalars(None) self._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in) @@ -976,15 +1007,14 @@ def plot_contour( import warnings warnings.simplefilter("ignore") - if isinstance(field_or_fields_container, (dpf.core.Field, dpf.core.FieldsContainer)): fields_container = None if isinstance(field_or_fields_container, dpf.core.Field): fields_container = dpf.core.FieldsContainer( server=field_or_fields_container._server ) - fields_container.add_label(DefinitionLabels.time) - fields_container.add_field({DefinitionLabels.time: 1}, field_or_fields_container) + fields_container.add_label("id") + fields_container.add_field({"id": 1}, field_or_fields_container) elif isinstance(field_or_fields_container, dpf.core.FieldsContainer): fields_container = field_or_fields_container else: @@ -1022,43 +1052,96 @@ def plot_contour( unit = field.unit break + # If ElementalNodal, first extend results to mid-nodes + if location == locations.elemental_nodal: + fields_container = dpf.core.operators.averaging.extend_to_mid_nodes_fc( + fields_container=fields_container + ).eval() + + location_data_len = mesh.location_data_len(location) if location == locations.nodal: mesh_location = mesh.nodes elif location == locations.elemental: mesh_location = mesh.elements elif location == locations.faces: mesh_location = mesh.faces + elif location == locations.elemental_nodal: + mesh_location = mesh.elements else: - raise ValueError("Only elemental, nodal or faces location are supported for plotting.") + raise ValueError( + "Only elemental, elemental nodal, nodal or faces location are supported for plotting." + ) # pre-loop: check if shell layers for each field, if yes, set the shell layers - changeOp = core.Operator("change_shellLayers") - for field in fields_container: - shell_layer_check = field.shell_layers - if shell_layer_check in [ - eshell_layers.topbottom, - eshell_layers.topbottommid, - ]: - changeOp.inputs.fields_container.connect(fields_container) - sl = eshell_layers.top - if shell_layers is not None: - if not isinstance(shell_layers, eshell_layers): - raise TypeError( - "shell_layer attribute must be a core.shell_layers instance." - ) - sl = shell_layers - changeOp.inputs.e_shell_layer.connect(sl.value) # top layers taken - fields_container = changeOp.get_output(0, core.types.fields_container) - break + changeOp = core.operators.utility.change_shell_layers() + if location == locations.elemental_nodal: + # change_shell_layers does not support elemental_nodal when given a fields_container + new_fields_container = dpf.core.FieldsContainer() + for l in fields_container.labels: + new_fields_container.add_label(l) + for i, field in enumerate(fields_container): + label_space_i = fields_container.get_label_space(i) + shell_layer_check = field.shell_layers + if shell_layer_check in [ + eshell_layers.topbottom, + eshell_layers.topbottommid, + ]: + changeOp.inputs.fields_container.connect(field) + changeOp.inputs.merge.connect(True) + sl = eshell_layers.top + if shell_layers is not None: + if not isinstance(shell_layers, eshell_layers): + raise TypeError( + "shell_layer attribute must be a core.shell_layers instance." + ) + sl = shell_layers + changeOp.inputs.e_shell_layer.connect(sl.value) # top layers taken + field = changeOp.get_output(0, core.types.field) + new_fields_container.add_field(label_space=label_space_i, field=field) + fields_container = new_fields_container + else: + for field in fields_container: + shell_layer_check = field.shell_layers + if shell_layer_check in [ + eshell_layers.topbottom, + eshell_layers.topbottommid, + ]: + changeOp.inputs.fields_container.connect(fields_container) + sl = eshell_layers.top + if shell_layers is not None: + if not isinstance(shell_layers, eshell_layers): + raise TypeError( + "shell_layer attribute must be a core.shell_layers instance." + ) + sl = shell_layers + changeOp.inputs.e_shell_layer.connect(sl.value) # top layers taken + fields_container = changeOp.get_output(0, core.types.fields_container) + break # Merge field data into a single array if component_count > 1: - overall_data = np.full((len(mesh_location), component_count), np.nan) + overall_data = np.full((location_data_len, component_count), np.nan) else: - overall_data = np.full(len(mesh_location), np.nan) + overall_data = np.full(location_data_len, np.nan) + + # field._data_pointer gives the first index of each entity data + # (should be of size nb_elements) for field in fields_container: ind, mask = mesh_location.map_scoping(field.scoping) + if location == locations.elemental_nodal: + # Rework ind and mask to take into account n_nodes per element + # entity_index_map = field._data_pointer + n_nodes_list = mesh.get_elemental_nodal_size_list().astype(np.int32) + first_index = np.insert(np.cumsum(n_nodes_list)[:-1], 0, 0).astype(np.int32) + mask_2 = np.asarray( + [mask_i for i, mask_i in enumerate(mask) for _ in range(n_nodes_list[ind[i]])] + ) + ind_2 = np.asarray( + [first_index[ind_i] + j for ind_i in ind for j in range(n_nodes_list[ind_i])] + ) + mask = mask_2 + ind = ind_2 overall_data[ind] = field.data[mask] # create the plotter and add the meshes @@ -1087,15 +1170,20 @@ def plot_contour( bound_method=self._internal_plotter._plotter.add_mesh, **kwargs ) as_linear = True + if location == locations.elemental_nodal: + as_linear = False if deform_by: grid = mesh._as_vtk(mesh.deform_by(deform_by, scale_factor), as_linear=as_linear) self._internal_plotter.add_scale_factor_legend(scale_factor, **kwargs) else: if as_linear != mesh.as_linear: grid = mesh._as_vtk(mesh.nodes.coordinates_field, as_linear=as_linear) + mesh._full_grid = grid mesh.as_linear = as_linear else: grid = mesh.grid + if location == locations.elemental_nodal: + grid = grid.shrink(1.0) grid.clear_data() self._internal_plotter._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in) diff --git a/tests/test_plotter.py b/tests/test_plotter.py index 70599d73945..5765989645a 100644 --- a/tests/test_plotter.py +++ b/tests/test_plotter.py @@ -31,6 +31,7 @@ from conftest import ( SERVERS_VERSION_GREATER_THAN_OR_EQUAL_TO_5_0, SERVERS_VERSION_GREATER_THAN_OR_EQUAL_TO_7_0, + SERVERS_VERSION_GREATER_THAN_OR_EQUAL_TO_10_0, running_docker, ) @@ -129,7 +130,7 @@ def test_plotter_on_fields_container_elemental(allkindofcomplexity): avg_op.inputs.fields_container.connect(stress.outputs.fields_container) fc = avg_op.outputs.fields_container() pl = Plotter(model.metadata.meshed_region) - cpos = pl.plot_contour(fc) + _ = pl.plot_contour(fc) @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") @@ -141,7 +142,7 @@ def test_plotter_on_fields_container_nodal(allkindofcomplexity): avg_op.inputs.fields_container.connect(stress.outputs.fields_container) fc = avg_op.outputs.fields_container() pl = Plotter(model.metadata.meshed_region) - cpos = pl.plot_contour(fc) + _ = pl.plot_contour(fc) @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") @@ -199,6 +200,70 @@ def test_field_nodal_plot(allkindofcomplexity): remove_picture(picture) +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_field_elemental_nodal_plot_simple(simple_bar): + model = Model(simple_bar) + stress = model.results.element_nodal_forces() + fc = stress.outputs.fields_container() + f = fc[0] + f.plot() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_field_elemental_nodal_plot_scoped(simple_bar): + model = Model(simple_bar) + mesh_scoping = dpf.core.mesh_scoping_factory.elemental_scoping( + element_ids=list(range(1501, 3001)) + ) + stress = model.results.element_nodal_forces.on_mesh_scoping(mesh_scoping) + fc = stress.eval() + f = fc[0] + f.plot() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_field_elemental_nodal_plot_multiple_solid_types(): + from ansys.dpf.core import examples + + model = dpf.core.Model(examples.download_hemisphere()) + stress = model.results.stress() + fc = stress.outputs.fields_container() + f = fc[0] + f.plot() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_field_elemental_nodal_plot_shells(): + from ansys.dpf.core import examples + + model = dpf.core.Model(examples.download_pontoon()) + stress = model.results.stress() + fc = stress.outputs.fields_container() + f = fc[0] + f.plot() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +@pytest.mark.skipif(not SERVERS_VERSION_GREATER_THAN_OR_EQUAL_TO_10_0, reason="Old bug before 25R2") +def test_field_elemental_nodal_plot_multi_shells(multishells): + fc = core.operators.result.stress(data_sources=core.DataSources(multishells)).eval() + from ansys.dpf.core.plotter import Plotter + + field = fc[0] + plt = Plotter(field.meshed_region) + plt.plot_contour(fc) + field.plot() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +@pytest.mark.skipif(not SERVERS_VERSION_GREATER_THAN_OR_EQUAL_TO_10_0, reason="Old bug before 25R2") +def test_dpf_plotter_add_field_elemental_nodal_multi_shells(multishells): + fc: core.FieldsContainer = core.operators.result.stress( + data_sources=core.DataSources(multishells), + ).eval() + fc.plot() + + @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") def test_field_solid_plot(allkindofcomplexity): model = Model(allkindofcomplexity) @@ -304,6 +369,54 @@ def test_dpf_plotter_add_field_change_shell_layer(multishells): plt.add_field(field=field) +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_dpf_plotter_add_field_elemental_nodal_plot_simple(simple_bar): + field: core.Field = core.operators.result.element_nodal_forces( + data_sources=core.DataSources(simple_bar), + ).eval()[0] + plt = DpfPlotter() + plt.add_field(field=field) + plt.show_figure() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_dpf_plotter_add_field_elemental_nodal_plot_scoped(simple_bar): + mesh_scoping = dpf.core.mesh_scoping_factory.elemental_scoping( + element_ids=list(range(1501, 3001)) + ) + field: core.Field = core.operators.result.element_nodal_forces( + data_sources=core.DataSources(simple_bar), + mesh_scoping=mesh_scoping, + ).eval()[0] + plt = DpfPlotter() + plt.add_field(field=field) + plt.show_figure() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_dpf_plotter_add_field_elemental_nodal_multiple_solids(): + from ansys.dpf.core import examples + + field: core.Field = core.operators.result.stress( + data_sources=core.DataSources(examples.download_hemisphere()), + ).eval()[0] + plt = DpfPlotter() + plt.add_field(field=field) + plt.show_figure() + + +@pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") +def test_dpf_plotter_add_field_elemental_nodal_shells(): + from ansys.dpf.core import examples + + field: core.Field = core.operators.result.stress( + data_sources=core.DataSources(examples.download_pontoon()), + ).eval()[0] + plt = DpfPlotter() + plt.add_field(field=field) + plt.show_figure() + + @pytest.mark.skipif(not HAS_PYVISTA, reason="Please install pyvista") def test_plot_fieldscontainer_on_mesh_scoping(multishells): model = core.Model(multishells)