Restruction; New main model; Fixes;

src/anisotropeCubic.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import math
+from . import geometry_utils
+
+import GEOM
+geompy = geometry_utils.getGeom()
+
+class StructuredGrains:
+    def __init__(self, radius, stackAngle, theta, layers):
+        self.pos = [0, 0, 0]
+        self.angle = [0, 0, 0]
+        self.radius = radius
+        self.theta = theta
+        self.layers = layers
+        
+        # Parameters and dependencies
+        R = self.radius / (1 - self.theta)
+        
+        C1 = 0.8 #fillet[0]
+        C2 = 0.4 #fillet[1]
+        self.theta1 = 0.01
+        self.theta2 = 0.28
+        
+        Cf = C1 + (C2 - C1) / (self.theta2 - self.theta1) * (self.theta - self.theta1)
+        R_fillet = Cf * (self.radius * math.sqrt(2) - R)
+        
+        ###
+        stackang = [
+            0.5 * math.pi - stackAngle[0], 
+            0.5 * math.pi - stackAngle[1], 
+            0.5 * math.pi - stackAngle[2]
+        ]
+        
+        xvec = geompy.MakeVector(
+            geompy.MakeVertex(0, 0, 0),
+            geompy.MakeVertex(1, 0, 0))
+        yvec = geometry_utils.rotate(xvec, [0.5 * math.pi, 0, 0])
+        zvec = geometry_utils.rotate(xvec, [0, 0.5 * math.pi, 0])
+        
+        grain = geompy.MakeSpherePntR(geompy.MakeVertex(pos[0], pos[1], pos[2]), R)
+        
+        xstack = geompy.MakeMultiTranslation1D(grain, xvec, 2 * self.radius, self.layers[0])
+        ystack = geompy.MakeMultiTranslation1D(xgrain, yvec, 2 * self.radius, self.layers[1])
+        zstack = geompy.MakeMultiTranslation1D(ygrain, zvec, 2 * self.radius, self.layers[2])
+        
+        # Correct position to zero
+        stack = geompy.MakeTranslation(zstack, -2 * self.radius, 0, 0)
+        
+        self.geometry = geompy.ExtractShapes(stack, geompy.ShapeType["SOLID"], True)
+        self.geometry = geompy.MakeFuseList(self.geometry, False, False)
+                
+        if not R_fillet == 0:
+            self.geometry = geompy.MakeFilletAll(self.geometry, R_fillet)
+
+
+class AnisotropeCubic:
+    def __init__(self, scale, grains, style):
+        self.pos = [0, 0, 0]
+        self.angle = [0, 0, 0]
+        self.scale = scale
+        self.grains = grains
+
+        # Bounding box
+        if style == 0:
+            # Square
+            profile = (
+                geompy.Sketcher3D()
+                .addPointAbsolute(0, 0, 0)
+                .addPointAbsolute(0, 0, self.scale[2])
+                .addPointAbsolute(0, self.scale[1], self.scale[2])
+                .addPointAbsolute(0, self.scale[1], 0)
+                .addPointAbsolute(0, 0, 0)
+            )
+            
+            face = geompy.MakeFaceWires([profile.wire()], 1)
+
+        elif style == 1:
+            # Rombus
+            profile = (
+                geompy.Sketcher3D()
+                .addPointAbsolute(self.scale[0], 0.5 * self.scale[1], 0)
+                .addPointAbsolute(0.5 * self.scale[0], 0, 0.5 * self.scale[2])
+                .addPointAbsolute(0, 0.5 * self.scale[1], self.scale[2])
+                .addPointAbsolute(0.5 * self.scale[0], self.scale[1], 0.5 * self.scale[2])
+                .addPointAbsolute(self.scale[0], 0.5 * self.scale[1], 0)
+            )
+
+            face = geompy.MakeFaceWires([profile.wire()], 1)
+            face = geompy.MakeTranslation(face, 
+                0.5 * self.scale[1], 0, 0)
+
+        self.boundingbox = geompy.MakePrismVecH(face,
+            geompy.MakeVectorDXDYDZ(1, 0, 0), 
+            self.scale[0])
+        
+        # Geometry
+        self.geometry = geompy.MakeCutList(box, [self.grains], True)
+
+
+
+

src/geometry_utils.py

