Destruction [4]: Mesh

run.py

@@ -17,10 +17,13 @@ from src import salome_utils
 from src import foam_utils
 
 def createTasks():
+    ###
+    #   Control variables
+    ##
     structures = [
-        "simpleCubic", 
+        "simple", 
-        "bodyCenteredCubic", 
+        #"bodyCentered", 
-        "faceCenteredCubic"
+        #"faceCentered"
     ]
     directions = [
         [1, 0, 0],
@@ -29,19 +32,22 @@ def createTasks():
     ]
     fillet = 1
 
+    ###
+    #   Tasks
+    ##
     Task = namedtuple("Task", ["structure", "coeff", "fillet", "direction", "saveto"])
     tasks = []
 
     for structure in structures:
-        if structure == "simpleCubic":
+        if structure == "simple":
-            #theta = [c * 0.01 for c in range(1, 28 + 1)]
+            theta = [c * 0.01 for c in range(1, 28 + 1)]
-            theta = [ 0.01, 0.28 ]
+            #theta = [ 0.01, 0.28 ]
 
-        elif structure == "faceCenteredCubic":
+        elif structure == "faceCentered":
             #theta = [c * 0.01 for c in range(1, 13 + 1)]
             theta = [ 0.01, 0.13 ]
 
-        elif structure == "bodyCenteredCubic":
+        elif structure == "bodyCentered":
             #theta = [c * 0.01 for c in range(1, 18 + 1)]
             theta = [ 0.01, 0.13, 0.14, 0.18 ]
 

samples/__init__.py

@@ -11,66 +11,101 @@ LOG = os.path.join(ROOT, "logs")
 
 import salome
 
-from simpleCubic import simpleCubic
+from simple import simpleCubic, simpleHexagonalPrism
-from faceCenteredCubic import faceCenteredCubic
+from faceCentered import faceCenteredCubic, faceCenteredHexagonalPrism
-from bodyCenteredCubic import bodyCenteredCubic
+from bodyCentered import bodyCenteredCubic, bodyCenteredHexagonalPrism
 
 from src import geometry_utils
 from src import mesh_utils
 
-def genMesh(stype, theta, fillet, flowdirection, saveto):
+def genMesh(stype, theta, fillet, direction, saveto):
-    _G = globals()
 
-    structure = _G.get(stype)
+    logging.info("""genMesh: 
+    structure type:\t{}
+    coefficient:\t{}
+    fillet:\t{}
+    flow direction:\t{}
+    export path:\t{}""".format(stype, theta, fillet, direction, saveto))
 
-    if structure:
+    params = (theta, fillet, direction)
-        salome.salome_init()
 
-        grains, geometry1, geometry2 = structure(theta, fillet)
+    salome.salome_init()
-        geometry = geometry1
+    
-        
+    ###
-        if flowdirection == [1, 1, 1]:
+    #   Structure and mesh configurations
-            geometry = geometry2
+    ##
-            norm = [-1, 1, 0]
+    if stype == "simple":
-            bcount = 6
+        if direction in [[1, 0, 0], [0, 0, 1]]:
+            structure = simpleCubic
 
-            # initial angle
+        elif direction == [1, 1, 1]:
-            angle = math.pi / 6
+            structure = simpleHexagonalPrism
-            v1 = Quaternion(axis = norm, angle = math.pi / 2).rotate(flowdirection)
-            normvec = Quaternion(axis = flowdirection, angle = angle).rotate(v1)
-            direction = [1, 1, 1]
 
-        if flowdirection == [1, 0, 0]:
+    elif stype == "faceCentered":
-            normvec = [0, 0, 1]
+        if direction in [[1, 0, 0], [0, 0, 1]]:
-            bcount = 4
+            structure = faceCenteredCubic
-            direction = [1, 1, 0]
 
-        if flowdirection == [0, 0, 1]:
+        elif direction == [1, 1, 1]:
-            normvec = [1, 1, 0]
+            structure = faceCenteredHexagonalPrism
-            bcount = 4
-            direction = [0, 0, 1]
-        
-        #
-        geometry = geometry_utils.geompy.RemoveExtraEdges(geometry, False) 
 
-        #
+    elif stype == "bodyCentered":
-        boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains)
+        if direction in [[1, 0, 0], [0, 0, 1]]:
+            structure = bodyCenteredCubic
 
