Fix: fillet radius on all structures

Move: StructureGeometry class to geometry.py

anisotropy/salomepl/geometry.py

@@ -4,22 +4,212 @@
 
 import logging
 
-logger = logging.getLogger("anisotropy")
+GEOM_IMPORTED = False
 
 try:
     import GEOM
     from salome.geom import geomBuilder
 
-except ImportError:
+    GEOM_IMPORTED = True
-    logger.debug("Trying to get SALOME geometry modules outside SALOME environment. Modules won't be imported.")
 
-if globals().get("geomBuilder"):
+except:
-    geompy = geomBuilder.New()
+    pass
-
-else:
-    geompy = None
-
-def getGeom():
-    return geompy
 
 
+class StructureGeometry(object):
+    def __init__(
+        self,
+        direction: list = None,
+        theta: float = None,
+        r0: float = 1,
+        #L: float = None,
+        #radius: float = None,
+        filletsEnabled: bool = False,
+        #fillets: float = None,
+        **kwargs
+    ):
+        """Constructor method.
+        
+        :param direction:
+            Flow vector that characterizes geometry.
+        
+        :param theta:
+            Spheres overlap parameter.
+        
+        :param r0:
+            Initial spheres radius.
+        
+        :param filletsEnabled:
+            Enable fillets beetween spheres.
+        """
+        #   Geometry parameters
+        self.direction = direction
+        self.theta = theta
+        self.r0 = r0
+        self.filletsEnabled = filletsEnabled
+        #self.fillets = fillets
+        
+        #   General attributes
+        self.shape = None
+        self.groups = []
+        self.shapeCell = None
+        self.shapeLattice = None
+
+        #   Geometry module
+        if not GEOM_IMPORTED:
+            raise ImportError("Cannot find the salome modules.")
+
+        else:
+            self.geo = geomBuilder.New() 
+    
+    @property
+    def name(self):
+        """(Override) Shape name.
+        """
+        pass
+    
+    @property
+    def L(self):
+        """(Override) Parameter depending on the ``r0``.
+        """
+        pass
+    
+    @property
+    def radius(self):
+        """Spheres radius
+        """
+        return self.r0 / (1 - self.theta)
+    
+    @property
+    def volumeCell(self):
+        """General volume of the cell.
+        """
+        return self.geo.BasicProperties(self.shapeCell, theTolerance = 1e-06)[2]
+    
+    @property
+    def volume(self):
+        """Volume of the structure.
+        """
+        return self.geo.BasicProperties(self.shape, theTolerance = 1e-06)[2]
+    
+    @property
+    def porosity(self):
+        """Porosity of the structure.
+        """
+        return self.volume / self.volumeCell
+    
+    def build(self):
+        """(Override) Construct shape and physical groups.
+        """
+        pass
+    
+    def isValid(self) -> (bool, str):
+        """Check a topology of the given shape.
+        
+        :return:
+            True, if the shape "seems to be valid" else False and description.
+        """
+        return self.geo.CheckShape(self.shape, theIsCheckGeom = True, theReturnStatus = 1)
+    
+    def heal(self):
+        """Try to heal the shape.
+        """
+        self.shape = self.geo.RemoveExtraEdges(self.shape, doUnionFaces = False)
+    
+    def createGroupAll(self, mainShape):
+        """Create group from all the shape faces.
+        
+        :param mainShape:
+            Input shape.
+        
+        :return:
+            Created group.
+        """
+        group = self.geo.CreateGroup(mainShape, self.geo.ShapeType["FACE"])
+        
+        self.geo.UnionIDs(
+            group,
+            self.geo.SubShapeAllIDs(mainShape, self.geo.ShapeType["FACE"])
+        )
+        
+        return group
+        
+    def createGroup(self, mainShape, subShape, name: str, cutShapes: list = None, scaleTransform: float = None):
+        """Create group from the sub shape.
+        
+        :param mainShape:
+            Input shape.
+        
+        :param subShape:
+            Input sub shape.
+        
+        :param name:
+            Name of the new group.
+        
+        :param cutShapes:
+            List of shapes for cut from the sub shape.
+            
+        :param scaleTransform:
+            Value of the scale transform regarding to the zero point.
+        
+        :return:
+            Created group.
+        """
+        group = self.geo.CreateGroup(mainShape, self.geo.ShapeType["FACE"], theName = name)
+        
+        if cutShapes:
+            subShape = self.geo.MakeCutList(subShape, cutShapes)
+        
+        if scaleTransform:
+            subShape = self.geo.MakeScaleTransform(subShape, self.geo.MakeVertex(0, 0, 0), scaleTransform)
+        
+        self.geo.UnionList(
+            group, 
+            self.geo.SubShapeAll(
+                self.geo.GetInPlace(mainShape, subShape, True), 
+                self.geo.ShapeType["FACE"]
+            )
+        )
+        
+        return group
+    
+    def export(
+        filename: str,
+        deflection: float = 0.001
+    ):
+        """Export a shape.
+        
+        Supported formats: step, vtk.
+        
+        :param filename:
+            Name of the file to store the given shape in.
+        
+        :param deflection:
+            vtk: Deflection of the given shape.
+        
+        :return:
+            Output, error messages and returncode
+        """
+        out, err, returncode = "", "", 0
+        ext = os.path.splitext(filename)[1][1: ]
+        
+        try:
+            if ext == "step":
+                self.geo.ExportSTEP(self.shape, filename) # theUnit = GEOM.LU_METER)
+            
+            elif ext == "vtk":
+                self.geo.ExportVTK(self.shape, filename, theDeflection = deflection)
+            
+            else:
+                raise NotImplementedError(f"{ ext } is not supported")
+            
+        except NotImplementedError as e:
+            err = e
+            returncode = 1
+        
+        except Exception as e:
+            err = e.details.text
+            returncode = 1
+        
+        return out, err, returncode
+ 

