2 directions for boundary faces

src/genmesh.py

@@ -2,52 +2,40 @@ import os
 import subprocess
 import multiprocessing
 
-src = os.getcwd()
+def salome(src_path, build_path, coefficient, direction):
-build = os.path.join(src, "../build")
+    subprocess.run(["salome", "start", "-t", src_path, "args:{},{}".format(build_path, coefficient, direction)])
 
-if not os.path.exists(build):
+if __name__ == "__main__":
-    os.makedirs(build)
+    # Get main paths
+    project = os.getcwd()
+    src = os.path.join(project, "src")
+    build = os.path.join(project, "build")
 
-###
+    if not os.path.exists(build):
+        os.makedirs(build)
 
-alpha = []
+    ###
+    processes = []
+    structures = ["simple-cubic"] #, "bc-cubic", "fc-cubic"]
+    directions = ["001", "100"]
+    coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ]
 
-for n in range(0.01, 0.13 + 0.01, 0.01):
+    for structure in structures:
-    alpha.append(n)
+        for direction in directions:
+            for coefficient in coefficients:
+                src_path = os.path.join(src, "{}.py".format(structure))
+                build_path = os.path.join(build, 
+                    structure, 
+                    "direction-{}".format(direction), 
+                    "alpha-{}".format(coefficient))
                 
-# Structures
+                if not os.path.exists(build_path):
-simpleCubic = os.path.join(src, "simple-cubic/main.py")
+                    os.makedirs(build_path)
-# Body-centered cubic
-#bcCubic = os.path.join(path, "bc-cubic/main.py")
-# Face-centered cubic
-#fcCubic = os.path.join(path, "fc-cubic/main.py")
 
-###
+                print("starting process")
+                p = multiprocessing.Process(target = salome, args = (src_path, build_path, coefficient, direction))
+                processes.append(p)
+                p.start()
 
-processes = []
+    for process in processes:
-structure = ["simple-cubic"] #, "bc-cubic", "fc-cubic"]
+        process.join()
-direction = ["1"]
-
-def salome(src_path, build_path, arg):
-    subprocess.run(["salome", "start", "-t", src_path, "args:{},{}".format(build_path, arg)])
-
-for s in structure:
-    s_path = os.path.join(build, s)
-    
-    for d in direction:
-        d_path = os.path.join(s_path, d)
-
-        for c in alpha:
-            src_path = os.path.join(src, "%s/main.py" % s)
-            build_path = os.path.join(d_path, str(c))
-            
-            if not os.path.exists(build_path):
-                os.makedirs(build_path)
-
-            print("starting process")
-            p = multiprocessing.Process(target = salome, args = (src_path, build_path, c))
-            processes.append(p)
-            p.start()
-
-        for process in processes:
-            process.join()

src/simple-cubic/geometry.py

@@ -1,83 +0,0 @@
-import GEOM
-from salome.geom import geomBuilder
-geompy = geomBuilder.New()
-
-import math
-
-def create(alpha):
-    # xyz axes
-    axes = [
-        geompy.MakeVectorDXDYDZ(1, 0, 0),
-        geompy.MakeVectorDXDYDZ(0, 1, 0),
-        geompy.MakeVectorDXDYDZ(0, 0, 1)
-    ]
-    
-    # Main box
-    box = geompy.MakeBoxDXDYDZ(2 * math.sqrt(2), 2 * math.sqrt(2), 2)
-    box = geompy.MakeRotation(box, axes[2], 45 * math.pi / 180.0)
-    box = geompy.MakeTranslation(box, 2, 0, 0)
-    
-    vtx = [
-        geompy.MakeVertex(2, 0, 0),
-        geompy.MakeVertex(2, 2, 0),
-        geompy.MakeVertex(2, 2, 2)
-    ]
-    
-    line = geompy.MakeLineTwoPnt(vtx[1], vtx[2])
-    
-    # Spheres for cutting
-    sphere = geompy.MakeSpherePntR(vtx[0], 1 / (1 - alpha))
-    sphere = geompy.MakeMultiTranslation1D(sphere, axes[1], 2, 3)
-    cut = geompy.MakeCutList(box, [sphere], True)
-    
-    sphere2 = geompy.MakeTranslation(sphere, 0, 0, 2)
-    cut2 = geompy.MakeCutList(cut, [sphere2], True)
-    
-    sphere3 = geompy.MakeRotation(sphere, line, 90 * math.pi / 180.0)
-    cut3 = geompy.MakeCutList(cut2, [sphere3], True)
-    
-    sphere4 = geompy.MakeTranslation(sphere3, 0, 0, 2)
-    
-    # Main geometry
-    Pore = geompy.MakeCutList(cut3, [sphere4], True)
-    geompy.addToStudy(Pore, "PORE")
-
-    #
-    #   D /\ B  A(1, 1, 0)  \vec{n}(1, 1, 0)
-    #    /  \   B(3, 3, 0)  \vec{n}(1, 1, 0)
-    #    \  /   C(3, 1, 0)  \vec{n}(-1, 1, 0)
-    #   A \/ C  D(1, 3, 0)  \vec{n}(-1, 1, 0)
-    #
-
-    # Prepare faces
-    vtx.append(geompy.MakeVertex(2, 2, 2))
-    vtx.append(geompy.MakeVertexWithRef(vtx[3], 0, 0, 1))
-    vec2 = geompy.MakeVector(vtx[3], vtx[4])
-    plane = geompy.MakePlane(vtx[3], vec2, 5)
-    common1 = geompy.MakeCommonList([Pore, plane], True)
-    plane = geompy.MakeTranslation(plane, 0, 0, -2)
-    common2 = geompy.MakeCommonList([Pore, plane], True)
-    
-    # TODO: make 4 planes A B C D
-    Pore = geompy.MakeFillet(Pore, 0.1, geompy.ShapeType["EDGE"], [24, 27, 35])
-
-
-    # Main groups (inlet, outlet, wall)
-    inlet = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"])
-    gip = geompy.GetInPlace(Pore, common1, True)
-    faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
-    geompy.UnionList(inlet, faces)
-    
-    outlet = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"])
-    gip = geompy.GetInPlace(Pore, common2, True)
-    faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"]) 
-    geompy.UnionList(outlet, faces)
-    
-    PoreGroup = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"])
-    faces = geompy.SubShapeAllIDs(Pore, geompy.ShapeType["FACE"]) 
-    geompy.UnionIDs(PoreGroup, faces)
-    
-    wall = geompy.CutListOfGroups([PoreGroup], [inlet, outlet])
-    
-    
-    return Pore, [inlet, outlet, wall]

