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
Also:
- Code cleaning in dual mesh
- Remove tmp folder of dual mesh
- Fix for Windows (thanks to Nabil)

src/SMESH_I/SMESH_Gen_i.cxx

@@ -2859,14 +2859,15 @@ SMESH::SMESH_Mesh_ptr SMESH_Gen_i::CreateDualMesh(SMESH::SMESH_IDSource_ptr mesh
     ats = "False";
 
   std::string cmd="import salome.smesh.smesh_tools as smt\n";
-  cmd +="smt.smesh_create_dual_mesh(\"" + mesh_ior + "\", \"" +
+  cmd +="smt.smesh_create_dual_mesh(\"" + mesh_ior + "\", r\"" +
         dual_mesh_file.string() + "\", mesh_name=\"" + mesh_name + "\", adapt_to_shape=" + ats + ")";
   MESSAGE(cmd);
 
   PyObject *py_main = PyImport_AddModule("__main__");
   PyObject *py_dict = PyModule_GetDict(py_main);
+  PyObject *local_dict = PyDict_New();
 
-  PyRun_String(cmd.c_str(), Py_file_input, py_dict, py_dict);
+  PyRun_String(cmd.c_str(), Py_file_input, py_dict, local_dict);
 
   if (PyErr_Occurred()) {
     // Restrieving python error
@@ -2931,6 +2932,10 @@ SMESH::SMESH_Mesh_ptr SMESH_Gen_i::CreateDualMesh(SMESH::SMESH_IDSource_ptr mesh
     MESSAGE("Dump of groups");
     SMESH::ListOfGroups_var groups = newMesh->GetGroups();
 
+#ifndef _DEBUG_
+  fs::remove_all(tmp_folder);
+#endif
+
   return newMesh._retn();
 }
 

src/SMESH_SWIG/smeshBuilder.py

@@ -804,7 +804,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)
 

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,60 +15,48 @@ smesh = smeshBuilder.New()
 
 from salome.kernel.logger import Logger
 logger = Logger("salome.smesh.smesh_tools")
-logger.setLevel("DEBUG")
+logger.setLevel("WARNING")
 
-def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name="MESH"):
+# prefix for groups with internal usage
-    """ Create a dual of the mesh in input_file into output_file
+# 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:
-        mesh_ior (string): corba Id of the Tetrahedron mesh
+        mc_mesh_file (MEDFileMesh): mesh on which to read the group
-        output_file (string): dual mesh file
+        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
     """
-    # Import mesh from file
+    try:
-    mesh = salome.orb.string_to_object(mesh_ior)
+        grp_tria = mc_mesh_file.getGroup(grp_level, grp_name)
-    if not mesh:
+    except:
-        raise Exception("Could not find mesh using id: ", mesh_ior)
+        logger.debug("""No group found for %s at level %i.
-
+                     It is normal behaviour for degenerated geom edge."""\
-    shape = mesh.GetShapeToMesh()
+                     %(grp_name, grp_level))
-
+        return None, None
-    # Creating output file
-    logger.debug("Creating file with mesh: "+mesh_name)
-    myfile = mc.MEDFileUMesh()
-    myfile.setName(mesh_name)
-
-
-    # We got a meshProxy so we need to convert pointer to MEDCoupling
-    int_ptr = mesh.ExportMEDCoupling(True, True)
-    dab = mc.FromPyIntPtrToDataArrayByte(int_ptr)
-    mc_mesh_file = mc.MEDFileMesh.New(dab)
-    tetras = mc_mesh_file[0]
-    # End of SMESH -> MEDCoupling part for dualmesh
-
-    tetras = mc.MEDCoupling1SGTUMesh(tetras)
-    polyh = tetras.computeDualMesh()
-
-    ## Adding skin + transfering groups on faces from tetras mesh
-    mesh2d = polyh.buildUnstructured().computeSkin()
-    mesh2d.setName(mesh_name)
-    myfile.setMeshAtLevel(-1, mesh2d)
-
-
-    for grp_name in mc_mesh_file.getGroupsOnSpecifiedLev(-1):
-        # This group is created by the export
-        if grp_name == "Group_Of_All_Faces":
-            logger.debug("Skipping group: "+ grp_name)
-            continue
-        logger.debug("Transferring group: "+ grp_name)
-
-        grp_tria = mc_mesh_file.getGroup(-1, grp_name)
     # Retrieve the nodes in group