anisotropy/samples/bodyCentered.py

@@ -2,8 +2,8 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from anisotropy.samples.structure import StructureGeometry
+from anisotropy.salomepl.geometry import StructureGeometry
-from math import pi, sqrt
+from numpy import pi, sqrt, cos, arccos, fix 
 import logging
 
 class BodyCentered(StructureGeometry):
@@ -27,19 +27,14 @@ class BodyCentered(StructureGeometry):
     
     @property
     def fillets(self):
-        if self.direction == [1.0, 1.0, 1.0]:
+        analytical = self.r0 * (sqrt(2) / sqrt(1 - cos(pi - 2 * arccos(sqrt(2 / 3)))) - 
-            C1, C2 = 0.3, 0.2
+            1 / (1 - self.theta))
-            Cf = C1 + (C2 - C1) / (self.thetaMax - self.thetaMin) * (self.theta - self.thetaMin)
+        # ISSUE: MakeFilletAll : Fillet can't be computed on the given shape with the given radius.
-            delta = 0.02
+        # Temporary solution: degrade the precision (minRound <= analytical).
-            
+        rTol = 3
-            return delta - Cf * (self.radius - self.r0)
+        minRound = fix(10 ** rTol * analytical) * 10 ** -rTol
-        
+
-        else:
+        return minRound
-            C1, C2 = 0.3, 0.2
-            Cf = C1 + (C2 - C1) / (self.thetaMax - self.thetaMin) * (self.theta - self.thetaMin)
-            delta = 0.02
-            
-            return delta - Cf * (self.radius - self.r0)
         
     def build(self):
         ###
@@ -172,7 +167,7 @@ class BodyCentered(StructureGeometry):
             raise Exception(f"Direction { self.direction } is not implemented")
         
         ###
-        #   Grains
+        #   Lattice
         ##
         ox = self.geo.MakeVectorDXDYDZ(1, 0, 0)
         oy = self.geo.MakeVectorDXDYDZ(0, 1, 0)
@@ -184,10 +179,7 @@ class BodyCentered(StructureGeometry):
         lattice1 = self.geo.MakeMultiTranslation2D(grain, ox, self.L, xn, oy, self.L, yn)
         lattice1 = self.geo.MakeMultiTranslation1D(lattice1, oz, self.L, zn)
 
-        #grain = self.geo.MakeSpherePntR(self.geo.MakeVertex(*spos2), radius)
+        lattice2 = self.geo.MakeTranslation(lattice1, -0.5 * self.L, -0.5 * self.L, -0.5 * self.L)
-        #lattice2 = self.geo.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
-        #lattice2 = self.geo.MakeMultiTranslation1D(lattice2, oz, L, zn)
-        lattice2 = self.geo.MakeTranslation(lattice1, 0.5 * self.L, 0.5 * self.L, 0.5 * self.L)
         
         grains = self.geo.ExtractShapes(lattice1, self.geo.ShapeType["SOLID"], True)
         grains += self.geo.ExtractShapes(lattice2, self.geo.ShapeType["SOLID"], True)

anisotropy/samples/faceCentered.py

@@ -2,8 +2,8 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from anisotropy.samples.structure import StructureGeometry
+from anisotropy.salomepl.geometry import StructureGeometry
-from math import pi, sqrt
+from numpy import pi, sqrt, fix
 import logging
 
 class FaceCentered(StructureGeometry):