src/simple-cubic/main.py

@@ -1,35 +0,0 @@
-import os, sys
-
-import salome
-salome.salome_init()
-
-import geometry, mesh
-
-#alpha = [ 0.1, 0.15, 0.2 ]
-build_path = str(sys.argv[1])
-coef = float(sys.argv[2])
-
-#for coef in alpha:
-print("alpha = {}".format(coef))
-print("Building geometry ...")
-
-Pore, bc = geometry.create(coef)
-#geometry.geompy.addToStudy(Pore, 'Pore {}'.format(coef))
-#geometry.geompy.addToStudyInFather(Pore, bc[0], "INLET")
-print("Building mesh ...")
-
-PoreMesh = mesh.create(Pore, bc)
-isDone = PoreMesh.Compute()
-
-status = "Succesfully" if isDone else "Mesh is not computed"
-print(status)
-
-try:
-    filename = os.path.join(build_path, "mesh.unv")
-    PoreMesh.ExportUNV(filename)
-    pass
-except:
-    print('ExportUNV() failed. Invalid file name?')
-
-if salome.sg.hasDesktop():
-    salome.sg.updateObjBrowser()

src/simple-cubic/mesh.py

@@ -1,34 +0,0 @@
-import SMESH, SALOMEDS
-from salome.smesh import smeshBuilder
-smesh = smeshBuilder.New()
-
-
-def create(geomObj, bc):
-    mesh = smesh.Mesh(geomObj)
-    netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
-
-    param = netgen.Parameters()
-    param.SetSecondOrder( 0 )
-    param.SetOptimize( 1 )
-    param.SetChordalError( -1 )
-    param.SetChordalErrorEnabled( 0 )
-    param.SetUseSurfaceCurvature( 1 )
-    param.SetFuseEdges( 1 )
-    param.SetCheckChartBoundary( 0 )
-    param.SetMinSize( 0.01 )
-    param.SetMaxSize( 0.1 )
-    param.SetFineness( 4 )
-    #param.SetGrowthRate( 0.1 )
-    #param.SetNbSegPerEdge( 5 )
-    #param.SetNbSegPerRadius( 10 )
-    param.SetQuadAllowed( 0 )
-
-    vlayer = netgen.ViscousLayers(0.025, 2, 1.1, [],
-        1, smeshBuilder.NODE_OFFSET)
-    
-    mesh.GroupOnGeom(bc[0], 'inlet', SMESH.FACE)
-    mesh.GroupOnGeom(bc[1], 'outlet', SMESH.FACE)
-    mesh.GroupOnGeom(bc[2], 'wall', SMESH.FACE)
-
-    return mesh
-

src/simpleCubic.py

@@ -15,6 +15,19 @@ class simpleCubic:
         salome.salome_init()
 
     def geometryCreate(self, alpha):