-        fineness = 1
+        elif direction == [1, 1, 1]:
-        viscousLayers = {
+            structure = bodyCenteredHexagonalPrism
-            "thickness": 0.0001,
+    
-            "number": 3,
+    ###
-            "stretch": 1.2
+    #   Shape
-        }
+    ##
-        mesh = mesh_utils.meshCreate(geometry, boundary, fineness, viscousLayers)
+    geompy = geometry_utils.getGeom()
-        mesh_utils.meshCompute(mesh)
+    shape, groups = structure(*params)
+    [length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06)
 
-        mesh_utils.meshExport(mesh, saveto)
+    logging.info("""shape:
+    surface area:\t{}
+    volume:\t{}""".format(surfaceArea, volume))
+    
+    ###
+    #   Mesh
+    ##
+    fineness = 1
+    parameters = mesh_utils.Parameters(
+        minSize = 0.001,
+        maxSize = 0.1,
+        growthRate = 0.1,
+        nbSegPerEdge = 5,
+        nbSegPerRadius = 10,
+        chordalErrorEnabled = False,
+        chordalError = -1,
+        secondOrder = False,
+        optimize = True,
+        quadAllowed = False,
+        useSurfaceCurvature = True,
+        fuseEdges = True,
+        checkChartBoundary = False
+    )
+    
+    facesToIgnore = []
+    for group in groups:
+        if group.GetName() in ["inlet", "outlet"]:
+            facesToIgnore.append(group)
 
-        salome.salome_close()
+    viscousLayers = mesh_utils.ViscousLayers(
+        thickness = 0.001,
+        numberOfLayers = 3,
+        stretchFactor = 1.2,
+        isFacesToIgnore = True,
+        facesToIgnore = facesToIgnore,
+        extrusionMethod = mesh_utils.smeshBuilder.NODE_OFFSET
+    )
+    
+    mesh = mesh_utils.meshCreate(shape, groups, fineness, parameters, viscousLayers)
+    mesh_utils.meshCompute(mesh)
 
-    else:
+    mesh_utils.meshExport(mesh, saveto)
-        raise Exception("Unknown sample function")
+
+    salome.salome_close()
 
 
 if __name__ == "__main__":
@@ -89,12 +124,5 @@ if __name__ == "__main__":
     flowdirection = [int(coord) for coord in sys.argv[4]]
     saveto = str(sys.argv[5])
 
-    logging.info("""genMesh: 
-    structure type:\t{}
-    coefficient:\t{}
-    fillet:\t{}
-    flow direction:\t{}
-    export path:\t{}""".format(stype, theta, fillet, flowdirection, saveto))
-
     genMesh(stype, theta, fillet, flowdirection, saveto)
 

samples/bodyCentered.py