@@ -0,0 +1,122 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import GEOM
+from salome.geom import geomBuilder
+geompy = geomBuilder.New()
+
+import math
+import logging
+
+def getGeom():
+    return geompy
+
+def rotate(gobj, ang):
+    x = geompy.MakeVectorDXDYDZ(1, 0, 0),
+    y = geompy.MakeVectorDXDYDZ(0, 1, 0),
+    z = geompy.MakeVectorDXDYDZ(0, 0, 1)
+
+    # yaw
+    rotated = geompy.MakeRotation(gobj, z, ang[2])
+    # pitch
+    rotated = geompy.MakeRotation(rotated, y, ang[1])
+    # roll
+    rotated = geompy.MakeRotation(rotated, x, ang[0])
+
+    return rotated
+
+def createGroup(gobj, planelist, grains, name):
+    gr = geompy.CreateGroup(gobj, geompy.ShapeType["FACE"], name)
+
+    grcomp = geompy.MakeCompound(planelist)
+    grcut = geompy.MakeCutList(grcomp, [grains], True)
+
+    gip = geompy.GetInPlace(gobj, grcut, True)
+    faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
+    geompy.UnionList(gr, faces)
+
+    return gr
+
+def boundaryCreate(gobj, dvec, grains):
+
+    xvec = geompy.MakeVector(
+        geompy.MakeVertex(0, 0, 0),
+        geompy.MakeVertex(dvec[0], dvec[1], dvec[2]))
+    xvec = rotate(dvec, self.angle)
+
+    yvec = rotate(xvec, [0.5 * math.pi, 0, 0])
+    zvec = rotate(xvec, [0, 0.5 * math.pi, 0])
+
+    logging.info("boundaryCreate: dvec = {}".format(dvec))
+
+    planes = geompy.ExtractShapes(gobj, geompy.ShapeType["FACE"], True)
+    inletplanes = []
+    outletplanes = []
+    uplanes = []
+
+    fwplanes = []
+    bwplanes = []
+    lplanes = []
+    rplanes = []
+
+    for plane in planes:
+        nvec = geompy.GetNormal(plane)
+        xang = geompy.GetAngle(nvec, xvec)
+        yang = geompy.GetAngle(nvec, yvec)
+        zang = geompy.GetAngle(nvec, zvec)
+
+        if xang == 0:
+            inletplanes.append(plane)
+
+        elif xang == 180:
+            outletplanes.append(plane)
+
+        elif yang == 0:
+            fwplanes.append(plane)
+
+        elif yang == 180:
+            bwplanes.append(plane)
+
+        elif zang == 0:
+            lplanes.append(plane)
+
+        elif zang == 180:
+            rplanes.append(plane)
+
+    logging.info(
+        "boundaryCreate: inletplanes = {}, outletplanes = {}, hplanes = {}".format(
+            len(inletplane), len(outletplane), len(hplanes)))
+
+    logging.info(
+        "boundaryCreate: fwplanes = {}, bwplanes = {}, lplanes = {}, rplanes = {}".format(
+            len(fwplanes), len(bwplanes), len(lplanes), len(rplanes)))
+
+    # Main groups
+    inlet = createGroup(gobj, inletplane, grains, "inlet")
+    outlet = createGroup(gobj, grains, outletplane, "outlet")
+
+    symetryPlaneFW = createGroup(gobj, fwplanes, grains, "symetryPlaneFW")
+    symetryPlaneBW = createGroup(gobj, bwplanes, grains, "symetryPlaneBW")
+    symetryPlaneL = createGroup(gobj, lplanes, grains, "symetryPlaneL")
+    symetryPlaneR = createGroup(gobj, rplanes, grains, "symetryPlaneR")
+
+    # wall
+    allgroup = geompy.CreateGroup(gobj, geompy.ShapeType["FACE"])
+    faces = geompy.SubShapeAllIDs(gobj, geompy.ShapeType["FACE"])
+    geompy.UnionIDs(allgroup, faces)
+    wall = geompy.CutListOfGroups([allgroup],
+        [inlet, outlet, symetryPlaneFW, symetryPlaneBW, symetryPlaneL, symetryPlaneR], "wall")
+
+    boundary = {
+        "inlet": inlet,
+        "outlet": outlet,
+        "symetryPlaneFW": symetryPlaneFW,
+        "symetryPlaneBW": symetryPlaneBW,
+        "symetryPlaneL": symetryPlaneL,
+        "symetryPlaneR": symetryPlaneR,
+        "wall": wall
+    }
+
+    return boundary
+
+