+        """
+        Create the simple cubic geometry.
+
+        Parameters:
+            alpha (float): Sphere intersection parameter which used for cutting spheres from box.
+                
+                Radius = R_0 / (1 - alpha)
+                Should be from 0.01 to 0.13
+
+        Returns:
+            Configured geometry.
+        """
+
         geompy = geomBuilder.New()
         
         #
@@ -58,15 +71,44 @@ class simpleCubic:
         
         geompy.addToStudy(self.geometry, self.name)
 
+        return self.geometry
+
     def boundaryCreate(self, direction):
-        geompy = geomBuilder.New()
+        """
+        Create the boundary faces from the geometry.
+
+        Parameters:
+            direction (str): Direction of the flow.
+
+                '001' for the flow with normal vector (0, 0, -1) to face.
+                '100' for the flow with normal vector (-1, 0, 0) to face.
+
+        Returns:
+            boundary (dict):
+
+            {
+                "inlet": <GEOM._objref_GEOM_Object>,
+                "outlet": <GEOM._objref_GEOM_Object>,
+                "symetryPlane": <GEOM._objref_GEOM_Object>,
+                "wall": <GEOM._objref_GEOM_Object>
+            }
+
+        """
         #
-        #   D /\ B  A(1, 1, 0)  \vec{n}(1, 1, 0)
+        #      _____    z      |
-        #    /  \   B(3, 3, 0)  \vec{n}(1, 1, 0)
+        #     //////|   |      | flow
-        #    \  /   C(3, 1, 0)  \vec{n}(-1, 1, 0)
+        #    ////// |   |___y  f
-        #   A \/ C  D(1, 3, 0)  \vec{n}(-1, 1, 0)
+        #    |    | /   /     
+        #    |____|/   /x      direction [0, 0, 1]
+        #
+        #      _____    z        f
+        #     /    /|   |       / flow
+        #    /____/ |   |___y  /
+        #    |||||| /   /      
+        #    ||||||/   /x      direction [1, 0, 0]
         #
         
+        geompy = geomBuilder.New() 
         center = geompy.MakeVertex(2, 2, 1)
         rot = [0, 0, 45]
 
@@ -78,7 +120,9 @@ class simpleCubic:
         
         elif direction == "100":
             norm = geompy.MakeVector(center, 
-                geompy.MakeVertexWithRef(center, 1, 0, 0))
+                geompy.MakeVertexWithRef(center, 
+                    math.cos((90 + rot[2]) * math.pi / 180.0), 
+                    -math.sin((90 + rot[2]) * math.pi / 180.0), 0))
             vstep = math.sqrt(2)
             hstep = 1
         
@@ -102,14 +146,13 @@ class simpleCubic:
         box = geompy.MakeBoxDXDYDZ(2 * math.sqrt(2), 2 * math.sqrt(2), 2)
         box = geompy.MakeRotation(box, axes[2], 45 * math.pi / 180.0)
         box = geompy.MakeTranslation(box, 2, 0, 0)
-        planes = geompy.ExtractShapes(box, 
+        planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True)
-            geompy.ShapeType["FACE"], True)
 
         vplanes = []
         hplanes = []
         for plane in planes:
             planeNorm = geompy.GetNormal(plane)
-            angle = geompy.GetAngle(planeNorm, norm)
+            angle = abs(geompy.GetAngle(planeNorm, norm))
 
             if angle == 0 or angle == 180:
                 vplanes.append(plane)
@@ -117,26 +160,29 @@ class simpleCubic:
             else:
                 hplanes.append(plane)
         
-        zvplane1 = vplanes[0]
-        zvplane2 = vplanes[1]
-
         if direction == "001":
-            if geompy.GetPosition(zvplane1)[3] > geompy.GetPosition(zvplane1)[3]:
+            z1 = geompy.GetPosition(vplanes[0])[3]
-                inletplane = zvplane1
+            z2 = geompy.GetPosition(vplanes[1])[3]
-                outletplane = zvplane2
+
+            if z1 > z2:
+                inletplane = vplanes[0]
+                outletplane = vplanes[1]
             
             else:
-                inletplane = zvplane2
+                inletplane = vplanes[1]
-                outletplane = zvplane1
+                outletplane = vplanes[0]
 
         elif direction == "100":
-            if geompy.GetPosition(zvplane1)[1] > geompy.GetPosition(zvplane1)[1]:
+            x1 = geompy.GetPosition(vplanes[0])[1]
-                inletplane = zvplane1
+            x2 = geompy.GetPosition(vplanes[1])[1]
-                outletplane = zvplane2
+
+            if x1 > x2:
+                inletplane = vplanes[0]
+                outletplane = vplanes[1]
             
             else:
-                inletplane = zvplane2
+                inletplane = vplanes[1]
-                outletplane = zvplane1
+                outletplane = vplanes[0]
 
 
         # inlet and outlet