@@ -0,0 +1,284 @@
+import salome
+salome.salome_init()
+
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+
+from math import pi, sqrt
+
+def bodyCenteredCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]):
+    
+    ###
+    #   Parameters
+    ##
+    L = 1.0
+    r0 = L * sqrt(3) / 4
+
+    radius = r0 / (1 - theta)
+    xn, yn, zn = 3, 3, 3 
+    
+    length = 2 * r0
+    width = L / 2
+    diag = L * sqrt(2)
+    height = L
+
+    point = []
+    xl = sqrt(diag ** 2 + diag ** 2) * 0.5
+    yw = xl
+    zh = height
+
+    C1, C2 = 0.8, 0.05
+    theta1, theta2 = 0.01, 0.18
+    Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1)
+    filletradius = Cf * (radius - r0)
+    
+    scale = 100
+    oo = geompy.MakeVertex(0, 0, 0)
+    spos1 = (0, 0, 0)
+    spos2 = (0, 0, 0)
+
+    ###
+    #   Bounding box
+    ##
+    if direction == [1, 0, 0]:
+        sk = geompy.Sketcher3D()
+        sk.addPointsAbsolute(xl, 0, 0)
+        sk.addPointsAbsolute(0, yw, 0)
+        sk.addPointsAbsolute(0, yw, zh)
+        sk.addPointsAbsolute(xl, 0, zh)
+        sk.addPointsAbsolute(xl, 0, 0)
+
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, diag)
+
+    elif direction == [0, 0, 1]:
+        sk = geompy.Sketcher3D()
+        sk.addPointsAbsolute(0, yw, 0)
+        sk.addPointsAbsolute(xl, 0, 0)
+        sk.addPointsAbsolute(2 * xl, yw, 0)
+        sk.addPointsAbsolute(xl, 2 * yw, 0)
+        sk.addPointsAbsolute(0, yw, 0)
+
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, zh)
+
+    else:
+        raise Exception("The direction is not implemented")
+    
+    inletface = geompy.MakeScaleTransform(inletface, oo, scale)
+    cubic = geompy.MakeScaleTransform(cubic, oo, scale)
+
+    faces = geompy.ExtractShapes(cubic, geompy.ShapeType["FACE"], False)
+    symetryface = []
+
+    for face in faces:
+        norm = geompy.GetNormal(face)
+        angle = round(geompy.GetAngle(norm, vecflow), 0)
+
+        if (angle == 0 or angle == 180) and not face == inletface:
+            outletface = face
+        
+        else:
+            symetryface.append(face)
+    
+    ###
+    #   Grains
+    ##
+    ox = geompy.MakeVectorDXDYDZ(1, 0, 0)
+    oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
+    oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
+    xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
+    xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius)
+    lattice1 = geompy.MakeMultiTranslation2D(grain, ox, L, xn, oy, L, yn)
+    lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn)
+
+    #grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius)
+    #lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
+    #lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn)
+    lattice2 = geompy.MakeTranslation(lattice1, 0.5 * L, 0.5 * L, 0.5 * L)
+    
+    grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
+    grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
+    grains = geompy.MakeFuseList(grains, False, False)
+
+    grains = geompy.MakeScaleTransform(grains, oo, scale)
+
+    if fillet:
+        grains = geompy.MakeFilletAll(grains, filletradius * scale)
+
+    ###
+    #   Groups
+    ##
+    shape = geompy.MakeCutList(cubic, [grains])
+    shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "bodyCenteredCubic")
+
+    sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
+    geompy.UnionIDs(sall,
+        geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
+
+    inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
+    inletshape = geompy.MakeCutList(inletface, [grains])
+    inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
+    geompy.UnionList(inlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
+
+    outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
+    outletshape = geompy.MakeCutList(outletface, [grains])
+    outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
+    geompy.UnionList(outlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
+    
+    symetry = []
+    for (n, face) in enumerate(symetryface):
+        name = "symetry" + str(n)
+        symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
+        symetryshape = geompy.MakeCutList(face, [grains])
+        symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
+        geompy.UnionList(symetry[n], geompy.SubShapeAll(
+            geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
+
+    groups = []
+    groups.append(inlet)
+    groups.append(outlet)
+    groups.extend(symetry)
+    wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
+
+    return shape, groups
+
+def bodyCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]):
+    
+    ###
+    #   Parameters
+    ##
+    L = 1.0
+    r0 = L * sqrt(3) / 4
+
+    radius = r0 / (1 - theta)
+    xn, yn, zn = 3, 3, 3 
+    
+    length = 2 * r0
+    width = L / 2
+    diag = L * sqrt(2)
+    height = diag / 3
+
+    point = []
+    xl, yw, zh = -L / 4, -L / 4, -L / 4
+    point.append((L / 3 + xl, L / 3 + yw, 4 * L / 3 + zh))
+    point.append((L + xl, 0 + yw, L + zh))
+    point.append((4 * L / 3 + xl, L / 3 + yw, L / 3 + zh))
+    point.append((L + xl, L + yw, 0 + zh))
+    point.append((L / 3 + xl, 4 * L / 3 + yw, L / 3 + zh))
+    point.append((0 + xl, L + yw, L + zh))
+    point.append((L / 3 + xl, L / 3 + yw, 4 * L / 3 + zh))
+
+    C1, C2 = 0.8, 0.05
+    theta1, theta2 = 0.01, 0.18
+    Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1)
+    filletradius = Cf * (radius - r0)
+    
+    scale = 100
+    oo = geompy.MakeVertex(0, 0, 0)
+    spos1 = (0, 0, 0)
+    spos2 = (0, 0, 0)
+
+    ###
+    #   Bounding box
+    ## 
+    sk = geompy.Sketcher3D()
+
+    for p in point:
+        sk.addPointsAbsolute(*p)
+    
+    inletface = geompy.MakeFaceWires([sk.wire()], False)
+    vecflow = geompy.GetNormal(inletface)
+    hexagonPrism = geompy.MakePrismVecH(inletface, vecflow, L * sqrt(3))
+
+    inletface = geompy.MakeScaleTransform(inletface, oo, scale)
+    hexagonPrism = geompy.MakeScaleTransform(hexagonPrism, oo, scale)
+
+    faces = geompy.ExtractShapes(hexagonPrism, geompy.ShapeType["FACE"], False)
+    symetryface = []
+
+    for face in faces:
+        norm = geompy.GetNormal(face)
+        angle = round(geompy.GetAngle(norm, vecflow), 0)
+
+        if (angle == 0 or angle == 180) and not face == inletface:
+            outletface = face
+        
+        else:
+            symetryface.append(face)
+    
+    ###
+    #   Grains
+    ##
+    ox = geompy.MakeVectorDXDYDZ(1, 0, 0)
+    oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
+    oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
+    xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
+    xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius)
+    lattice1 = geompy.MakeMultiTranslation2D(grain, ox, L, xn, oy, L, yn)
+    lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn)
+
+    #grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius)
+    #lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
+    #lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn)
+    lattice2 = geompy.MakeTranslation(lattice1, 0.5 * L, 0.5 * L, 0.5 * L)
+    
+    grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
+    grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
+    grains = geompy.MakeFuseList(grains, False, False)
+
+    grains = geompy.MakeScaleTransform(grains, oo, scale)
+
+    if fillet:
+        grains = geompy.MakeFilletAll(grains, filletradius * scale)
+
+    ###
+    #   Groups
+    ##
+    shape = geompy.MakeCutList(hexagonPrism, [grains])
+    shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "bodyCenteredCubic")
+
+    sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
+    geompy.UnionIDs(sall,
+        geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
+
+    inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
+    inletshape = geompy.MakeCutList(inletface, [grains])
+    inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
+    geompy.UnionList(inlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
+
+    outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
+    outletshape = geompy.MakeCutList(outletface, [grains])
+    outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
+    geompy.UnionList(outlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
+    
+    symetry = []
+    for (n, face) in enumerate(symetryface):
+        name = "symetry" + str(n)
+        symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
+        symetryshape = geompy.MakeCutList(face, [grains])
+        symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
+        geompy.UnionList(symetry[n], geompy.SubShapeAll(
+            geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
+
+    groups = []
+    groups.append(inlet)
+    groups.append(outlet)
+    groups.extend(symetry)
+    wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
+
+    return shape, groups
+

samples/faceCentered.py

@@ -0,0 +1,282 @@
+import salome
+salome.salome_init()
+
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+
+from math import pi, sqrt
+
+def faceCenteredCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]):
+    
+    ###
+    #   Parameters
+    ##
+    L = 1.0
+    r0 = L * sqrt(2) / 4
+
+    radius = r0 / (1 - theta)
+    xn, yn, zn = 3, 3, 3 
+    
+    length = 2 * r0
+    width = L / 2
+    diag = L * sqrt(3)
+    height = diag / 3
+
+    point = []
+    xl = sqrt(length ** 2 + length ** 2) * 0.5
+    yw = xl
+    zh = width
+
+    C1, C2 = 0.8, 0.05
+    theta1, theta2 = 0.01, 0.13
+    Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1)
+    filletradius = Cf * (radius - r0)
+    
+    scale = 100
+    oo = geompy.MakeVertex(0, 0, 0)
+    spos1 = (-width * (xn - 1), 0, -width * (zn - 2))
+    spos2 = (-width * xn, 0, -width * (zn - 1))
+
+    ###
+    #   Bounding box
+    ##
+    if direction == [1, 0, 0]:
+        sk = geompy.Sketcher3D()
+        sk.addPointsAbsolute(0, 0, -zh)
+        sk.addPointsAbsolute(-xl, yw, -zh)
+        sk.addPointsAbsolute(-xl, yw, zh)
+        sk.addPointsAbsolute(0, 0, zh)
+        sk.addPointsAbsolute(0, 0, -zh)
+
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, length)
+
+    elif direction == [0, 0, 1]:
+        sk = geompy.Sketcher3D()
+        sk.addPointsAbsolute(0, 0, -zh)
+        sk.addPointsAbsolute(xl, yw, -zh)
+        sk.addPointsAbsolute(0, 2 * yw, -zh)
+        sk.addPointsAbsolute(-xl, yw, -zh)
+        sk.addPointsAbsolute(0, 0, -zh)
+
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, 2 * zh)
+
+    else:
+        raise Exception("The direction is not implemented")
+    
+    inletface = geompy.MakeScaleTransform(inletface, oo, scale)
+    cubic = geompy.MakeScaleTransform(cubic, oo, scale)
+
+    faces = geompy.ExtractShapes(cubic, geompy.ShapeType["FACE"], False)
+    symetryface = []
+
+    for face in faces:
+        norm = geompy.GetNormal(face)
+        angle = round(geompy.GetAngle(norm, vecflow), 0)
+
+        if (angle == 0 or angle == 180) and not face == inletface:
+            outletface = face
+        
+        else:
+            symetryface.append(face)
+    
+    ###
+    #   Grains
+    ##
+    ox = geompy.MakeVectorDXDYDZ(1, 0, 0)
+    oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
+    oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
+    xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
+    xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius)
+    lattice1 = geompy.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn)
+    lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn - 1)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius)
+    lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
+    lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn)
+    
+    grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
+    grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
+    grains = geompy.MakeFuseList(grains, False, False)
+
+    grains = geompy.MakeScaleTransform(grains, oo, scale)
+
+    if fillet:
+        grains = geompy.MakeFilletAll(grains, filletradius * scale)
+
+    ###
+    #   Groups
+    ##
+    shape = geompy.MakeCutList(cubic, [grains])
+    shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "faceCenteredCubic")
+
+    sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
+    geompy.UnionIDs(sall,
+        geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
+
+    inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
+    inletshape = geompy.MakeCutList(inletface, [grains])
+    inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
+    geompy.UnionList(inlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
+
+    outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
+    outletshape = geompy.MakeCutList(outletface, [grains])
+    outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
+    geompy.UnionList(outlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
+    
+    symetry = []
+    for (n, face) in enumerate(symetryface):
+        name = "symetry" + str(n)
+        symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
+        symetryshape = geompy.MakeCutList(face, [grains])
+        symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
+        geompy.UnionList(symetry[n], geompy.SubShapeAll(
+            geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
+
+    groups = []
+    groups.append(inlet)
+    groups.append(outlet)
+    groups.extend(symetry)
+    wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
+
+    return shape, groups
+
+def faceCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]):
+    
+    ###
+    #   Parameters
+    ##
+    L = 1.0
+    r0 = L * sqrt(2) / 4
+
+    radius = r0 / (1 - theta)
+    xn, yn, zn = 3, 3, 3 
+    
+    length = 2 * r0
+    width = L / 2
+    diag = L * sqrt(3)
+    height = diag / 3
+
+    point = []
+    xl, yw, zh = -(xn - 2) * L / 3, -(yn - 2) * L / 3, -(zn - 2) * L / 3
+    point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh))
+    point.append((0 + xl, -width + yw, 0 + zh))
+    point.append((width / 3 + xl, -2 * width / 3 + yw, -2 * width / 3 + zh))
+    point.append((0 + xl, 0 + yw, -width + zh))
+    point.append((-2 * width / 3 + xl, width / 3 + yw, -2 * width / 3 + zh))
+    point.append((-width + xl, 0 + yw, 0 + zh))
+    point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh))
+
+    C1, C2 = 0.8, 0.05
+    theta1, theta2 = 0.01, 0.13
+    Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1)
+    filletradius = Cf * (radius - r0)
+    
+    scale = 100
+    oo = geompy.MakeVertex(0, 0, 0)
+    spos1 = (-width * (xn - 1), 0, -width * (zn - 2))
+    spos2 = (-width * xn, 0, -width * (zn - 1))
+
+    ###
+    #   Bounding box
+    ## 
+    sk = geompy.Sketcher3D()
+
+    for p in point:
+        sk.addPointsAbsolute(*p)
+    
+    inletface = geompy.MakeFaceWires([sk.wire()], False)
+    vecflow = geompy.GetNormal(inletface)
+    hexagonPrism = geompy.MakePrismVecH(inletface, vecflow, L * sqrt(3))
+
+    inletface = geompy.MakeScaleTransform(inletface, oo, scale)
+    hexagonPrism = geompy.MakeScaleTransform(hexagonPrism, oo, scale)
+
+    faces = geompy.ExtractShapes(hexagonPrism, geompy.ShapeType["FACE"], False)
+    symetryface = []
+
+    for face in faces:
+        norm = geompy.GetNormal(face)
+        angle = round(geompy.GetAngle(norm, vecflow), 0)
+
+        if (angle == 0 or angle == 180) and not face == inletface:
+            outletface = face
+        
+        else:
+            symetryface.append(face)
+    
+    ###
+    #   Grains
+    ##
+    ox = geompy.MakeVectorDXDYDZ(1, 0, 0)
+    oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
+    oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
+    xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
+    xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius)
+    lattice1 = geompy.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn)
+    lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn - 1)
+
+    grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius)
+    lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
+    lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn)
+    
+    grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
+    grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
+    grains = geompy.MakeFuseList(grains, False, False)
+
+    grains = geompy.MakeScaleTransform(grains, oo, scale)
+
+    if fillet:
+        grains = geompy.MakeFilletAll(grains, filletradius * scale)
+
+    ###
+    #   Groups
+    ##
+    shape = geompy.MakeCutList(hexagonPrism, [grains])
+    shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "faceCenteredCubic")
+
+    sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
+    geompy.UnionIDs(sall,
+        geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
+
+    inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
+    inletshape = geompy.MakeCutList(inletface, [grains])
+    inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
+    geompy.UnionList(inlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
+
+    outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
+    outletshape = geompy.MakeCutList(outletface, [grains])
+    outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
+    geompy.UnionList(outlet, geompy.SubShapeAll(
+        geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
+    
+    symetry = []
+    for (n, face) in enumerate(symetryface):
+        name = "symetry" + str(n)
+        symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
+        symetryshape = geompy.MakeCutList(face, [grains])
+        symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
+        geompy.UnionList(symetry[n], geompy.SubShapeAll(
+            geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
+
+    groups = []
+    groups.append(inlet)
+    groups.append(outlet)
+    groups.extend(symetry)
+    wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
+
+    return shape, groups
+

samples/simple.py

@@ -7,7 +7,7 @@ geompy = geomBuilder.New()
 
 from math import pi, sqrt
 
-def simpleCubic(theta = 0.01, fillet = False):
+def simpleCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]):
     
     ###
     #   Parameters
@@ -37,17 +37,33 @@ def simpleCubic(theta = 0.01, fillet = False):
     ###
     #   Bounding box
     ##
-    sk = geompy.Sketcher3D()
+    if direction == [1, 0, 0]:
-    sk.addPointsAbsolute(xl, 0, 0)
+        sk = geompy.Sketcher3D()
-    sk.addPointsAbsolute(0, yw, 0)
+        sk.addPointsAbsolute(xl, 0, 0)
-    sk.addPointsAbsolute(0, yw, zh)
+        sk.addPointsAbsolute(0, yw, 0)
-    sk.addPointsAbsolute(xl, 0, zh)
+        sk.addPointsAbsolute(0, yw, zh)
-    sk.addPointsAbsolute(xl, 0, 0)
+        sk.addPointsAbsolute(xl, 0, zh)
-    
+        sk.addPointsAbsolute(xl, 0, 0)
-    inletface = geompy.MakeFaceWires([sk.wire()], True)
-    vecflow = geompy.GetNormal(inletface)
-    cubic = geompy.MakePrismVecH(inletface, vecflow, width)
 
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, width)
+
+    elif direction == [0, 0, 1]:
+        sk = geompy.Sketcher3D()
+        sk.addPointsAbsolute(0, yw, 0)
+        sk.addPointsAbsolute(xl, 0, 0)
+        sk.addPointsAbsolute(2 * xl, yw, 0)
+        sk.addPointsAbsolute(xl, 2 * yw, 0)
+        sk.addPointsAbsolute(0, yw, 0)
+
+        inletface = geompy.MakeFaceWires([sk.wire()], True)
+        vecflow = geompy.GetNormal(inletface)
+        cubic = geompy.MakePrismVecH(inletface, vecflow, height)
+
+    else:
+        raise Exception("The direction is not implemented")
+    
     inletface = geompy.MakeScaleTransform(inletface, oo, scale)
     cubic = geompy.MakeScaleTransform(cubic, oo, scale)
 
@@ -119,11 +135,12 @@ def simpleCubic(theta = 0.01, fillet = False):
     groups.append(outlet)
     groups.extend(symetry)
     wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
 
-    return shape
+    return shape, groups
     
 
-def simpleHexagonalPrism(theta = 0.01, fillet = False):
+def simpleHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]):
     
     ###
     #   Parameters
@@ -239,6 +256,7 @@ def simpleHexagonalPrism(theta = 0.01, fillet = False):
     groups.append(outlet)
     groups.extend(symetry)
     wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
+    groups.append(wall)
+
+    return shape, groups
 
-    return shape
-    

src/mesh_utils.py

@@ -3,12 +3,38 @@ from salome.smesh import smeshBuilder
 smesh = smeshBuilder.New()
 
 import logging
+from collections import namedtuple
+
+Parameters = namedtuple("Parameters", [
+    "minSize",
+    "maxSize",
+    "growthRate",
+    "nbSegPerEdge",
+    "nbSegPerRadius",
+    "chordalErrorEnabled",
+    "chordalError",
+    "secondOrder",
+    "optimize",
+    "quadAllowed",
+    "useSurfaceCurvature",
+    "fuseEdges",
+    "checkChartBoundary"
+])
+
+ViscousLayers = namedtuple("ViscousLayers", [
+    "thickness",
+    "numberOfLayers",
+    "stretchFactor",
+    "isFacesToIgnore",
+    "facesToIgnore",
+    "extrusionMethod"
+])
 
 def getSmesh():
     return smesh
 
 
-def meshCreate(gobj, boundary, fineness, viscousLayers=None):
+def meshCreate(shape, groups, fineness, parameters, viscousLayers = None):
     """
     Creates a mesh from a geometry.
 
@@ -21,19 +47,13 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None):
             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.
 
     """
+    ###
+    #   Netgen
+    ##
     Fineness = {
         0: "Very coarse",
         1: "Coarse",
@@ -43,31 +63,34 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None):
         5: "Custom"
     }[fineness]
 