@@ -27,11 +27,13 @@ class FaceCentered(StructureGeometry):
     
     @property
     def fillets(self):
-        if self.direction == [1.0, 1.0, 1.0]:
+        analytical = self.r0 * (2 * sqrt(3) / 3 - 1 / (1 - self.theta))
-            return self.filletsScale * (2 * sqrt(3) / 3 - 1 / (1 - self.theta)) * self.r0 
+        # ISSUE: MakeFilletAll : Fillet can't be computed on the given shape with the given radius.
-        
+        # Temporary solution: degrade the precision (minRound <= analytical).
-        else:
+        rTol = 3
-            return self.filletsScale * (2 * sqrt(3) / 3 - 1 / (1 - self.theta)) * self.r0 
+        minRound = fix(10 ** rTol * analytical) * 10 ** -rTol
+
+        return minRound
         
     def build(self):
         ###

anisotropy/samples/simple.py

@@ -2,8 +2,9 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from anisotropy.samples.structure import StructureGeometry
+from anisotropy.salomepl.geometry import StructureGeometry
-from math import pi, sqrt
+from numpy import pi, sqrt, fix
+import logging
 
 class Simple(StructureGeometry):
     @property
@@ -26,20 +27,13 @@ class Simple(StructureGeometry):
     
     @property
     def fillets(self):
-        if self.direction == [1.0, 1.0, 1.0]:
+        analytical = self.r0 * (sqrt(2) - 1 / (1 - self.theta))
-            C1, C2 = 0.8, 0.5
+        # ISSUE: MakeFilletAll : Fillet can't be computed on the given shape with the given radius.
-            Cf = C1 + (C2 - C1) / (self.thetaMax - self.thetaMin) * (self.theta - self.thetaMin)
+        # Temporary solution: degrade the precision (minRound <= analytical).
-            delta = 0.2
+        rTol = 3
-            
+        minRound = fix(10 ** rTol * analytical) * 10 ** -rTol
-            return delta - Cf * (self.radius - self.r0)
+
-        
+        return minRound
-        else:
-            C1, C2 = 0.8, 0.5
-            Cf = C1 + (C2 - C1) / (self.thetaMax - self.thetaMin) * (self.theta - self.thetaMin)
-            delta = 0.2
-            
-            return delta - Cf * (self.radius - self.r0)
-            
 
     def build(self):
         ###
@@ -192,9 +186,10 @@ class Simple(StructureGeometry):
         self.shape = self.geo.MakeCutList(poreCell, [grains])
         self.shape = self.geo.MakeScaleTransform(self.shape, oo, 1 / scale, theName = self.name)
         
-        isValid, _ = self.isValid()
+        isValid, msg = self.isValid()
         
         if not isValid:
+            logging.warning(msg)
             self.heal()
             
         ###

anisotropy/samples/structure.py

