Adding code for Python parallel mesh with test + removing nbMesherThreads

idl/SMESH_Mesh.idl

@@ -905,11 +905,6 @@ module SMESH
     void SetNbThreads(in long nbThreads);
     /*!
 
-     * \brief Set Number of Threads for mesher
-     */
-    void SetMesherNbThreads(in long nbThreads);
-
-
     /*!
      * Get mesh description
      */

src/SMESH/SMESH_Gen.cxx

@@ -279,7 +279,7 @@ bool SMESH_Gen::parallelComputeSubMeshes(
 
   TopAbs_ShapeEnum previousShapeType = TopAbs_VERTEX;
   int nbThreads = aMesh.GetNbThreads();
-  MESSAGE("Compute submeshes with threads: " << nbThreads << " mesher: " << aMesh.GetMesherNbThreads());
+  MESSAGE("Compute submeshes with threads: " << nbThreads);
 
 
   smIt = shapeSM->getDependsOnIterator(includeSelf, !complexShapeFirst);

src/SMESH/SMESH_Mesh.cxx

@@ -87,7 +87,7 @@ namespace fs=boost::filesystem;
 #define MAX_MED_GROUP_NAME_LENGTH 80
 
 #ifdef _DEBUG_
-static int MYDEBUG = 1;
+static int MYDEBUG = 0;
 #else
 static int MYDEBUG = 0;
 #endif

src/SMESH/SMESH_Mesh.hxx

@@ -392,9 +392,6 @@ class SMESH_EXPORT SMESH_Mesh
   int GetNbThreads(){return _NbThreads;};
   void SetNbThreads(int nbThreads){_NbThreads=nbThreads;};
 
-  int GetMesherNbThreads(){return _MesherNbThreads;};
-  void SetMesherNbThreads(int nbThreads){_MesherNbThreads=nbThreads;};
-
   void InitPoolThreads(){_pool = new boost::asio::thread_pool(_NbThreads);};
   void DeletePoolThreads(){delete _pool;};
 
@@ -456,7 +453,6 @@ protected:
   // Mutex for multhitreading write in SMESH_Mesh
   std::mutex _my_lock;
   int _NbThreads=0;
-  int _MesherNbThreads=0;
 
 protected:
   SMESH_Mesh();

src/SMESH_I/SMESH_2smeshpy.cxx

@@ -2213,7 +2213,7 @@ bool _pyMesh::NeedMeshAccess( const Handle(_pyCommand)& theCommand )
         "GetElemFaceNodes", "GetFaceNormal", "FindElementByNodes",
         "IsPoly","IsQuadratic","BaryCenter","GetHypothesisList", "SetAutoColor", "GetAutoColor",
         "Clear", "ConvertToStandalone", "GetMeshOrder", "SetMeshOrder",
-        "SetNbThreads", "SetMesherNbThreads"
+        "SetNbThreads"
         ,"" }; // <- mark of end
     sameMethods.Insert( names );
   }

src/SMESH_I/SMESH_Mesh_i.cxx

@@ -7046,15 +7046,6 @@ void SMESH_Mesh_i::SetNbThreads(int nbThreads){
   _impl->SetNbThreads(nbThreads);
 }
 