-    mesh = smesh.Mesh(gobj)
+    # Mesh
+    mesh = smesh.Mesh(shape)
     netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
 
+    # Parameters
     param = netgen.Parameters()
-    param.SetMinSize( 0.001 )
+    param.SetMinSize(parameters.minSize)
-    param.SetMaxSize( 0.1 )
+    param.SetMaxSize(parameters.maxSize)
-    
-    param.SetSecondOrder( 0 )
-    param.SetOptimize( 1 )
-    param.SetQuadAllowed( 0 )
-    param.SetChordalError( -1 )
-    param.SetChordalErrorEnabled( 0 )
-    param.SetUseSurfaceCurvature( 1 )
-    param.SetFuseEdges( 1 )
-    param.SetCheckChartBoundary( 0 )
-    
     param.SetFineness(fineness)
     
-    # TODO: add customization
     if fineness == 5:
-        param.SetGrowthRate( 0.1 )
+        param.SetGrowthRate(parameters.growthRate)
-        param.SetNbSegPerEdge( 5 )
+        param.SetNbSegPerEdge(parameters.nbSegPerEdge)
-        param.SetNbSegPerRadius( 10 )
+        param.SetNbSegPerRadius(parameters.nbSegPerRadius)
     
     
+    param.SetChordalErrorEnabled(parameters.chordalErrorEnabled)
+    param.SetChordalError(parameters.chordalError)
+
+    param.SetSecondOrder(parameters.secondOrder)
+    param.SetOptimize(parameters.optimize)
+    param.SetQuadAllowed(parameters.quadAllowed)
+    
+    param.SetUseSurfaceCurvature(parameters.useSurfaceCurvature)
+    param.SetFuseEdges(parameters.fuseEdges)
+    param.SetCheckChartBoundary(parameters.checkChartBoundary)
+    
+        
     logging.info("""meshCreate:
     fineness:\t{}
     min size:\t{}
@@ -86,29 +109,41 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None):
         True if param.GetSecondOrder() else False, 
         True if param.GetOptimize() else False))
 
+    
+    ###
+    #   Groups
+    ##
+    for group in groups:
+        mesh.GroupOnGeom(group, "{}_".format(group.GetName()), SMESH.FACE)
 
-
+    ###
+    #   Viscous layers
+    ##
     if not viscousLayers is None:
+        vlayer = netgen.ViscousLayers(
+            viscousLayers.thickness,
+            viscousLayers.numberOfLayers,
+            viscousLayers.stretchFactor,
+            viscousLayers.facesToIgnore,
+            viscousLayers.isFacesToIgnore, 
+            viscousLayers.extrusionMethod
+        )
+
         logging.info("""meshCreate:
     viscous layers: 
         thickness:\t{}
         number:\t{}
         stretch factor:\t{}""".format(
-            viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"]))
+            viscousLayers.thickness, 
+            viscousLayers.numberOfLayers, 
+            viscousLayers.stretchFactor))
 
-        vlayer = netgen.ViscousLayers(
+        
-            viscousLayers["thickness"],
-            viscousLayers["number"],
-            viscousLayers["stretch"],
-            [boundary["inlet"], boundary["outlet"]],
-            1, smeshBuilder.NODE_OFFSET)
 
     else:
         logging.info("""meshCreate: 
     viscous layers: disabled""")
 
-    for name, b in boundary.items():
-        mesh.GroupOnGeom(b, "{}_".format(name), SMESH.FACE)
 
     return mesh
 
@@ -175,7 +210,7 @@ def meshExport(mobj, path):
         mobj.ExportUNV(path)
 
         logging.info("""meshExport:
-        format:\t{}""".format("unv"))
+    format:\t{}""".format("unv"))
 
     except:
         logging.error("""meshExport: Cannot export.""")

