Improve projection of dual mesh nodes:

- first project nodes on edges, then project on faces
- keep the relation from the sub-shapes to the mesh by groups

src/SMESH_SWIG/smesh_tools.py

@@ -22,6 +22,94 @@ from salome.kernel.logger import Logger
 logger = Logger("salome.smesh.smesh_tools")
 logger.setLevel("DEBUG")
 
+# 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)
+
+    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
+        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.
+        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))
+        return None, None
+    # Retrieve the nodes in group
+    orig_grp_nodes = grp_tria.computeFetchedNodeIds()
+    # Find all the cells lying on one of the nodes
+    id_grp_poly = mesh2d.getCellIdsLyingOnNodes(orig_grp_nodes, False)
+
+    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))
+        return None, None
+
+    # 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.iota()
+
+    # cells that are really in the group
+    id_to_keep = id_poly.buildSubstraction(id_poly_border)
+
+    id_grp_poly = id_grp_poly[id_to_keep]
+
+    id_grp_poly.setName(grp_name.strip())
+
+    # get nodes added on this group
+    grp_poly = mesh2d[id_grp_poly]
+    grp_nodes_poly = grp_poly.computeFetchedNodeIds()
+    nodes_added = grp_nodes_poly.buildSubstraction(orig_grp_nodes)
+
+    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
+
+    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.
+    """
+    for i in nodes:
+        x, y, z = coords[i].getValues()
+        vertex = geompy.MakeVertex(x, y, z)
+        try:
+            prj = geompy.MakeProjection(vertex, subshape)
+        except:
+            logger.warning("Projection failed for %.5f %.5f %.5f but we continue with next node"%(x, y, z))
+            continue
+        new_coor = geompy.PointCoordinates(prj)
+        # update its coordinates in the mesh
+        coords[i] = new_coor
+    pass
+
+def __deleteObj(theObj):
+    """ Delete object from a Study
+
+    Args:
+        theObj (GEOM or SMESH object): object to remove from the study
+    """
+    aStudy = salome.myStudy
+    aStudyBuilder = aStudy.NewBuilder()
+    SO = aStudy.FindObjectIOR(aStudy.ConvertObjectToIOR(theObj))
+    if SO is not None:
+        aStudyBuilder.RemoveObjectWithChildren(SO)
+    pass
+
 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
 
@@ -36,11 +124,33 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
 
     shape = mesh.GetShapeToMesh()
 
-    # Creating output file
+    if adapt_to_shape:
-    logger.debug("Creating file with mesh: "+mesh_name)
+        faces = geompy.SubShapeAll(shape, geompy.ShapeType["FACE"])
-    myfile = mc.MEDFileUMesh()
+        faces_ids = geompy.GetSubShapesIDs(shape, faces)
-    myfile.setName(mesh_name)
 
+        # Create group with each face
+        # 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)
+            id2face[face_id] = face
+            mesh_groups.append(gr_mesh)
+            # get the edges bounding this face
+            # so that we can project the nodes on edges before nodes on faces
+            face_edges = geompy.SubShapeAll(face, geompy.ShapeType["EDGE"])
+            face_edges_ids = geompy.GetSubShapesIDs(shape, face_edges)
+            id2face_edges_ids[face_id] = face_edges_ids
+
+        edges = geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"])
+        edges_ids = geompy.GetSubShapesIDs(shape, edges)
+
+        id2edge = {}
+        for edge, edge_id in zip(edges, edges_ids):
+            gr_mesh = mesh.CreateGroupFromGEOM(SMESH.EDGE, '%sedge_%i'%(__prefix, edge_id), edge)
+            id2edge[edge_id] = edge
+            mesh_groups.append(gr_mesh)
 
     # We got a meshProxy so we need to convert pointer to MEDCoupling
     int_ptr = mesh.ExportMEDCoupling(True, True)
@@ -50,14 +160,21 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
     # End of SMESH -> MEDCoupling part for dualmesh
 
     tetras = mc.MEDCoupling1SGTUMesh(tetras)
+
+    # Create the polyhedra from the tetrahedra (main goal of this function)
     polyh = tetras.computeDualMesh()
 
     ## Adding skin + transfering groups on faces from tetras mesh
     mesh2d = polyh.buildUnstructured().computeSkin()
     mesh2d.setName(mesh_name)
-    myfile.setMeshAtLevel(-1, mesh2d)
 
+    polyh_coords = polyh.getCoords()
 
+    edges_group_names = mc_mesh_file.getGroupsOnSpecifiedLev(-2)
+
+    treated_edges = []
+
+    mc_groups = []
     for grp_name in mc_mesh_file.getGroupsOnSpecifiedLev(-1):
         # This group is created by the export
         if grp_name == "Group_Of_All_Faces":
@@ -65,63 +182,59 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
             continue
         logger.debug("Transferring group: "+ grp_name)
 
-        grp_tria = mc_mesh_file.getGroup(-1, grp_name)
+        # get the polygons ids on the dual mesh from the triangles group
-        # Retrieve the nodes in group
+        id_grp_poly, nodes_added_on_tri = __getIdsGrpDualFromOrig(mc_mesh_file, grp_name, mesh2d, -1)
-        grp_nodes = grp_tria.computeFetchedNodeIds()
-        # Find all the cells lying on one of the nodes
-        id_grp_poly = mesh2d.getCellIdsLyingOnNodes(grp_nodes, False)
 
-        grp_poly = mesh2d[id_grp_poly]
+        if id_grp_poly is not None and grp_name[:4] == __prefix:
+            # This group is on a specific geom face
+            face_id = grp_name.split("_")[-1]
+            face_id = int(face_id)
+            face = id2face[face_id]
 
-        # find the extra face cells, on the border of the group (lying on nodes, but outside the group)
+            # for each face, get the edges bounding it
-        id_poly_border = grp_poly.findCellIdsOnBoundary()
+            grp_poly = mesh2d[id_grp_poly]
+            mesh1d = grp_poly.computeSkin()
 
-        # ids of the cells in grp_poly
+            face_edges_id = id2face_edges_ids[face_id]
-        id_poly=mc.DataArrayInt64.New(grp_poly.getNumberOfCells(), 1)
+            for edge_id in face_edges_id:
-        id_poly.iota()
+                grp_seg_name = "%sedge_%i"%(__prefix, edge_id)
 
-        # cells that are really in the group
+                # get the segs on the dual mesh from the segs group
-        id_to_keep = id_poly.buildSubstraction(id_poly_border)
+                id_grp_seg, nodes_added_on_segs = __getIdsGrpDualFromOrig(mc_mesh_file, grp_seg_name, mesh1d, -2)
 
-        id_grp_poly = id_grp_poly[id_to_keep]
+                # project new nodes on its geom edge
-        id_grp_poly.setName(grp_name.strip())
+                # (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)
+                    else:
+                        treated_edges.append(edge_id)
 
-        myfile.addGroup(-1, id_grp_poly)
+                    # remove these nodes from the nodes to project on face
+                    nodes_added_on_tri = nodes_added_on_tri.buildSubstraction(nodes_added_on_segs)
 
-    # Getting list of new points added on the skin
+            # project new nodes on its geom face
-    skin = tetras.buildUnstructured().computeSkin()
+            __projectNodeOnSubshape(nodes_added_on_tri, face, polyh_coords)
-    skin_polyh = polyh.buildUnstructured().computeSkin()
+        else:
-    allNodesOnSkinPolyh = skin_polyh.computeFetchedNodeIds()
+            # add the group to write it
-    allNodesOnSkin = skin.computeFetchedNodeIds()
+            mc_groups.append(id_grp_poly)
-    ptsAdded = allNodesOnSkinPolyh.buildSubstraction(allNodesOnSkin)
-    ptsAddedMesh = mc.MEDCouplingUMesh.Build0DMeshFromCoords( skin_polyh.getCoords()[ptsAdded] )
-
-    if adapt_to_shape:
-        logger.debug("Adapting to shape")
-        ptsAddedCoo = ptsAddedMesh.getCoords()
-        ptsAddedCooModified = ptsAddedCoo[:]
-
-        # Matching faces with their ids
-        faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True)
-        id2face = {}
-        for face in faces:
-            id2face[face.GetSubShapeIndices()[0]] = face
-
-        ## Projecting each points added by the dual mesh on the surface it is
-        # associated with
-        for i, tup in enumerate(ptsAddedCooModified):
-            vertex = geompy.MakeVertex(*tuple(tup))
-            shapes = geompy.GetShapesNearPoint(shape, vertex,
-                                               geompy.ShapeType["FACE"])
-            prj = geompy.MakeProjection(vertex,
-                                        id2face[shapes.GetSubShapeIndices()[0]])
-            new_coor = geompy.PointCoordinates(prj)
-            ptsAddedCooModified[i] = new_coor
-
-        polyh.getCoords()[ptsAdded] = ptsAddedCooModified
 
+    # Creating output file
+    logger.debug("Creating file with mesh: "+mesh_name)
+    myfile = mc.MEDFileUMesh()
+    myfile.setName(mesh_name)
     polyh.setName(mesh_name)
     myfile.setMeshAtLevel(0, polyh)
+    myfile.setMeshAtLevel(-1, mesh2d)
 
-    logger.debug("Writting dual mesh in :"+output_file)
+    for group in mc_groups:
+        myfile.addGroup(-1, group)
+
+    logger.debug("Writing dual mesh in: "+output_file)
     myfile.write(output_file, 2)
+
+    if adapt_to_shape:
+        # delete temporary groups
+        for gr in mesh_groups:
+            __deleteObj(gr)

test/SMESH_create_dual_mesh_tpipe.py

@@ -83,11 +83,13 @@ Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_3D)
 isDone = Mesh_1.Compute()
 
 # Create groups
+d_geom_groups = {}
 d_mesh_groups = {}
 for geom_group in geom_groups:
   gr = Mesh_1.Group(geom_group)
   gr_name = gr.GetName()
   d_mesh_groups[gr_name] = gr
+  d_geom_groups[gr_name] = geom_group
 
 # Check tetra mesh volume
 mesh_volume = smesh.GetVolume(Mesh_1)
@@ -110,25 +112,42 @@ dual_Mesh_raw_group_names = dual_Mesh_raw_1.GetGroupNames()
 
 assert len(dual_Mesh_raw_groups) == 4
 
-for gr_dual, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names):
+## Create dual mesh with projection on faces
+dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True)
+
+# Check dual mesh volume
+dual_volume = dual_Mesh_1.GetVolume()
+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
+
+# Check groups
+dual_Mesh_groups = dual_Mesh_1.GetGroups()
+dual_Mesh_group_names = dual_Mesh_1.GetGroupNames()
+d_mesh_dual_groups = {}
+for gr_name, gr in zip(dual_Mesh_group_names, dual_Mesh_groups):
+    gr_name = gr_name.strip()
+    d_mesh_dual_groups[gr_name] = gr
+
+for gr_dual_raw, gr_name in zip(dual_Mesh_raw_groups, dual_Mesh_raw_group_names):
   gr_name = gr_name.strip()
   gr_tri = d_mesh_groups[gr_name]
+  gr_dual = d_mesh_dual_groups[gr_name]
+  gr_geom = d_geom_groups[gr_name]
+  area_gr_geom = geompy.BasicProperties(gr_geom)[1]
   area_gr_tri = smesh.GetArea(gr_tri)
+  area_gr_dual_raw = smesh.GetArea(gr_dual_raw)
   area_gr_dual = smesh.GetArea(gr_dual)
   print(gr_name)
+  print("Area geom: ", area_gr_geom)
   print("Area tri: ", area_gr_tri)
+  print("Area dual raw:", area_gr_dual_raw)
   print("Area dual:", area_gr_dual)
-  assert abs(area_gr_tri-area_gr_dual)/area_gr_tri < 1e-3
+  assert abs(area_gr_geom-area_gr_dual)/area_gr_geom < 0.015
+  assert abs(area_gr_tri-area_gr_dual_raw)/area_gr_tri < 1e-3
 
-#dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True)
-
-
-#dual_volume = dual_Mesh_1.GetVolume()
-#dual_raw_volume = dual_Mesh_raw_1.GetVolume()
-#print("dual_volume: ", dual_volume)
-#print("dual_raw_volume: ", dual_raw_volume)
-
-#assert (dual_volume >= dual_raw_volume)
 
 if salome.sg.hasDesktop():
   salome.sg.updateObjBrowser()