diff --git a/examples/09-averaging/00-compute_and_average.py b/examples/09-averaging/00-compute_and_average.py index 1c3c4cd9e8..fe0fc29cd5 100644 --- a/examples/09-averaging/00-compute_and_average.py +++ b/examples/09-averaging/00-compute_and_average.py @@ -99,8 +99,6 @@ def compute_von_mises_then_average(analysis): min_max.inputs.field.connect(avg_von_mises) max_val = min_max.outputs.field_max() - mesh.plot(avg_von_mises) - return max_val.data[0] @@ -139,8 +137,6 @@ def average_then_compute_von_mises(analysis): min_max.inputs.field.connect(avg_von_mises) max_val = min_max.outputs.field_max() - mesh.plot(avg_von_mises) - return max_val.data[0] diff --git a/src/ansys/dpf/core/meshed_region.py b/src/ansys/dpf/core/meshed_region.py index 05d6e84eec..c988e14ac7 100644 --- a/src/ansys/dpf/core/meshed_region.py +++ b/src/ansys/dpf/core/meshed_region.py @@ -108,6 +108,7 @@ def __init__(self, num_nodes=None, num_elements=None, mesh=None, server=None): self._full_grid = None self._elements = None self._nodes = None + self.as_linear = None def _get_scoping(self, loc=locations.nodal): """ @@ -457,20 +458,10 @@ def _as_vtk(self, coordinates=None, as_linear=True, include_ids=False): # consider adding this when scoping request is faster if include_ids: - self._nodeids = self.elements.scoping.ids - self._elementids = self.nodes.scoping.ids - grid["node_ids"] = self._elementids - grid["element_ids"] = self._nodeids - - # Quick fix required to hold onto the data as PyVista does not make a copy. - # All of those now return DPFArrays - if coordinates is None: - coordinates_field = self.nodes.coordinates_field - coordinates = self.nodes.coordinates_field.data - else: - coordinates_field = coordinates - coordinates = coordinates.data - setattr(grid, "_dpf_cache", [coordinates, coordinates_field]) + self._nodeids = self.nodes.scoping.ids + self._elementids = self.elements.scoping.ids + grid["node_ids"] = self._nodeids + grid["element_ids"] = self._elementids return grid diff --git a/src/ansys/dpf/core/plotter.py b/src/ansys/dpf/core/plotter.py index cfa043a8c8..10de6861a1 100644 --- a/src/ansys/dpf/core/plotter.py +++ b/src/ansys/dpf/core/plotter.py @@ -125,7 +125,7 @@ def add_plane(self, plane, field=None, **kwargs): plane[f"{field.name}"] = field.data self._plotter.add_mesh(plane_plot, **kwargs) - def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, **kwargs): + def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, as_linear=True, **kwargs): kwargs = self._set_scalar_bar_title(kwargs) @@ -144,9 +144,17 @@ def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, **kwargs): # 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 + 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 else: - grid = meshed_region._as_vtk(meshed_region.deform_by(deform_by, scale_factor)) + grid = meshed_region._as_vtk( + meshed_region.deform_by(deform_by, scale_factor), as_linear=as_linear + ) # show axes show_axes = kwargs.pop("show_axes", None) @@ -228,6 +236,7 @@ def add_field( deform_by=None, scale_factor=1.0, scale_factor_legend=None, + as_linear=True, **kwargs, ): # Get the field name @@ -276,7 +285,9 @@ def add_field( if not deform_by: grid = meshed_region.grid else: - grid = meshed_region._as_vtk(meshed_region.deform_by(deform_by, scale_factor)) + grid = meshed_region._as_vtk( + meshed_region.deform_by(deform_by, scale_factor), as_linear + ) grid.set_active_scalars(None) self._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in) @@ -481,6 +492,7 @@ def add_mesh(self, meshed_region, deform_by=None, scale_factor=1.0, **kwargs): meshed_region=meshed_region, deform_by=deform_by, scale_factor=scale_factor, + as_linear=True, **kwargs, ) @@ -543,6 +555,7 @@ def add_field( label_point_size=label_point_size, deform_by=deform_by, scale_factor=scale_factor, + as_linear=True, **kwargs, ) @@ -859,11 +872,16 @@ def plot_contour( kwargs_in = _sort_supported_kwargs( bound_method=self._internal_plotter._plotter.add_mesh, **kwargs ) + as_linear = True if deform_by: - grid = mesh._as_vtk(mesh.deform_by(deform_by, scale_factor)) + 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: - grid = mesh.grid + if as_linear != mesh.as_linear: + grid = mesh._as_vtk(mesh.nodes.coordinates_field, as_linear=as_linear) + mesh.as_linear = as_linear + else: + grid = mesh.grid grid.clear_data() self._internal_plotter._plotter.add_mesh(grid, scalars=overall_data, **kwargs_in) diff --git a/src/ansys/dpf/core/vtk_helper.py b/src/ansys/dpf/core/vtk_helper.py index 7b7d956311..174de24083 100644 --- a/src/ansys/dpf/core/vtk_helper.py +++ b/src/ansys/dpf/core/vtk_helper.py @@ -176,7 +176,7 @@ def dpf_mesh_to_vtk_py(mesh, nodes, as_linear): Meshed Region to export to pyVista format nodes : dpf.Field - Field containing the nodes of the mesh. + Field containing the node coordinates of the mesh. as_linear : bool Export quadratic surface elements as linear. @@ -189,9 +189,11 @@ def dpf_mesh_to_vtk_py(mesh, nodes, as_linear): etypes = mesh.elements.element_types_field.data connectivity = mesh.elements.connectivities_field if nodes is None: - nodes = mesh.nodes.coordinates_field.data + coordinates_field = mesh.nodes.coordinates_field + node_coordinates = coordinates_field.data else: - nodes = nodes.data + coordinates_field = nodes + node_coordinates = nodes.data elem_size = np.ediff1d(np.append(connectivity._data_pointer, connectivity.shape)) @@ -207,22 +209,21 @@ def dpf_mesh_to_vtk_py(mesh, nodes, as_linear): insert_ind = np.cumsum(elem_size) insert_ind = np.hstack(([0], insert_ind))[:-1] - # Handle semiparabolic elements - nullmask = connectivity.data == -1 - connectivity.data[nullmask] = 0 - if nullmask.any(): - nodes[0] = np.nan + # if not as_linear: + # # Pre-handle semi-parabolic elements + # semi_mask = connectivity.data == -1 + # if semi_mask.any(): + # # Modify -1 connectivity values + # repeated_data_pointers = connectivity._data_pointer.repeat(repeats=elem_size) + # connectivity.data[semi_mask] = connectivity.data[repeated_data_pointers[semi_mask]] - # For each polyhedron, cell = [nCellFaces, nFace0pts, i, j, k, ..., nFace1pts, i, j, k, ...] - # polys_ind = insert_ind[polyhedron_mask] - # cells = np.take(cells, sorted(set(insert_ind)-set(polys_ind))) # partition cells in vtk format cells = np.insert(connectivity.data, insert_ind, elem_size) # Check if polyhedrons are present if element_types.Polyhedron.value in etypes: cells = np.array(cells) - nodes = np.array(nodes) + nodes = np.array(node_coordinates) insert_ind = insert_ind + np.asarray(list(range(len(insert_ind)))) # Replace in cells values for polyhedron format # [NValuesToFollow, @@ -252,45 +253,93 @@ def compute_offset(): """Return the starting point of a cell in the cells array""" return insert_ind + np.arange(insert_ind.size) + cells_insert_ind = compute_offset() # convert kAns to VTK cell type offset = None if as_linear: + # Map the vtk_cell_type to linear versions of the initial elements types vtk_cell_type = VTK_LINEAR_MAPPING[etypes] - # visualization bug within VTK with quadratic surf cells - ansquad8_mask = etypes == 6 - if np.any(ansquad8_mask): # kAnsQuad8 - - # simply copy the edge node indices to the midside points - offset = compute_offset() - cell_pos = offset[ansquad8_mask] - cells[cell_pos + 5] = cells[cell_pos + 1] - cells[cell_pos + 6] = cells[cell_pos + 2] - cells[cell_pos + 7] = cells[cell_pos + 3] - cells[cell_pos + 8] = cells[cell_pos + 4] - - anstri6_mask = etypes == 4 # kAnsTri6 = 4 - if np.any(anstri6_mask): - if offset is None: - offset = compute_offset() - cell_pos = offset[anstri6_mask] - cells[cell_pos + 4] = cells[cell_pos + 1] - cells[cell_pos + 5] = cells[cell_pos + 2] - cells[cell_pos + 6] = cells[cell_pos + 3] + # Create a global mask of connectivity values to take + mask = np.full(cells.shape, True) + + # Get a mask of quad8 elements in etypes + quad8_mask = etypes == 6 + # If any quad8 + if np.any(quad8_mask): # kAnsQuad8 + # Get the starting indices of quad8 elements in cells + insert_ind_quad8 = cells_insert_ind[quad8_mask] + # insert_ind_quad8 += np.arange(insert_ind_quad8.size) + mask[insert_ind_quad8 + 5] = False + mask[insert_ind_quad8 + 6] = False + mask[insert_ind_quad8 + 7] = False + mask[insert_ind_quad8 + 8] = False + cells[insert_ind_quad8] //= 2 + + tri6_mask = etypes == 4 # kAnsTri6 = 4 + if np.any(tri6_mask): + insert_ind_tri6 = cells_insert_ind[tri6_mask] + # insert_ind_tri6 += np.arange(insert_ind_tri6.size) + mask[insert_ind_tri6 + 4] = False + mask[insert_ind_tri6 + 5] = False + mask[insert_ind_tri6 + 6] = False + cells[insert_ind_tri6] //= 2 + cells = cells[mask] else: vtk_cell_type = VTK_MAPPING[etypes] + # Handle semi-parabolic elements + semi_mask = cells == -1 + if semi_mask.any(): + cells_insert_ind = compute_offset() + # Create a global mask of connectivity values to take + mask = np.full(cells.shape, True) + # Build a map of size cells with repeated element beginning index + repeated_insert_ind = cells_insert_ind.repeat(repeats=elem_size + 1) + # Apply the semi-mask to get a unique set of indices of semi-parabolic elements in cells + semi_indices_in_cells = np.array(list(set(repeated_insert_ind[semi_mask]))) + semi_sizes = cells[semi_indices_in_cells] + semi_quad8 = semi_sizes == 8 + if semi_quad8.any(): + mask[semi_indices_in_cells[semi_quad8] + 5] = False + mask[semi_indices_in_cells[semi_quad8] + 6] = False + mask[semi_indices_in_cells[semi_quad8] + 7] = False + mask[semi_indices_in_cells[semi_quad8] + 8] = False + cells[semi_indices_in_cells[semi_quad8]] //= 2 + + quad8_mask = etypes == 6 + semi_quad8_mask = (cells[cells_insert_ind] == 4) & quad8_mask + vtk_cell_type[semi_quad8_mask] = VTK_LINEAR_MAPPING[6] + semi_tri6 = semi_sizes == 6 + if semi_tri6.any(): + mask[semi_indices_in_cells[semi_tri6] + 4] = False + mask[semi_indices_in_cells[semi_tri6] + 5] = False + mask[semi_indices_in_cells[semi_tri6] + 6] = False + cells[semi_indices_in_cells[semi_tri6]] //= 2 + + tri6_mask = etypes == 4 + semi_tri6_mask = (cells[cells_insert_ind] == 3) & tri6_mask + vtk_cell_type[semi_tri6_mask] = VTK_LINEAR_MAPPING[4] + # Update cells with the mask + cells = cells[mask] + # different treatment depending on the version of vtk if VTK9: # compute offset array when < VTK v9 - return pv.UnstructuredGrid(cells, vtk_cell_type, nodes) + grid = pv.UnstructuredGrid(cells, vtk_cell_type, node_coordinates) + + # Quick fix required to hold onto the data as PyVista does not make a copy. + # All of those now return DPFArrays + setattr(grid, "_dpf_cache", [node_coordinates, coordinates_field]) + + return grid # might be computed when checking for VTK quadratic bug if offset is None: offset = compute_offset() - return pv.UnstructuredGrid(offset, cells, vtk_cell_type, nodes) + return pv.UnstructuredGrid(offset, cells, vtk_cell_type, node_coordinates) def dpf_mesh_to_vtk(mesh, nodes=None, as_linear=True): diff --git a/tests/test_meshregion.py b/tests/test_meshregion.py index 98acbeeaf9..4bb6b6e750 100644 --- a/tests/test_meshregion.py +++ b/tests/test_meshregion.py @@ -65,7 +65,7 @@ def test_get_set_unit_meshedregion(simple_bar_model): def test_get_node_meshedregion(simple_bar_model): mesh = simple_bar_model.metadata.meshed_region node = mesh.nodes.node_by_index(1) - scop = mesh._get_scoping(dpf.core.locations.nodal) + scop = mesh.nodes.scoping assert node.id == scop.id(1) assert node.index == 1