temp/samples.old.2/__init__.py

@@ -0,0 +1,100 @@
+from collections import namedtuple
+import os, sys
+import logging
+from pyquaternion import Quaternion
+import math
+
+ROOT = "/home/nafaryus/projects/anisotrope-cube"
+sys.path.append(ROOT)
+
+LOG = os.path.join(ROOT, "logs")
+
+import salome
+
+from simpleCubic import simpleCubic
+from faceCenteredCubic import faceCenteredCubic
+from bodyCenteredCubic import bodyCenteredCubic
+
+from src import geometry_utils
+from src import mesh_utils
+
+def genMesh(stype, theta, fillet, flowdirection, saveto):
+    _G = globals()
+
+    structure = _G.get(stype)
+
+    if structure:
+        salome.salome_init()
+
+        grains, geometry1, geometry2 = structure(theta, fillet)
+        geometry = geometry1
+        
+        if flowdirection == [1, 1, 1]:
+            geometry = geometry2
+            norm = [-1, 1, 0]
+            bcount = 6
+
+            # initial angle
+            angle = math.pi / 6
+            v1 = Quaternion(axis = norm, angle = math.pi / 2).rotate(flowdirection)
+            normvec = Quaternion(axis = flowdirection, angle = angle).rotate(v1)
+            direction = [1, 1, 1]
+
+        if flowdirection == [1, 0, 0]:
+            normvec = [0, 0, 1]
+            bcount = 4
+            direction = [1, 1, 0]
+
+        if flowdirection == [0, 0, 1]:
+            normvec = [1, 1, 0]
+            bcount = 4
+            direction = [0, 0, 1]
+        
+        #
+        #geometry = geometry_utils.geompy.RemoveExtraEdges(geometry, False) 
+
+        #
+        #boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains)
+
+        fineness = 1
+        viscousLayers = {
+            "thickness": 0.0001,
+            "number": 3,
+            "stretch": 1.2
+        }
+        mesh = mesh_utils.meshCreate(geometry, boundary, fineness, viscousLayers)
+        mesh_utils.meshCompute(mesh)
+
+        mesh_utils.meshExport(mesh, saveto)
+
+        salome.salome_close()
+
+    else:
+        raise Exception("Unknown sample function")
+
+
+if __name__ == "__main__":
+    
+    logging.basicConfig(
+        level=logging.INFO, 
+        format="%(levelname)s: %(message)s",
+        handlers = [
+            logging.StreamHandler(),
+            logging.FileHandler("{}/cubic.log".format(LOG))
+        ])
+    
+    stype = str(sys.argv[1])
+    theta = float(sys.argv[2])
+    fillet = True if int(sys.argv[3]) == 1 else False
+    flowdirection = [int(coord) for coord in sys.argv[4]]
+    saveto = str(sys.argv[5])
+
+    logging.info("""genMesh: 
+    structure type:\t{}
+    coefficient:\t{}
+    fillet:\t{}
+    flow direction:\t{}
+    export path:\t{}""".format(stype, theta, fillet, flowdirection, saveto))
+
+    genMesh(stype, theta, fillet, flowdirection, saveto)
+