-        grp_nodes = grp_tria.computeFetchedNodeIds()
+    orig_grp_nodes = grp_tria.computeFetchedNodeIds()
     # Find all the cells lying on one of the nodes
-        id_grp_poly = mesh2d.getCellIdsLyingOnNodes(grp_nodes, False)
+    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)
+    # 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
@@ -84,44 +67,188 @@ def smesh_create_dual_mesh(mesh_ior, output_file, adapt_to_shape=True, mesh_name
     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())
 
-        myfile.addGroup(-1, id_grp_poly)
+    # 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)
 
-    # Getting list of new points added on the skin
+    return id_grp_poly, nodes_added
-    skin = tetras.buildUnstructured().computeSkin()
+
-    skin_polyh = polyh.buildUnstructured().computeSkin()
+def __projectNodeOnSubshape(nodes, subshape, coords):
-    allNodesOnSkinPolyh = skin_polyh.computeFetchedNodeIds()
+    """ Project nodes on a sub-shape (face or edge) and update the mesh
-    allNodesOnSkin = skin.computeFetchedNodeIds()
+        coordinates
-    ptsAdded = allNodesOnSkinPolyh.buildSubstraction(allNodesOnSkin)
+
-    ptsAddedMesh = mc.MEDCouplingUMesh.Build0DMeshFromCoords( skin_polyh.getCoords()[ptsAdded] )
+    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
+
+    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)
+    if not mesh:
+        raise Exception("Could not find mesh using id: ", mesh_ior)
+
+    shape = mesh.GetShapeToMesh()
 
     if adapt_to_shape:
-        logger.debug("Adapting to shape")
+        faces = geompy.SubShapeAll(shape, geompy.ShapeType["FACE"])
-        ptsAddedCoo = ptsAddedMesh.getCoords()
+        faces_ids = geompy.GetSubShapesIDs(shape, faces)
-        ptsAddedCooModified = ptsAddedCoo[:]
 
-        # Matching faces with their ids
+        # Create group with each face
-        faces = geompy.ExtractShapes(shape, geompy.ShapeType["FACE"], True)
+        # so that we don't need GetFaceNearPoint to get the face to project the
+        # point to
         id2face = {}
-        for face in faces:
+        id2face_edges_ids = {}
-            id2face[face.GetSubShapeIndices()[0]] = face
+        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
 
-        ## Projecting each points added by the dual mesh on the surface it is
+        edges = geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"])
-        # associated with
+        edges_ids = geompy.GetSubShapesIDs(shape, edges)
-        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
+        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)
+    dab = mc.FromPyIntPtrToDataArrayByte(int_ptr)
+    mc_mesh_file = mc.MEDFileMesh.New(dab)
+    tetras = mc_mesh_file[0]
+    # 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)
+
+    polyh_coords = polyh.getCoords()
+
+    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":
+            logger.debug("Skipping group: "+ grp_name)
+            continue
+        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)
+
+        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]
+
+            # for each face, get the edges bounding it
+            grp_poly = mesh2d[id_grp_poly]
+            mesh1d = grp_poly.computeSkin()
+
+            face_edges_id = id2face_edges_ids[face_id]
+            for edge_id in face_edges_id:
+                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)
+
+                # project new nodes on its geom edge
+                # (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)
+
+                    # remove these nodes from the nodes to project on face
+                    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)
+        else:
+            # add the group to write it
+            mc_groups.append(id_grp_poly)
+
+    # 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 grp in mesh_groups:
+            __deleteObj(grp)

test/SMESH_create_dual_mesh_adapt.py

@@ -85,6 +85,8 @@ dual_Mesh_1 = smesh.CreateDualMesh(Mesh_1, 'dual_Mesh_1', True)
 #Comparing volumes
 dual_volume = dual_Mesh_1.GetVolume()
 dual_raw_volume = dual_Mesh_raw_1.GetVolume()
+tetra_volume = Mesh_1.GetVolume()
+print("tetra_volume: ", tetra_volume)
 print("dual_volume: ", dual_volume)
 print("dual_raw_volume: ", dual_raw_volume)
 

test/SMESH_create_dual_mesh_tpipe.py

@@ -79,20 +79,24 @@ 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
+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)
+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)
 
@@ -102,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.14
 
 # Check groups
 dual_Mesh_raw_groups = dual_Mesh_raw_1.GetGroups()
@@ -110,25 +114,44 @@ 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("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.14
+
+# 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()