@@ -1,197 +0,0 @@
-# -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
-
-from anisotropy.salomepl import geometry
-
-class StructureGeometry(object):
-    def __init__(
-        self,
-        direction: list = None,
-        theta: float = None,
-        r0: float = 1,
-        #L: float = None,
-        #radius: float = None,
-        filletsEnabled: bool = False,
-        #fillets: float = None,
-        **kwargs
-    ):
-        """Constructor method.
-        
-        :param direction:
-            Flow vector that characterizes geometry.
-        
-        :param theta:
-            Spheres overlap parameter.
-        
-        :param r0:
-            Initial spheres radius.
-        
-        :param filletsEnabled:
-            Enable fillets beetween spheres.
-        """
-        #    Geometry parameters
-        self.direction = direction
-        self.theta = theta
-        self.r0 = r0
-        self.filletsEnabled = filletsEnabled
-        #self.fillets = fillets
-        
-        #    General attributes
-        self.geo = geometry.getGeom()
-        self.shape = None
-        self.groups = []
-        self.shapeCell = None
-        self.shapeLattice = None
-    
-    @property
-    def name(self):
-        """(Override) Shape name.
-        """
-        pass
-    
-    @property
-    def L(self):
-        """(Override) Parameter depending on the ``r0``.
-        """
-        pass
-    
-    @property
-    def radius(self):
-        """Spheres radius
-        """
-        return self.r0 / (1 - self.theta)
-    
-    @property
-    def volumeCell(self):
-        """General volume of the cell.
-        """
-        return self.geo.BasicProperties(self.shapeCell, theTolerance = 1e-06)[2]
-    
-    @property
-    def volume(self):
-        """Volume of the structure.
-        """
-        return self.geo.BasicProperties(self.shape, theTolerance = 1e-06)[2]
-    
-    @property
-    def porosity(self):
-        """Porosity of the structure.
-        """
-        return self.volume / self.volumeCell
-    
-    def build(self):
-        """(Override) Construct shape and physical groups.
-        """
-        pass
-    
-    def isValid(self) -> (bool, str):
-        """Check a topology of the given shape.
-        
-        :return:
-            True, if the shape "seems to be valid" else False and description.
-        """
-        return self.geo.CheckShape(self.shape, theIsCheckGeom = True, theReturnStatus = 1)
-    
-    def heal(self):
-        """Try to heal the shape.
-        """
-        self.shape = self.geo.RemoveExtraEdges(self.shape, doUnionFaces = False)
-    
-    def createGroupAll(self, mainShape):
-        """Create group from all the shape faces.
-        
-        :param mainShape:
-            Input shape.
-        
-        :return:
-            Created group.
-        """
-        group = self.geo.CreateGroup(mainShape, self.geo.ShapeType["FACE"])
-        
-        self.geo.UnionIDs(
-            group,
-            self.geo.SubShapeAllIDs(mainShape, self.geo.ShapeType["FACE"])
-        )
-        
-        return group
-        
-    def createGroup(self, mainShape, subShape, name: str, cutShapes: list = None, scaleTransform: float = None):
-        """Create group from the sub shape.
-        
-        :param mainShape:
-            Input shape.
-        
-        :param subShape:
-            Input sub shape.
-        
-        :param name:
-            Name of the new group.
-        
-        :param cutShapes:
-            List of shapes for cut from the sub shape.
-            
-        :param scaleTransform:
-            Value of the scale transform regarding to the zero point.
-        
-        :return:
-            Created group.
-        """
-        group = self.geo.CreateGroup(mainShape, self.geo.ShapeType["FACE"], theName = name)
-        
-        if cutShapes:
-            subShape = self.geo.MakeCutList(subShape, cutShapes)
-        
-        if scaleTransform:
-            subShape = self.geo.MakeScaleTransform(subShape, self.geo.MakeVertex(0, 0, 0), scaleTransform)
-        
-        self.geo.UnionList(
-            group, 
-            self.geo.SubShapeAll(
-                self.geo.GetInPlace(mainShape, subShape, True), 
-                self.geo.ShapeType["FACE"]
-            )
-        )
-        
-        return group
-    
-    def export(
-        filename: str,
-        deflection: float = 0.001
-    ):
-        """Export a shape.
-        
-        Supported formats: step, vtk.
-        
-        :param filename:
-            Name of the file to store the given shape in.
-        
-        :param deflection:
-            vtk: Deflection of the given shape.
-        
-        :return:
-            Output, error messages and returncode
-        """
-        out, err, returncode = "", "", 0
-        ext = os.path.splitext(filename)[1][1: ]
-        
-        try:
-            if ext == "step":
-                self.geo.ExportSTEP(self.shape, filename) # theUnit = GEOM.LU_METER)
-            
-            elif ext == "vtk":
-                self.geo.ExportVTK(self.shape, filename, theDeflection = deflection)
-            
-            else:
-                raise NotImplementedError(f"{ ext } is not supported")
-            
-        except NotImplementedError as e:
-            err = e
-            returncode = 1
-        
-        except Exception as e:
-            err = e.details.text
-            returncode = 1
-        
-        return out, err, returncode
-        

playground/database.ipynb

@@ -454,7 +454,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.6"
+   "version": "3.10.0"
   }
  },
  "nbformat": 4,

playground/mesh.py

@@ -20,7 +20,7 @@ class Netgen3D(MeshAlgorithm):
         self.hypo.SetMinSize(value)
 
 class MEFISTO(MeshAlgorithm):
-    
+    pass 
     
 class Mesh(object):
     def __init__(self, geom):
@@ -116,4 +116,4 @@ class Mesh(object):
             returncode = 1
         
         return out, err, returncode
-        
+        

playground/salome-run.py

@@ -0,0 +1,5 @@
+# -*- coding: utf-8 -*-
+
+import sys; path = "/home/nafaryus/projects/anisotropy"; sys.path.extend([path, path + "/env/lib/python3.10/site-packages"])
+from anisotropy.samples import Simple, FaceCentered
+s = Simple([1, 0, 0], 0.28, filletsEnabled = True)

playground/test_smesh.py

@@ -17,3 +17,6 @@ mesh = smesh.Mesh(SHAPE)
 
 algo = mesh.Tetrahedron(algo = ALGO)
 
+###
+#   fc 111 0.12 N3D N2D(0.1, 0.0001, Mod)
+#   3min 321k

playground/training-models.ipynb

@@ -295,7 +295,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.6"
+   "version": "3.10.0"
   }
  },
  "nbformat": 4,