-//=============================================================================
-/*!
- * \brief Set the number of threads for the mesher for a parallel computation
- */
-//=============================================================================
-void SMESH_Mesh_i::SetMesherNbThreads(int nbThreads){
-  _impl->SetMesherNbThreads(nbThreads);
-}
-
 
 //=============================================================================
 /*!

src/SMESH_SWIG/smeshBuilder.py

@@ -462,6 +462,21 @@ class smeshBuilder( SMESH._objref_SMESH_Gen, object ):
             obj,name = name,obj
         return Mesh(self, self.geompyD, obj, name)
 
+    def ParallelMesh(self, obj, param, nbThreads, name=0):
+        """
+        Create a parallel mesh.
+
+        Parameters:
+            obj: geometrical object for meshing
+            name: the name for the new mesh.
+            param: full mesh parameters
+            nbThreads: Number of threads for parallelisation.
+
+        Returns:
+            an instance of class :class:`ParallelMesh`.
+        """
+        return ParallelMesh(self, self.geompyD, obj, param, nbThreads, name)
+
     def RemoveMesh( self, mesh ):
         """
         Delete a mesh
@@ -1863,31 +1878,6 @@ class Mesh(metaclass = MeshMeta):
                 geom = self.geom
         return self.smeshpyD.Evaluate(self.mesh, geom)
 
-    def ParallelCompute(self, nbThreads, mesherNbThreads=1, geom=0, discardModifs=False, refresh=False):
-        """
-        Parallel computation of the mesh and return the status of the computation
-        The mesh must contains have be constructed using create_parallel_mesh
-
-        Parameters:
-                nbThreads: Number of threads to use for a parallel computation
-                geom: geomtrical shape on which mesh data should be computed
-                discardModifs: if True and the mesh has been edited since
-                        a last total re-compute and that may prevent successful partial re-compute,
-                        then the mesh is cleaned before Compute()
-                refresh: if *True*, Object Browser is automatically updated (when running in GUI)
-
-        Returns:
-                True or False
-        """
-        if (nbThreads <= 1):
-            raise ValueError("nbThreads must be strictly greater than 1")
-        if (mesherNbThreads < 1):
-            raise ValueError("nbThreads must be greater than 1")
-
-        self.mesh.SetMesherNbThreads(mesherNbThreads)
-        self.mesh.SetNbThreads(nbThreads)
-        return self.Compute(geom=geom, discardModifs=discardModifs, refresh=refresh)
-
     def Compute(self, geom=0, discardModifs=False, refresh=False):
         """
         Compute the mesh and return the status of the computation
@@ -7511,6 +7501,121 @@ class Mesh(metaclass = MeshMeta):
     pass # end of Mesh class
 
 
+def _copy_netgen_param(dim, local_param, global_param):
+    if dim==1:
+        #TODO: Try to identify why we need to substract 1
+        local_param.NumberOfSegments(int(global_param.GetNbSegPerEdge())-1)
+    elif dim==2:
+        local_param.SetMaxSize(global_param.GetMaxSize())
+        local_param.SetMinSize(global_param.GetMinSize())
+        local_param.SetOptimize(global_param.GetOptimize())
+        local_param.SetFineness(global_param.GetFineness())
+        local_param.SetNbSegPerEdge(global_param.GetNbSegPerEdge())
+        local_param.SetNbSegPerRadius(global_param.GetNbSegPerRadius())
+        local_param.SetGrowthRate(global_param.GetGrowthRate()*0.9)
+        local_param.SetChordalError(global_param.GetChordalError())
+        local_param.SetChordalErrorEnabled(global_param.GetChordalErrorEnabled())
+        local_param.SetUseSurfaceCurvature(global_param.GetUseSurfaceCurvature())
+        local_param.SetUseDelauney(global_param.GetUseDelauney())
+        local_param.SetQuadAllowed(global_param.GetQuadAllowed())
+        local_param.SetWorstElemMeasure(global_param.GetWorstElemMeasure())
+        local_param.SetCheckChartBoundary(global_param.GetCheckChartBoundary())
+        local_param.SetNbThreads(global_param.GetNbThreads())
+    else:
+        local_param.SetMaxSize(global_param.GetMaxSize())
+        local_param.SetMinSize(global_param.GetMinSize())
+        local_param.SetOptimize(global_param.GetOptimize())
+        local_param.SetCheckOverlapping(global_param.GetCheckOverlapping())
+        local_param.SetCheckChartBoundary(global_param.GetCheckChartBoundary())
+        local_param.SetFineness(global_param.GetFineness())
+        local_param.SetNbSegPerEdge(global_param.GetNbSegPerEdge())
+        local_param.SetNbSegPerRadius(global_param.GetNbSegPerRadius())
+        local_param.SetGrowthRate(global_param.GetGrowthRate())
+        local_param.SetNbThreads(global_param.GetNbThreads())
+
+class ParallelMesh(Mesh):
+    """
+    Surcharge on Mesh for parallel computation of a mesh
+    """
+
+    def __init__(self, smeshpyD, geompyD, geom, param, nbThreads, name=0):
+        """
+        Create a parallel mesh.
+
+        Parameters:
+            geom: geometrical object for meshing
+            param: full mesh parameters
+            nbThreads: Number of threads for parallelisation.
+            name: the name for the new mesh.
+
+        Returns:
+            an instance of class :class:`ParallelMesh`.
+        """
+
+        if not isinstance(geom, geomBuilder.GEOM._objref_GEOM_Object):
+            raise ValueError("geom argument must be a geometry")
+
+        if not isinstance(param, NETGENPlugin._objref_NETGENPlugin_Hypothesis):
+            raise ValueError("param must come from NETGENPlugin")
+
+        if nbThreads < 1:
+            raise ValueError("Number of threads must be stricly greater than 1")
+
+        # Splitting geometry into 3D elements and all the 2D/1D into one compound
+        object_solids = geompyD.ExtractShapes(geom, geompyD.ShapeType["SOLID"],
+                                              True)
+
+        solids = []
+        isolid = 0
+        for solid in object_solids:
+            isolid += 1
+            geompyD.addToStudyInFather( geom, solid, 'Solid_{}'.format(isolid) )
+            solids.append(solid)
+
+        faces = []
+        iface = 0
+        for isolid, solid in enumerate(solids):
+            solid_faces = geompyD.ExtractShapes(solid, geompyD.ShapeType["FACE"],
+                                                True)
+            for face in solid_faces:
+                faces.append(face)
+                iface += 1
+                geompyD.addToStudyInFather(solid, face,
+                                           'Face_{}'.format(iface))
+
+        # Creating submesh for edges 1D/2D part
+
+        all_faces = geompyD.MakeCompound(faces)
+        geompyD.addToStudy(all_faces, 'Compound_1')
+        all_faces = geompyD.MakeGlueEdges(all_faces, 1e-07)
+        all_faces = geompyD.MakeGlueFaces(all_faces, 1e-07)
+        geompyD.addToStudy(all_faces, 'global2D')
+
+        super(ParallelMesh, self).__init__(smeshpyD, geompyD, geom, name)
+
+        self.mesh.SetNbThreads(nbThreads)
+
+        self.UseExistingSegments()
+        self.UseExistingFaces()
+
+        algo2d = self.Triangle(geom=all_faces, algo="NETGEN_2D")
+        param2d = algo2d.Parameters()
+
+        _copy_netgen_param(2, param2d, param)
+
+        for solid_id, solid in enumerate(solids):
+            name = "Solid_{}".format(solid_id)
+            self.UseExistingSegments(geom=solid)
+            self.UseExistingFaces(geom=solid)
+            algo3d = self.Tetrahedron(geom=solid, algo="NETGEN_3D_Remote")
+
+            param3d = algo3d.Parameters()
+
+            _copy_netgen_param(3, param3d, param)
+
+    pass # End of ParallelMesh
+
+
 class meshProxy(SMESH._objref_SMESH_Mesh):
     """
     Private class used to compensate change of CORBA API of SMESH_Mesh for backward compatibility

test/SMESH_ParallelCompute.py

@@ -0,0 +1,131 @@
+# contains function to compute a mesh in parallel
+from platform import java_ver
+import sys
+from tkinter import W
+import salome
+
+import time
+
+
+salome.salome_init()
+import salome_notebook
+notebook = salome_notebook.NoteBook()
+
+###
+### GEOM component
+###
+
+import GEOM
+from salome.geom import geomBuilder
+from salome.smesh import smeshBuilder
+import math
+import SALOMEDS
+
+import numpy as np
+
+geompy = geomBuilder.New()
+
+smesh = smeshBuilder.New()
+
+
+def build_seq_mesh(nbox, boxsize, offset):
+    # Create 3D faces
+    boxes = []
+    # First creating all the boxes
+    for i in range(nbox):
+        for j in range(nbox):
+            for k in range(nbox):
+
+                x_orig = i*(boxsize+offset)
+                y_orig = j*(boxsize+offset)
+                z_orig = k*(boxsize+offset)
+
+                tmp_box = geompy.MakeBoxDXDYDZ(boxsize, boxsize, boxsize)
+
+                if not i == j == k == 0:
+                    box = geompy.MakeTranslation(tmp_box, x_orig,
+                                                 y_orig, z_orig)
+                else:
+                    box = tmp_box
+
+                geompy.addToStudy(box, 'box_{}:{}:{}'.format(i, j, k))
+
+                boxes.append(box)
+
+    # Create fuse of all boxes
+    all_boxes = geompy.MakeCompound(boxes)
+    geompy.addToStudy(all_boxes, 'Compound_1')
+
+    # Removing duplicates faces and edges
+    all_boxes = geompy.MakeGlueFaces(all_boxes, 1e-07)
+    geompy.addToStudy(all_boxes, 'Glued_Faces_1')
+
+    all_boxes = geompy.MakeGlueEdges(all_boxes, 1e-07)
+    geompy.addToStudy(all_boxes, 'rubik_cube')
+
+
+    # Building sequetial mesh
+    print("Creating mesh")
+    all_box_mesh = smesh.Mesh(all_boxes, "seq_mesh")
+
+    print("Adding algo")
+    algo3d = all_box_mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
+
+    netgen_parameters = algo3d.Parameters()
+    netgen_parameters.SetMaxSize(34.641)
+    netgen_parameters.SetMinSize(0.141421)
+    netgen_parameters.SetOptimize(1)
+    netgen_parameters.SetCheckOverlapping(0)
+    netgen_parameters.SetCheckChartBoundary(0)
+    netgen_parameters.SetFineness(5)
+    netgen_parameters.SetNbSegPerEdge(16*(boxsize//100))
+    netgen_parameters.SetNbSegPerRadius(1.5)
+    netgen_parameters.SetGrowthRate(0.15)
+    netgen_parameters.SetChordalError(-1)
+    netgen_parameters.SetChordalErrorEnabled(0)
+    netgen_parameters.SetUseSurfaceCurvature(1)
+    netgen_parameters.SetQuadAllowed(0)
+    netgen_parameters.SetCheckOverlapping(False)
+    netgen_parameters.SetNbThreads(2)
+
+    return all_boxes, all_box_mesh, netgen_parameters
+
+def run_test(nbox=2, boxsize=100):
+    """ Run sequential mesh and parallel version of it
+
+    nbox: NUmber of boxes
+    boxsize: Size of each box
+    """
+    geom, seq_mesh, netgen_parameters = build_seq_mesh(nbox, boxsize, 0)
+
+    par_mesh = smesh.ParallelMesh(geom, netgen_parameters, 6, name="par_mesh")
+
+    start = time.monotonic()
+    is_done = seq_mesh.Compute()
+    assert is_done
+    stop = time.monotonic()
+    time_seq = stop-start
+
+    start = time.monotonic()
+    is_done = par_mesh.Compute()
+    assert is_done
+    stop = time.monotonic()
+    time_par = stop-start
+
+    print("  Tetrahedron: ", seq_mesh.NbTetras(), par_mesh.NbTetras())
+    print("  Triangle: ", seq_mesh.NbTriangles(), par_mesh.NbTriangles())
+    print("  edge: ", seq_mesh.NbEdges(), par_mesh.NbEdges())
+
+    assert  par_mesh.NbTetras() > 0
+    assert  par_mesh.NbTriangles() > 0
+    assert  par_mesh.NbEdges() > 0
+
+    print("Time elapsed (seq, par): ", time_seq, time_par)
+
+def main():
+    nbox = 2
+    boxsize = 100
+    run_test(nbox, boxsize)
+
+main()
+

test/netgen_runner.py

@@ -22,7 +22,7 @@ import medcoupling as mc
 
 def create_param_file(param_file):
     """ Create a parameter file for runner """
-    param="""1
+    param = """1
 34.64
 0.14
 16
@@ -43,8 +43,10 @@ def create_param_file(param_file):
 0
 0
 2
+2
 0
 
+0
 0
 0"""
     with open(param_file, "w") as ffile:
@@ -57,69 +59,71 @@ def test_netgen3d():
     box = geompy.MakeBoxDXDYDZ(200, 200, 200)
 
     geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True)
-    Groupe_1 = geompy.CreateGroup(box, geompy.ShapeType["FACE"])
+    groupe_1 = geompy.CreateGroup(box, geompy.ShapeType["FACE"])
-    geompy.UnionIDs(Groupe_1, [3, 13, 23, 27, 31, 33])
+    geompy.UnionIDs(groupe_1, [3, 13, 23, 27, 31, 33])
 
-    # TODO: useful ?
+    [_, _, _, _, _, _, groupe_1] = geompy.GetExistingSubObjects(box, False)
-    [_, _, _, _, _, _, Groupe_1] = geompy.GetExistingSubObjects(box, False)
 
     # Creating 2D mesh
-    NETGEN_2D_Parameters_1 = smesh.CreateHypothesisByAverageLength(
+    netgen_2d_parameters_1 = smesh.CreateHypothesisByAverageLength(
         'NETGEN_Parameters_2D', 'NETGENEngine', 34.641, 0)
-    Mesh2D = smesh.Mesh(Groupe_1, 'Maillage_1')
+    mesh_2d = smesh.Mesh(groupe_1, 'Maillage_1')
-    status = Mesh2D.AddHypothesis(Groupe_1, NETGEN_2D_Parameters_1)
+    mesh_2d.AddHypothesis(groupe_1, netgen_2d_parameters_1)
-    NETGEN_1D_2D = Mesh2D.Triangle(algo=smeshBuilder.NETGEN_1D2D)
+    mesh_2d.Triangle(algo=smeshBuilder.NETGEN_1D2D)
-    isDone = Mesh2D.Compute()
+    is_done = mesh_2d.Compute()
-    smesh.SetName(Mesh2D, 'Maillage_1')
+    assert is_done
+    smesh.SetName(mesh_2d, 'Maillage_1')
 
-    # tmp_dir = tempfile.mkdtemp()
     with tempfile.TemporaryDirectory() as tmp_dir:
-      mesh_file = path.join(tmp_dir, "mesh.med")
+        mesh_file = path.join(tmp_dir, "mesh.med")
-      shape_file = path.join(tmp_dir, "shape.step")
+        shape_file = path.join(tmp_dir, "shape.step")
-      param_file = path.join(tmp_dir, "param.txt")
+        param_file = path.join(tmp_dir, "param.txt")
-      output_mesh = path.join(tmp_dir, "mesh3D.med")
+        output_mesh = path.join(tmp_dir, "mesh3D.med")
 
-      print("Running in folder: ", tmp_dir)
+        print("Running in folder: ", tmp_dir)
-      create_param_file(param_file)
+        create_param_file(param_file)
 
-      Mesh2D.ExportMED(mesh_file, 0, 41, 1, Mesh2D, 1, [], '', -1, 1)
+        mesh_2d.ExportMED(mesh_file, 0, 41, 1, mesh_2d, 1, [], '', -1, 1)
-      geompy.ExportSTEP(box, shape_file, GEOM.LU_METER)
+        geompy.ExportSTEP(box, shape_file, GEOM.LU_METER)
 
-      runner = path.join("${NETGENPLUGIN_ROOT_DIR}",
+        runner = path.join("${NETGENPLUGIN_ROOT_DIR}",
-                          "bin",
+                           "bin",
-                          "salome",
+                           "salome",
-                          "NETGENPlugin_Runner")
+                           "NETGENPlugin_Runner")
 
-      cmd = "{runner} NETGEN3D {mesh_file} {shape_file} "\
+        if sys.platform == 'win32':
-            "{param_file} NONE 2 NONE {output_mesh}"\
+            runner += ".exe"
-            .format(runner=runner,
-                    mesh_file=mesh_file,
-                    shape_file=shape_file,
-                    param_file=param_file,
-                    output_mesh=output_mesh)
-      print(cmd)
-      subprocess.check_call(cmd, shell=True)
 
-      meshRead = mc.ReadUMeshFromFile (output_mesh, "MESH", 0)
+        cmd = "{runner} NETGEN3D {mesh_file} {shape_file} "\
+              "{param_file} NONE NONE {output_mesh}"\
+              .format(runner=runner,
+                      mesh_file=mesh_file,
+                      shape_file=shape_file,
+                      param_file=param_file,
+                      output_mesh=output_mesh)
+        print(cmd)
+        subprocess.check_call(cmd, shell=True)
 
-      nbTetras = meshRead.getNumberOfCellsWithType(mc.NORM_TETRA4)
+        mesh_read = mc.ReadUMeshFromFile(output_mesh, "MESH", 0)
-      nbPoints = meshRead.getNumberOfNodes()
 
-      meshRead = mc.ReadUMeshFromFile (output_mesh, "MESH", -1)
+        nb_tetras = mesh_read.getNumberOfCellsWithType(mc.NORM_TETRA4)
-      nbTriangles = meshRead.getNumberOfCellsWithType(mc.NORM_TRI3)
+        nb_points = mesh_read.getNumberOfNodes()
 
-      meshRead = mc.ReadUMeshFromFile (output_mesh, "MESH", -2)
+        mesh_read = mc.ReadUMeshFromFile(output_mesh, "MESH", -1)
-      nbSegments = meshRead.getNumberOfCellsWithType(mc.NORM_SEG2)
+        nb_triangles = mesh_read.getNumberOfCellsWithType(mc.NORM_TRI3)
 
-      print("Nb Tetras:", nbTetras)
+        mesh_read = mc.ReadUMeshFromFile(output_mesh, "MESH", -2)
-      print("Nb Triangles:", nbTriangles)
+        nb_segments = mesh_read.getNumberOfCellsWithType(mc.NORM_SEG2)
-      print("Nb Segments:", nbSegments)
-      print("Nb Points:", nbPoints)
 
-      assert(nbPoints > 0)
+        print("Nb Tetras:", nb_tetras)
-      assert(nbSegments > 0)
+        print("Nb Triangles:", nb_triangles)
-      assert(nbTriangles > 0)
+        print("Nb Segments:", nb_segments)
-      assert(nbTetras > 0)
+        print("Nb Points:", nb_points)
+
+        assert nb_points > 0
+        assert nb_segments > 0
+        assert nb_triangles > 0
+        assert nb_tetras > 0
 
 if __name__ == "__main__":
     test_netgen3d()

test/tests.set

@@ -64,6 +64,7 @@ SET(BAD_TESTS
   SMESH_test4.py
   SMESH_create_dual_mesh_adapt.py
   netgen_runner.py
+  SMESH_ParallelCompute.py
   )
 IF(NOT WIN32)
   LIST(APPEND BAD_TESTS