From 11d9d2377ca634335272c0d69eda9809171e8cf3 Mon Sep 17 00:00:00 2001 From: Christophe Bourcier Date: Wed, 16 Nov 2022 10:44:02 +0100 Subject: [PATCH] Code cleaning in dual mesh --- src/SMESH_SWIG/smeshBuilder.py | 1 - src/SMESH_SWIG/smesh_tools.py | 70 +++++++++++++++++----------- test/SMESH_create_dual_mesh_tpipe.py | 14 ++++-- 3 files changed, 51 insertions(+), 34 deletions(-) diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index 62db88951..99615c745 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -789,7 +789,6 @@ class smeshBuilder( SMESH._objref_SMESH_Gen, object ): """ if isinstance( mesh, Mesh ): mesh = mesh.GetMesh() - print("calling createdualmesh from Python") dualMesh = SMESH._objref_SMESH_Gen.CreateDualMesh(self, mesh, meshName, adaptToShape) return Mesh(self, self.geompyD, dualMesh) diff --git a/src/SMESH_SWIG/smesh_tools.py b/src/SMESH_SWIG/smesh_tools.py index 563d10ba8..19d9fb6fd 100644 --- a/src/SMESH_SWIG/smesh_tools.py +++ b/src/SMESH_SWIG/smesh_tools.py @@ -1,12 +1,7 @@ #!/usr/bin/env python3 -import sys import salome import medcoupling as mc -from math import pi -import numpy as np - -#salome.salome_init() import GEOM from salome.geom import geomBuilder @@ -20,30 +15,34 @@ smesh = smeshBuilder.New() from salome.kernel.logger import Logger logger = Logger("salome.smesh.smesh_tools") -logger.setLevel("DEBUG") +logger.setLevel("WARNING") # prefix for groups with internal usage # i.e. used to transfer the faces and edges sub-shapes ids to the mesh __prefix = "____" def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level): - """ Identify the polygonal cells ids matching the original group on the original mesh (before dual mesh) + """ Identify the polygonal cells ids matching the original group on the + original mesh (before dual mesh) Args: mc_mesh_file (MEDFileMesh): mesh on which to read the group grp_name (string): name of the group to read - mesh2d (MEDCouplingUMesh): mesh at lower level (-1 or -2) containing faces or segments cells + mesh2d (MEDCouplingUMesh): mesh at lower level (-1 or -2) containing + faces or segments cells grp_level (int): level on which to load the group (-1 or -2) Returns: - id_grp_poly (DataArrayInt64): ids of cells mathcing the group. None if the group has not been found. + id_grp_poly (DataArrayInt64): ids of cells mathcing the group. None if + the group has not been found. nodes_added (DataArrayInt64): new nodes added on the dual mesh """ try: grp_tria = mc_mesh_file.getGroup(grp_level, grp_name) except: logger.debug("""No group found for %s at level %i. - It is normal behaviour for degenerated geom edge."""%(grp_name, grp_level)) + It is normal behaviour for degenerated geom edge."""\ + %(grp_name, grp_level)) return None, None # Retrieve the nodes in group orig_grp_nodes = grp_tria.computeFetchedNodeIds() @@ -52,14 +51,16 @@ def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level): grp_poly = mesh2d[id_grp_poly] if grp_poly.getNumberOfCells() == 0: - logger.debug("""No poly cell found for %s at level %i."""%(grp_name, grp_level)) + logger.debug("""No poly cell found for %s at level %i."""\ + %(grp_name, grp_level)) return None, None - # find the extra face cells, on the border of the group (lying on nodes, but outside the group) + # find the extra face cells, on the border of the group (lying on nodes, + # but outside the group) id_poly_border = grp_poly.findCellIdsOnBoundary() # ids of the cells in grp_poly - id_poly=mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) + id_poly = mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1) id_poly.iota() # cells that are really in the group @@ -77,12 +78,14 @@ def __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, grp_level): return id_grp_poly, nodes_added def __projectNodeOnSubshape(nodes, subshape, coords): - """ Project nodes on a sub-shape (face or edge) and update the mesh coordinates + """ Project nodes on a sub-shape (face or edge) and update the mesh + coordinates Args: nodes (DataArrayInt): nodes ids to project subshape (GEOM object): face or edge on which to project the nodes - coords (DataArrayDouble): coordinates of the mesh to update. These coordinates are modified inside this function. + coords (DataArrayDouble): coordinates of the mesh to update. These + coordinates are modified inside this function. """ for i in nodes: x, y, z = coords[i].getValues() @@ -110,12 +113,15 @@ def __deleteObj(theObj): aStudyBuilder.RemoveObjectWithChildren(SO) pass -def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"): +def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, + mesh_name="MESH"): """ Create a dual of the mesh in input_file into output_file Args: mesh_ior (string): corba Id of the Tetrahedron mesh output_file (string): dual mesh file + adapt_to_shape (bool): If True will project boundary points on shape + mesh_name (string): Name of the dual Mesh """ # Import mesh from file mesh = salome.orb.string_to_object(mesh_ior) @@ -129,12 +135,15 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name faces_ids = geompy.GetSubShapesIDs(shape, faces) # Create group with each face - # so that we don't need GetFaceNearPoint to get the face to project the point to + # so that we don't need GetFaceNearPoint to get the face to project the + # point to id2face = {} id2face_edges_ids = {} mesh_groups = [] for face, face_id in zip(faces, faces_ids): - gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE, '%sface_%i'%(__prefix, face_id), face) + gr_mesh = mesh.CreateGroupFromGEOM(SMESH.FACE, + '%sface_%i'%(__prefix, face_id), + face) id2face[face_id] = face mesh_groups.append(gr_mesh) # get the edges bounding this face @@ -148,7 +157,9 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name id2edge = {} for edge, edge_id in zip(edges, edges_ids): - gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, '%sedge_%i'%(__prefix, edge_id), edge) + gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, + '%sedge_%i'%(__prefix, edge_id), + edge) id2edge[edge_id] = edge mesh_groups.append(gr_mesh) @@ -170,8 +181,6 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name polyh_coords = polyh.getCoords() - edges_group_names = mc_mesh_file.getGroupsOnSpecifiedLev(-2) - treated_edges = [] mc_groups = [] @@ -183,7 +192,8 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name logger.debug("Transferring group: "+ grp_name) # get the polygons ids on the dual mesh from the triangles group - id_grp_poly, nodes_added_on_tri = __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1) + id_grp_poly, nodes_added_on_tri = \ + __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1) if id_grp_poly is not None and grp_name[:4] == __prefix: # This group is on a specific geom face @@ -200,19 +210,23 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name grp_seg_name = "%sedge_%i"%(__prefix, edge_id) # get the segs on the dual mesh from the segs group - id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(mc_mesh_file, grp_seg_name, mesh1d, -2) + id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(\ + mc_mesh_file, grp_seg_name, mesh1d, -2) # project new nodes on its geom edge - # (if the group exists on this edge and it has not already been treated) + # (if the group exists on this edge and it has not already been + # treated) if id_grp_seg is not None: if edge_id not in treated_edges: edge = id2edge[edge_id] - __projectNodeOnSubshape(nodes_added_on_segs, edge, polyh_coords) + __projectNodeOnSubshape(nodes_added_on_segs, + edge, polyh_coords) else: treated_edges.append(edge_id) # remove these nodes from the nodes to project on face - nodes_added_on_tri = nodes_added_on_tri.buildSubstraction(nodes_added_on_segs) + nodes_added_on_tri = \ + nodes_added_on_tri.buildSubstraction(nodes_added_on_segs) # project new nodes on its geom face __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords) @@ -236,5 +250,5 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name if adapt_to_shape: # delete temporary groups - for gr in mesh_groups: - __deleteObj(gr) + for grp in mesh_groups: + __deleteObj(grp) diff --git a/test/SMESH_create_dual_mesh_tpipe.py b/test/SMESH_create_dual_mesh_tpipe.py index 3d121876d..3e57e86b1 100644 --- a/test/SMESH_create_dual_mesh_tpipe.py +++ b/test/SMESH_create_dual_mesh_tpipe.py @@ -79,7 +79,9 @@ Mesh_1 = smesh.Mesh(tpipe, "Mesh_coarse") Mesh_1.Triangle(algo=smeshBuilder.NETGEN_1D2D) -Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) +algo3d = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D) +algo3d.MaxElementVolume(0.0002) + isDone = Mesh_1.Compute() # Create groups @@ -92,9 +94,9 @@ for geom_group in geom_groups: d_geom_groups[gr_name] = geom_group # Check tetra mesh volume -mesh_volume = smesh.GetVolume(Mesh_1) +tetra_volume = smesh.GetVolume(Mesh_1) shape_volume = geompy.BasicProperties(tpipe)[2] -assert abs(mesh_volume-shape_volume)/shape_volume < 0.05 +assert abs(tetra_volume-shape_volume)/shape_volume < 0.05 dual_Mesh_raw_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_raw_1', False) @@ -104,7 +106,7 @@ assert dual_Mesh_raw_1.NbPolyhedrons() == dual_Mesh_raw_1.NbVolumes() # Check dual mesh volume dual_raw_volume = dual_Mesh_raw_1.GetVolume() -assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.11 +assert abs(dual_raw_volume-shape_volume)/shape_volume < 0.12 # Check groups dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups() @@ -117,11 +119,13 @@ dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True) # Check dual mesh volume dual_volume = dual_Mesh_1.GetVolume() +print("shape_volume: ", shape_volume) +print("tetra_volume: ", tetra_volume) print("dual_volume: ", dual_volume) print("dual_raw_volume: ", dual_raw_volume) assert (dual_volume >= dual_raw_volume) -assert abs(dual_volume-shape_volume)/shape_volume < 0.11 +assert abs(dual_volume-shape_volume)/shape_volume < 0.12 # Check groups dual_Mesh_groups = dual_Mesh_1.GetGroups() -- 2.39.2