src/mesh_utils.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import SMESH
+from salome.smesh import smeshBuilder
+smesh = smeshBuilder.New()
+
+import logging
+
+def getSmesh():
+    return smesh
+
+def meshCreate(gobj, boundary, 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.
+
+    """
+    Fineness = {
+        0: "Very coarse",
+        1: "Coarse",
+        2: "Moderate",
+        3: "Fine",
+        4: "Very fine"
+    }[fineness]
+
+    logging.info("meshCreate: mesh fineness - {}".format(Fineness))
+
+    mesh = smesh.Mesh(gobj)
+    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(fineness)
+    #param.SetGrowthRate( 0.1 )
+    #param.SetNbSegPerEdge( 5 )
+    #param.SetNbSegPerRadius( 10 )
+    param.SetQuadAllowed( 0 )
+
+    if not viscousLayers is None:
+        logging.info("meshCreate: viscous layers params - thickness = {}, number = {}, stretch factor = {}".format(
+            viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"]))
+
+        vlayer = netgen.ViscousLayers(
+            viscousLayers["thickness"],
+            viscousLayers["number"],
+            viscousLayers["stretch"],
+            [boundary["inlet"], boundary["outlet"]],
+            1, smeshBuilder.NODE_OFFSET)
+
+    else:
+        logging.info("meshCreate: viscous layers are disabled")
+
+    for name, b in boundary.items():
+        mesh.GroupOnGeom(b, "{}_".format(name), SMESH.FACE)
+
+    return mesh
+
+def meshCompute(mobj):
+    """Compute the mesh."""
+    status = mobj.Compute()
+
+    if status:
+        logging.info("Mesh succesfully computed.")
+
+    else:
+        logging.warning("Mesh is not computed.")
+
+def meshExport(mobj, 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(mobj.name))
+
+    try:
+        mobj.ExportUNV(exportpath)
+
+    except:
+        logging.error("Cannot export mesh to '{}'".format(exportpath))
+
+

src/salome_utils.py

@@ -0,0 +1,10 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import salome
+
+def hasDesktop() -> bool:
+    return salome.sg.hasDesktop()
+
+def execute():
+    pass

src/samples.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from . import geometry_utils
+import GEOM
+
+from . import mesh_utils
+import SMESH
+
+from . import anisotropeCubic
+
+def simpleCubic(theta, flowdirection):
+    radius = 1
+    stackAngle = [0.5 * math.pi, 0.5 * math.pi, 0.5 * math.pi]
+    theta = theta if theta else 0.1
+    layers = [3, 3, 3]
+    grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
+    
+    scale = [1, 1, 1]
+    flowdirection = flowdirection if flowdirection else [1, 0, 0]
+    style = 0
+    cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
+    boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
+    
+    fineness = 3
+    mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
+    mesh_utils.meshCompute(mesh)
+
+
+def bodyCenteredCubic(theta, flowdirection):
+    radius = 1
+    stackAngle = [0.5 * math.pi, 0.25 * math.pi, 0.25 * math.pi]
+    theta = theta if theta else 0.1
+    layers = [3, 3, 3]
+    grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
+    
+    scale = [1, 1, 1]
+    flowdirection = flowdirection if flowdirection else [1, 0, 0]
+    style = 0
+    cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
+    boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
+    
+    fineness = 3
+    mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
+    mesh_utils.meshCompute(mesh)
+
+
+def faceCenteredCubic(theta, flowdirection):
+    radius = 1
+    stackAngle = [0.5 * math.pi, 0.5 * math.pi, 0.5 * math.pi]
+    theta = theta if theta else 0.1
+    layers = [3, 3, 3]
+    grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
+    
+    scale = [1, 1, 1]
+    flowdirection = flowdirection if flowdirection else [1, 0, 0]
+    style = 0
+    cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
+    boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
+    
+    fineness = 3
+    mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
+    mesh_utils.meshCompute(mesh)
+
+
+def genMesh():
+    pass