@@ -146,9 +192,6 @@ class simpleCubic:
         common2 = geompy.MakeCommonList([self.geometry, outletplane], True)
         outlet = createGroup(common2, "outlet")
 
-        geompy.addToStudy(inletplane, "inletplane")
-        geompy.addToStudy(outletplane, "outletplane")
-        
         # symetryPlane(s)
         symetryPlane = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], "symetryPlane")
         
@@ -173,10 +216,35 @@ class simpleCubic:
 
         return self.boundary
 
-    def meshCreate(self):
+    def meshCreate(self, fineness, viscousLayers=None):
+        """
+        Creates a mesh from a geometry.
+
+        Parameters:
+            fineness (int): Fineness of mesh.
+                
+                0 - Very coarse,
+                1 - Coarse,
+                2 - Moderate,
+                3 - Fine,
+                4 - Very fine.
+
+            viscousLayers (dict or None): Defines viscous layers for mesh. 
+                By default, inlets and outlets specified without layers.
+                
+                {
+                    "thickness": float,
+                    "number": int,
+                    "stretch": float
+                }
+
+        Returns:
+            Configured instance of class <SMESH.SMESH_Mesh>, containig the parameters and boundary groups.
+
+        """
         smesh = smeshBuilder.New()
 
-        mesh = smesh.Mesh(geomObj)
+        mesh = smesh.Mesh(self.geometry)
         netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
 
         param = netgen.Parameters()
@@ -189,25 +257,28 @@ class simpleCubic:
         param.SetCheckChartBoundary( 0 )
         param.SetMinSize( 0.01 )
         param.SetMaxSize( 0.1 )
-        param.SetFineness( 4 )
+        param.SetFineness(fineness)
         #param.SetGrowthRate( 0.1 )
         #param.SetNbSegPerEdge( 5 )
         #param.SetNbSegPerRadius( 10 )
         param.SetQuadAllowed( 0 )
         
-        vlayer = netgen.ViscousLayers(0.025, 2, 1.1, [],
+        if not viscousLayers is None:
-            1, smeshBuilder.NODE_OFFSET)
+            vlayer = netgen.ViscousLayers(viscousLayers["thickness"], 
+                                          viscousLayers["number"], 
+                                          viscousLayers["stretch"], 
+                                          [self.boundary["inlet"], self.boundary["outlet"]],
+                                          1, smeshBuilder.NODE_OFFSET)
         
-        mesh.GroupOnGeom(self.boundary["inlet"], "inlet", SMESH.FACE)
+        for name, boundary in self.boundary.items():
-        mesh.GroupOnGeom(self.boundary["outlet"], "outlet", SMESH.FACE)
+            mesh.GroupOnGeom(boundary, name, SMESH.FACE)
-        mesh.GroupOnGeom(self.boundary["symetryPlane"], "symetryPlane", SMESH.FACE)
-        mesh.GroupOnGeom(self.boundary["wall"], "wall", SMESH.FACE)
 
         self.mesh = mesh
 
         return self.mesh
 
     def meshCompute(self):
+        """Compute the mesh."""
         status = self.mesh.Compute()
         
         if status:
@@ -217,6 +288,12 @@ class simpleCubic:
             print("Mesh is not computed.")
 
     def meshExport(self, path):
+        """
+        Export the mesh in a file in UNV format.
+
+        Parameters:
+            path (string): full path to the expected directory.
+        """
         exportpath = os.path.join(path, "{}.unv".format(self.name))
 
         try:
@@ -230,12 +307,15 @@ if __name__ == "__main__":
     alpha = float(sys.argv[2])
     direction = str(sys.argv[3])
 
-    #salome.salome_init()
     sc = simpleCubic("simpleCubic-{}-{}".format(direction, alpha))
     sc.geometryCreate(alpha)
     sc.boundaryCreate(direction)
-    #sc.meshCreate()
+    sc.meshCreate(2, {
-    #sc.meshCompute()
+        "thickness": 0.02,
+        "number": 2,
+        "stretch": 1.1
+    })
+    sc.meshCompute()
     #sc.meshExport(build)
 
     if salome.sg.hasDesktop():

src/structure.gv

@@ -1,14 +0,0 @@
-digraph structure {
-    node [shape=plaintext]
-    s1 [label="src/simple-cubic/geometry.py"]
-    s2 [label="src/simple-cubic/mesh.py"]
-    s3 [label="src/simple-cubic/main.py"]
-    s1 -> s3
-    s2 -> s3
-    s4 [label="src/genmesh.py"]
-    s3 -> s4
-    b1 [label="build/simple-cubic/direction-1/0.1/"]
-    s4 -> b1 [label="[alpha=0.1, direction=1, mesh.unv]"]
-    s5 [label="src/prefoam.py"]
-    s5 -> b1 [label="[0, constant, system (/src/foam/*)]"]
-}