Mod: shaping is stable now

Mod: shaping improved documentation
Remove: not necessary modules or moved

anisotropy/meshing/__init__.py

@@ -1,17 +1,16 @@
 # -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
 
 from . import utils
 from . import conversion
 from . import metrics
 
-from .mesh import Mesh
+from .mesh import Mesh, MeshingParameters
 
 
 __all__ = [
     "utils",
     "conversion",
     "metrics",
-    "Mesh"
+    "Mesh",
+    "MeshingParameters"
 ]

anisotropy/meshing/mesh.py

@@ -196,7 +196,7 @@ class Mesh:
             mesh.write(path)
 
     @property
-    def volumes(self) -> list[ndarray]: #   delete?
+    def volumes(self) -> list[ndarray]: 
         """Volumes.
 
         :return:
@@ -220,7 +220,7 @@ class Mesh:
         return np.sum([ metrics.volume(cell) for cell in self.volumes ]) 
 
     @property
-    def faces(self) -> list[ndarray]: #    delete?
+    def faces(self) -> list[ndarray]: 
         """Boundary faces.
 
         :return:

anisotropy/shaping/__init__.py

@@ -1,9 +1,21 @@
 # -*- coding: utf-8 -*-
-# This file is part of anisotropy.
+"""
-# License: GNU GPL version 3, see the file "LICENSE" for details.
+Shaping is a library for using OCC shapes, provides more convient
+functionality with power NumPy and Python OOP and contains interesing
+primitives.
+"""
 
-from .shape import ShapeError, Shape
+from . import utils
+
+from .shape import Shape
 from .periodic import Periodic
-from .simple import Simple
+
-from .faceCentered import FaceCentered
+from . import primitives
-from .bodyCentered import BodyCentered
+
+
+__all__ = [
+    "utils",
+    "primitives",
+    "Shape",
+    "Periodic"
+]

anisotropy/shaping/bodyCentered.py

@@ -1,173 +0,0 @@
-# -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
-
-from netgen.occ import *
-import numpy
-from numpy import pi, sqrt, arccos
-
-from .occExtended import *
-from . import Periodic
-from . import ShapeError
-
-class BodyCentered(Periodic):
-    def __init__(
-        self,
-        direction: list = None,
-        **kwargs
-    ):
-        Periodic.__init__(
-            self,
-            alpha = kwargs.get("alpha", 0.01),
-            r0 = kwargs.get("r0", 1),
-            filletsEnabled = kwargs.get("filletsEnabled", True),
-            gamma = pi - 2 * arccos(sqrt(2 / 3))
-        )
-
-        #   Parameters
-        self.direction = direction
-        self.alphaMin = 0.01
-        self.alphaMax = 0.18
-
-        #   Objects
-        self.lattice = None
-        self.cell = None
-        self.shape = None
-
-    
-    @property
-    def L(self):
-        return self.r0 * 4 / sqrt(3)
-
-
-    def build(self):
-        #   
-        zero = Pnt(0, 0, 0)
-
-        #   Lattice
-        spheres = numpy.array([], dtype = object)
-
-        for zn in range(3):
-            z = zn * self.L
-            z2 = z - 0.5 * self.L
-
-            for yn in range(3):
-                y = yn * self.L
-                y2 = y - 0.5 * self.L
-
-                for xn in range(3):
-                    x = xn * self.L
-                    x2 = x - 0.5 * self.L
-
-                    spheres = numpy.append(spheres, Sphere(Pnt(x, y, z), self.radius))
-                    # shifted
-                    spheres = numpy.append(spheres, Sphere(Pnt(x2, y2, z2), self.radius))
-
-        lattice = spheres.sum()
-        lattice = lattice.Scale(zero, self.maximize)
-
-        if self.filletsEnabled:
-            lattice = lattice.MakeFillet(lattice.edges, self.fillets * self.maximize)
-
-        self.lattice = lattice.Scale(zero, self.minimize)
-        
-        #   Inlet face
-        if (self.direction == numpy.array([1., 0., 0.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(2)
-            height = self.L
-
-            xl = sqrt(diag ** 2 + diag ** 2) * 0.5
-            yw = xl
-            zh = height
-
-            vertices = numpy.array([
-                (xl, 0, 0),
-                (0, yw, 0),
-                (0, yw, zh),
-                (xl, 0, zh)
-            ])
-            extr = diag
-
-        elif (self.direction == numpy.array([0., 0., 1.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(2)
-            height = self.L
-
-            xl = sqrt(diag ** 2 + diag ** 2) * 0.5
-            yw = xl
-            zh = height
-
-            vertices = numpy.array([
-                (0, yw, 0),
-                (xl, 0, 0),
-                (2 * xl, yw, 0),
-                (xl, 2 * yw, 0)
-            ])
-            extr = height
-
-        elif (self.direction == numpy.array([1., 1., 1.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(2)
-            height = diag / 3
-
-            xl = -self.L / 4
-            yw = -self.L / 4
-            zh = -self.L / 4
-
-            vertices = numpy.array([
-                (self.L / 3 + xl, self.L / 3 + yw, 4 * self.L / 3 + zh),
-                (self.L + xl, 0 + yw, self.L + zh),
-                (4 * self.L / 3 + xl, self.L / 3 + yw, self.L / 3 + zh),
-                (self.L + xl, self.L + yw, 0 + zh),
-                (self.L / 3 + xl, 4 * self.L / 3 + yw, self.L / 3 + zh),
-                (0 + xl, self.L + yw, self.L + zh)
-            ])
-            extr = self.L * sqrt(3) 
-
-        else:
-            raise Exception(f"Direction { self.direction } is not implemented")
-
-        #   Cell
-        circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ])
-        inletface = Face(circuit)
-        inletface.name = "inlet"
-
-        vecFlow = self.normal(inletface)
-        # ISSUE: don't use face.Extrude(length), only face.Extrude(length, vector) 
-        self.cell = inletface.Extrude(extr, Vec(*vecFlow))
-        
-        #   Boundaries
-        symetryId = 0
-
-        for face in self.cell.faces:
-            fNorm = self.normal(face)
-            fAngle = self.angle(vecFlow, fNorm)
-            
-            if fAngle == 0 or fAngle == numpy.pi:
-                if (face.center.pos() == inletface.center.pos()).prod():
-                    face.name = "inlet"
-
-                else:
-                    face.name = "outlet"
-
-            else:
-                face.name = f"symetry{ symetryId }"
-                symetryId += 1
-
-        #   Main shape
-        self.shape = self.cell - self.lattice
-
-        if not len(self.shape.solids) == 1:
-            raise ShapeError("Expected single solid shape")
-
-        else:
-            self.shape = self.shape.solids[0]
-
-        #   Boundaries (walls)
-        for face in self.shape.faces:
-            if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
-                face.name = "wall"

anisotropy/shaping/faceCentered.py

@@ -1,182 +0,0 @@
-# -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
-
-from netgen.occ import *
-import numpy
-from numpy import pi, sqrt
-
-from .occExtended import *
-from . import Periodic
-from . import ShapeError
-
-class FaceCentered(Periodic):
-    def __init__(
-        self,
-        direction: list = None,
-        **kwargs
-    ):
-        Periodic.__init__(
-            self,
-            alpha = kwargs.get("alpha", 0.01),
-            r0 = kwargs.get("r0", 1),
-            filletsEnabled = kwargs.get("filletsEnabled", True),
-            gamma = 2 / 3 * pi
-        )
-
-        #   Parameters
-        self.direction = direction
-        self.alphaMin = 0.01
-        self.alphaMax = 0.13
-
-        #   Objects
-        self.lattice = None
-        self.cell = None
-        self.shape = None
-
-
-    @property
-    def L(self):
-        return self.r0 * 4 / sqrt(2)
-        
-
-    def build(self):
-        #   
-        zero = Pnt(0, 0, 0)
-
-        #   Lattice
-        spheres = numpy.array([], dtype = object)
-        x0 = 0#-0.5 * self.L * (3 - 1)
-        x20 = 0#-0.5 * self.L * 3
-        z0 = -0.5 * self.L * (3 - 2)
-        z20 = -0.5 * self.L * (3 - 1)
-
-        for zn in range(3):
-            z = z0 + zn * self.L 
-            z2 = z20 + zn * self.L 
-
-            for yn in range(3):
-                y = yn * 2 * self.r0
-                y2 = yn * 2 * self.r0 + self.r0
-
-                for xn in range(3):
-                    x = x0 + xn * 2 * self.r0
-                    x2 = x20 + xn * 2 * self.r0 + self.r0
-
-                    # TODO: fix rotations (arcs intersection -> incorrect boolean operations
-                    spheres = numpy.append(spheres, Sphere(Pnt(x, y, z), self.radius).Rotate(Axis(Pnt(x, y, z), X), 45).Rotate(Axis(Pnt(x, y, z), Z), 45))
-                    # shifted
-                    spheres = numpy.append(spheres, Sphere(Pnt(x2, y2, z2), self.radius).Rotate(Axis(Pnt(x2, y2, z2), X), 45).Rotate(Axis(Pnt(x2, y2, z2), Z), 45))
-
-
-        lattice = spheres.sum()
-        lattice = lattice.Move(Vec(-self.r0 * 2, -self.r0 * 2, 0)).Rotate(Axis(zero, Z), 45)
-        lattice = lattice.Scale(zero, self.maximize)
-
-        if self.filletsEnabled:
-            lattice = lattice.MakeFillet(lattice.edges, self.fillets * self.maximize)
-
-        self.lattice = lattice.Scale(zero, self.minimize)
-
-        #   Inlet face
-        if (self.direction == numpy.array([1., 0., 0.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(3)
-            height = diag / 3
-
-            xl = sqrt(length ** 2 + length ** 2) * 0.5
-            yw = xl
-            zh = width
-
-            vertices = numpy.array([
-                (0, 0, -zh),
-                (-xl, yw, -zh),
-                (-xl, yw, zh),
-                (0, 0, zh)
-            ])
-            extr = length
-
-        elif (self.direction == numpy.array([0., 0., 1.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(3)
-            height = diag / 3
-
-            xl = sqrt(length ** 2 + length ** 2) * 0.5
-            yw = xl
-            zh = width
-
-            vertices = numpy.array([
-                (0, 0, -zh),
-                (xl, yw, -zh),
-                (0, 2 * yw, -zh),
-                (-xl, yw, -zh)
-            ])
-            extr = 2 * width
-
-        elif (self.direction == numpy.array([1., 1., 1.])).prod():
-            length = 2 * self.r0
-            width = self.L / 2
-            diag = self.L * sqrt(3)
-            height = diag / 3
-
-            xl = -(3 - 2) * self.L / 3
-            yw = -(3 - 2) * self.L / 3
-            zh = -(3 - 2) * self.L / 3
-
-            vertices = numpy.array([
-                (-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh),
-                (0 + xl, -width + yw, 0 + zh),
-                (width / 3 + xl, -2 * width / 3 + yw, -2 * width / 3 + zh),
-                (0 + xl, 0 + yw, -width + zh),
-                (-2 * width / 3 + xl, width / 3 + yw, -2 * width / 3 + zh),
-                (-width + xl, 0 + yw, 0 + zh)
-            ])
-            extr = sqrt(3) * self.L
-            
-        else:
-            raise Exception(f"Direction { self.direction } is not implemented")
-        
-        #   Cell
-        circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ])
-        inletface = Face(circuit)
-        inletface.name = "inlet"
-
-        vecFlow = self.normal(inletface)
-        # ISSUE: don't use face.Extrude(length), only face.Extrude(length, vector) 
-        self.cell = inletface.Extrude(extr, Vec(*vecFlow))
-        
-        #   Boundaries
-        symetryId = 0
-
-        for face in self.cell.faces:
-            fNorm = self.normal(face)
-            fAngle = self.angle(vecFlow, fNorm)
-            
-            if fAngle == 0 or fAngle == numpy.pi:
-                if (face.center.pos() == inletface.center.pos()).prod():
-                    face.name = "inlet"
-
-                else:
-                    face.name = "outlet"
-
-            else:
-                face.name = f"symetry{ symetryId }"
-                symetryId += 1
-
-        #   Main shape
-        self.shape = self.cell - self.lattice
-
-        if not len(self.shape.solids) == 1:
-            raise ShapeError("Expected single solid shape")
-
-        else:
-            self.shape = self.shape.solids[0]
-
-        #   Boundaries (walls)
-        for face in self.shape.faces:
-            if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
-                face.name = "wall"
-
-

anisotropy/shaping/occExtended.py

@@ -1,41 +0,0 @@
-# -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
-
-from functools import wraps
-from netgen import occ
-import numpy
-
-
-def add_method(cls):
-    """Add method to existing class. Use it as decorator.
-    """
-    def decorator(func):
-        @wraps(func) 
-        def wrapper(self, *args, **kwargs): 
-            return func(self, *args, **kwargs)
-
-        setattr(cls, func.__name__, wrapper)
-        
-        return func
-    return decorator
-
-
-@add_method(occ.gp_Pnt)
-def pos(self) -> numpy.array:
-    return numpy.array([self.x, self.y, self.z])
-
-
-# ISSUE: netgen.occ.Face.Extrude: the opposite face has the same name and normal vector as an initial face.
-#def reconstruct(shape):
-#    """Reconstruct shape with new objects. 
-#    """
-#    faces = []
-#
-#    for face in shape.faces:
-#        faceNew = occ.Face(face.wires[0])
-#        faceNew.name = face.name
-#        faces.append(faceNew)
-#
-#    return occ.Solid(faces)
-

anisotropy/shaping/periodic.py

@@ -1,9 +1,6 @@
 # -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from netgen.occ import *
+import numpy as np
-from numpy import pi, sqrt, cos, arccos, fix 
 
 from . import Shape
 
@@ -13,85 +10,60 @@ class Periodic(Shape):
         self,
         alpha: float = None,
         r0: float = 1,
-        #L: float = None,
+        L: float = None,
-        #radius: float = None,
         filletsEnabled: bool = True,
-        #fillets: float = None,
         gamma = None,
         **kwargs
     ):
-        """Constructor method.
+        """A Periodic object contains periodic structure.
         
         :param alpha:
             Spheres overlap parameter.
-        
         :param r0:
             Initial spheres radius.
-        
+        :param L:
+            Side length of periodic cell. Depends on r0.
+        :param gamma:
+            Angle between lines that connect centers of spheres.
         :param filletsEnabled:
-            Enable fillets beetween spheres.
+            If True, calculate fillets beetween spheres.
         """
         Shape.__init__(self)
         
         #   Parameters
         self.alpha = alpha
         self.r0 = r0
+        self.L = L
 
-        self.theta = 0.5 * pi
+        #   for lattice
+        # self.theta = 0.5 * pi
 
-        self.gamma = gamma or pi - 2 * (0.5 * self.theta)
+        self.gamma = gamma or np.pi - 2 * (0.5 * 0.5 * np.pi)
         self.filletsEnabled = filletsEnabled
         self.filletsScale = 0.8
         
-        # scale parameters for precision boolean operations
+        #   scale parameters for precision boolean operations
         self.maximize = 1e+2
         self.minimize = 1e-2
 
-    
     @property
-    def L(self):
+    def radius(self) -> float:
-        """(Override) Parameter depending on the ``r0``.
-        """
-        pass
-
-    @property
-    def radius(self):
         """Spheres radius
+
+        :return:
+            Radius.
         """
         return self.r0 / (1 - self.alpha)
 
     @property
     def fillets(self):
-        analytical = self.r0 * sqrt(2) / sqrt(1 - cos(self.gamma)) - self.radius
+        """Fillets radius beetween spheres.
-        # ISSUE: MakeFilletAll : Fillet can't be computed on the given shape with the given radius.
+
-        # Temporary solution: degrade the precision (minRound <= analytical).
+        :return:
+            Radius.
+        """
+        analytical = self.r0 * np.sqrt(2) / np.sqrt(1 - np.cos(self.gamma)) - self.radius
         rTol = 3
-        minRound = fix(10 ** rTol * analytical) * 10 ** -rTol
+        minRound = np.fix(10 ** rTol * analytical) * 10 ** -rTol
 
         return minRound * self.filletsScale
-
-    def lattice(self, 
-        theta = 0.5 * pi
-    ):
-        zero = Pnt(0, 0, 0)
-        maximize = 1e+2
-        minimize = 1e-2
-
-        #   Lattice
-        spheres = numpy.array([], dtype = object)
-
-        for zn in range(3):
-            z = zn * self.L
-            for yn in range(3):
-                y = yn * self.L
-                for xn in range(3):
-                    x = xn * self.L
-                    spheres = numpy.append(spheres, Sphere(Pnt(x, y, z), self.radius))
-
-        lattice = spheres.sum()
-        lattice = lattice.Scale(zero, maximize)
-
-        if self.filletsEnabled:
-            lattice = lattice.MakeFillet(lattice.edges, self.fillets * maximize)
-
-        self.lattice = lattice.Scale(zero, minimize)

anisotropy/shaping/primitives.py

@@ -0,0 +1,499 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import annotations
+from numpy import ndarray
+
+import numpy as np
+import netgen.occ as ng_occ
+
+from . import Periodic
+from . import utils
+
+
+def simple(alpha: float, direction: list | ndarray, **kwargs) -> Periodic:
+    """Simple periodic structure.
+
+    :param alpha:
+        Spheres overlap parameter.
+    :param direction:
+        Flow direction vector. This parameter affects the geometry type
+        and boundary (faces) names.
+    :return:
+        Periodic object.
+    """
+    #   base object
+    pc = Periodic(
+        alpha = alpha,
+        gamma = np.pi - 2 * (0.5 * 0.5 * np.pi)
+    )
+    
+    #   additional parameters
+    pc.__dict__.update(kwargs)
+    pc.L = 2 * pc.r0
+    pc.direction = np.asarray(direction)
+
+    #   additional attributes
+    pc.cell: ng_occ.Solid = None
+    pc.lattice: ng_occ.Solid = None
+
+    #   constants
+    zero = ng_occ.Pnt(0, 0, 0)
+
+    #   Lattice
+    spheres = [] 
+
+    for zn in range(3):
+        z = zn * pc.L
+
+        for yn in range(3):
+            y = yn * pc.L
+
+            for xn in range(3):
+                x = xn * pc.L
+
+                spheres.append(ng_occ.Sphere(ng_occ.Pnt(x, y, z), pc.radius))
+
+    lattice = np.sum(spheres)
+    lattice = lattice.Scale(zero, pc.maximize)
+
+    if pc.filletsEnabled:
+        lattice = lattice.MakeFillet(lattice.edges, pc.fillets * pc.maximize)
+
+    pc.lattice = lattice.Scale(zero, pc.minimize)
+    
+    #   Inlet face
+    if np.all(pc.direction == [1., 0., 0.]):
+        length = pc.L * np.sqrt(2)
+        width = pc.L * np.sqrt(2)
+        height = pc.L
+
+        xl = np.sqrt(length ** 2 * 0.5)
+        yw = xl
+        zh = height
+
+        vertices = np.array([
+            (xl, 0, 0),
+            (0, yw, 0),
+            (0, yw, zh),
+            (xl, 0, zh)
+        ])
+        extr = width
+
+    elif np.all(pc.direction == [0., 0., 1.]):
+        length = pc.L * np.sqrt(2)
+        width = pc.L * np.sqrt(2)
+        height = pc.L
+
+        xl = np.sqrt(length ** 2 * 0.5)
+        yw = xl
+        zh = height
+
+        vertices = np.array([
+            (0, yw, 0),
+            (xl, 0, 0),
+            (2 * xl, yw, 0),
+            (xl, 2 * yw, 0)
+        ])
+        extr = height
+
+    elif np.all(pc.direction == [1., 1., 1.]):
+        length = pc.L * np.sqrt(2)
+        width = pc.L * np.sqrt(2)
+        height = pc.L
+
+        xl = -pc.L - pc.L / 6 
+        yw = -pc.L - pc.L / 6
+        zh = -pc.L / 6
+
+        vertices = np.array([
+            (pc.L + xl, pc.L + yw, pc.L + zh),
+            (5 * pc.L / 3 + xl, 2 * pc.L / 3 + yw, 2 * pc.L / 3 + zh),
+            (2 * pc.L + xl, pc.L + yw, 0 + zh),
+            (5 * pc.L / 3 + xl, 5 * pc.L / 3 + yw, -pc.L / 3 + zh),
+            (pc.L + xl, 2 * pc.L + yw, 0 + zh),
+            (2 * pc.L / 3 + xl, 5 * pc.L / 3 + yw, 2 * pc.L / 3 + zh)
+        ])
+        extr = pc.L * np.sqrt(3) 
+
+    else:
+        raise Exception(f"Direction { pc.direction } is not implemented")
+    
+    #   Cell
+    circuit = ng_occ.Wire([ 
+        ng_occ.Segment(ng_occ.Pnt(*v1), ng_occ.Pnt(*v2))
+        for v1, v2 in zip(vertices, np.roll(vertices, -1, axis = 0)) 
+    ])
+    inletface = ng_occ.Face(circuit)
+    inletface.name = "inlet"
+
+    vecFlow = utils.normal(inletface)
+    pc.cell = inletface.Extrude(extr, ng_occ.Vec(*vecFlow))
+    
+    #   Boundary faces
+    symetryId = 0
+
+    for face in pc.cell.faces:
+        fNorm = utils.normal(face)
+        fAngle = utils.angle(vecFlow, fNorm)
+        
+        if fAngle == 0 or fAngle == np.pi:
+            if np.all(utils.pos(face.center) == utils.pos(inletface.center)):
+                face.name = "inlet"
+                
+            else:
+                face.name = "outlet"
+
+        else:
+            face.name = f"symetry{ symetryId }"
+            symetryId += 1
+
+    #   Main shape
+    pc.shape = pc.cell - pc.lattice
+
+    assert len(pc.shape.solids) == 1, "Expected single solid shape"
+
+    pc.shape = pc.shape.solids[0]
+
+    #   Boundary faces (walls)
+    for face in pc.shape.faces:
+        if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
+            face.name = "wall"
+
+    return pc
+
+
+def body_centered(alpha: float, direction: list | ndarray, **kwargs) -> Periodic:
+    """Body centered periodic structure.
+
+    :param alpha:
+        Spheres overlap parameter.
+    :param direction:
+        Flow direction vector. This parameter affects the geometry type
+        and boundary (faces) names.
+    :return:
+        Periodic object.
+    """
+    #   base object 
+    pc = Periodic(
+        alpha = alpha,
+        gamma = np.pi - 2 * np.arccos(np.sqrt(2 / 3))
+    )
+
+    #   additional parameters
+    pc.__dict__.update(kwargs)
+    pc.L = 4 / np.sqrt(3) * pc.r0
+    pc.direction = np.asarray(direction)
+
+    #   additional attributes
+    pc.cell: ng_occ.Solid = None
+    pc.lattice: ng_occ.Solid = None
+
+    #   constants   
+    zero = ng_occ.Pnt(0, 0, 0)
+
+    #   Lattice
+    spheres = [] 
+
+    for zn in range(3):
+        z = zn * pc.L
+        z2 = z - 0.5 * pc.L
+
+        for yn in range(3):
+            y = yn * pc.L
+            y2 = y - 0.5 * pc.L
+
+            for xn in range(3):
+                x = xn * pc.L
+                x2 = x - 0.5 * pc.L
+
+                spheres.append(ng_occ.Sphere(ng_occ.Pnt(x, y, z), pc.radius))
+                # shifted
+                spheres.append(ng_occ.Sphere(ng_occ.Pnt(x2, y2, z2), pc.radius))
+
+    lattice = np.sum(spheres)
+    lattice = lattice.Scale(zero, pc.maximize)
+
+    if pc.filletsEnabled:
+        lattice = lattice.MakeFillet(lattice.edges, pc.fillets * pc.maximize)
+
+    pc.lattice = lattice.Scale(zero, pc.minimize)
+    
+    #   Inlet face
+    if np.all(pc.direction == [1., 0., 0.]):
+        # length = 2 * pc.r0
+        # width = pc.L / 2
+        diag = pc.L * np.sqrt(2)
+        height = pc.L
+
+        xl = np.sqrt(diag ** 2 + diag ** 2) * 0.5
+        yw = xl
+        zh = height
+
+        vertices = np.array([
+            (xl, 0, 0),
+            (0, yw, 0),
+            (0, yw, zh),
+            (xl, 0, zh)
+        ])
+        extr = diag
+
+    elif np.all(pc.direction == [0., 0., 1.]):
+        # length = 2 * pc.r0
+        # width = pc.L / 2
+        diag = pc.L * np.sqrt(2)
+        height = pc.L
+
+        xl = np.sqrt(diag ** 2 + diag ** 2) * 0.5
+        yw = xl
+        zh = height
+
+        vertices = np.array([
+            (0, yw, 0),
+            (xl, 0, 0),
+            (2 * xl, yw, 0),
+            (xl, 2 * yw, 0)
+        ])
+        extr = height
+
+    elif np.all(pc.direction == [1., 1., 1.]):
+        # length = 2 * pc.r0
+        # width = pc.L / 2
+        diag = pc.L * np.sqrt(2)
+        height = diag / 3
+
+        xl = -pc.L / 4
+        yw = -pc.L / 4
+        zh = -pc.L / 4
+
+        vertices = np.array([
+            (pc.L / 3 + xl, pc.L / 3 + yw, 4 * pc.L / 3 + zh),
+            (pc.L + xl, 0 + yw, pc.L + zh),
+            (4 * pc.L / 3 + xl, pc.L / 3 + yw, pc.L / 3 + zh),
+            (pc.L + xl, pc.L + yw, 0 + zh),
+            (pc.L / 3 + xl, 4 * pc.L / 3 + yw, pc.L / 3 + zh),
+            (0 + xl, pc.L + yw, pc.L + zh)
+        ])
+        extr = pc.L * np.sqrt(3) 
+
+    else:
+        raise Exception(f"Direction { pc.direction } is not implemented")
+
+    #   Cell
+    circuit = ng_occ.Wire([ 
+        ng_occ.Segment(ng_occ.Pnt(*v1), ng_occ.Pnt(*v2)) 
+        for v1, v2 in zip(vertices, np.roll(vertices, -1, axis = 0)) 
+    ])
+    inletface = ng_occ.Face(circuit)
+    inletface.name = "inlet"
+
+    vecFlow = utils.normal(inletface)
+    pc.cell = inletface.Extrude(extr, ng_occ.Vec(*vecFlow))
+    
+    #   Boundary faces
+    symetryId = 0
+
+    for face in pc.cell.faces:
+        fNorm = utils.normal(face)
+        fAngle = utils.angle(vecFlow, fNorm)
+        
+        if fAngle == 0 or fAngle == np.pi:
+            if np.all(utils.pos(face.center) == utils.pos(inletface.center)):
+                face.name = "inlet"
+
+            else:
+                face.name = "outlet"
+
+        else:
+            face.name = f"symetry{ symetryId }"
+            symetryId += 1
+
+    #   Main shape
+    pc.shape = pc.cell - pc.lattice
+
+    assert len(pc.shape.solids) == 1, "Expected single solid shape"
+
+    pc.shape = pc.shape.solids[0]
+
+    #   Boundary faces (walls)
+    for face in pc.shape.faces:
+        if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
+            face.name = "wall"
+
+    return pc
+
+
+def face_centered(alpha: float, direction: list | ndarray, **kwargs) -> Periodic:
+    """Face centered periodic structure.
+
+    :param alpha:
+        Spheres overlap parameter.
+    :param direction:
+        Flow direction vector. This parameter affects the geometry type
+        and boundary (faces) names.
+    :return:
+        Periodic object.
+    """
+    #   base class
+    pc = Periodic(
+        alpha = alpha,
+        gamma = 2 / 3 * np.pi 
+    )
+
+    #   additional parameters
+    pc.__dict__.update(kwargs)
+    pc.L = 4 / np.sqrt(2) * pc.r0
+    pc.direction = np.asarray(direction)
+
+    #   additional attributes
+    pc.cell: ng_occ.Solid = None
+    pc.lattice: ng_occ.Solid = None
+
+    #   constants   
+    zero = ng_occ.Pnt(0, 0, 0)
+
+    #   Lattice
+    spheres = [] 
+    x0 = 0  
+    x20 = 0  
+    z0 = -0.5 * pc.L * (3 - 2)
+    z20 = -0.5 * pc.L * (3 - 1)
+
+    for zn in range(3):
+        z = z0 + zn * pc.L 
+        z2 = z20 + zn * pc.L 
+
+        for yn in range(3):
+            y = yn * 2 * pc.r0
+            y2 = yn * 2 * pc.r0 + pc.r0
+
+            for xn in range(3):
+                x = x0 + xn * 2 * pc.r0
+                x2 = x20 + xn * 2 * pc.r0 + pc.r0
+
+                # TODO: fix rotations (arcs intersection -> incorrect boolean operations
+                spheres.append(
+                    ng_occ.Sphere(ng_occ.Pnt(x, y, z), pc.radius)
+                    .Rotate(ng_occ.Axis(ng_occ.Pnt(x, y, z), ng_occ.X), 45)
+                    .Rotate(ng_occ.Axis(ng_occ.Pnt(x, y, z), ng_occ.Z), 45)
+                )
+                # shifted
+                spheres.append(
+                    ng_occ.Sphere(ng_occ.Pnt(x2, y2, z2), pc.radius)
+                    .Rotate(ng_occ.Axis(ng_occ.Pnt(x2, y2, z2), ng_occ.X), 45)
+                    .Rotate(ng_occ.Axis(ng_occ.Pnt(x2, y2, z2), ng_occ.Z), 45)
+                )
+
+    lattice = np.sum(spheres) 
+    lattice = (
+        lattice.Move(ng_occ.Vec(-pc.r0 * 2, -pc.r0 * 2, 0))
+        .Rotate(ng_occ.Axis(zero, ng_occ.Z), 45)
+    )
+    lattice = lattice.Scale(zero, pc.maximize)
+
+    if pc.filletsEnabled:
+        lattice = lattice.MakeFillet(lattice.edges, pc.fillets * pc.maximize)
+
+    pc.lattice = lattice.Scale(zero, pc.minimize)
+
+    #   Inlet face
+    if np.all(pc.direction == [1., 0., 0.]):
+        length = 2 * pc.r0
+        width = pc.L / 2
+        # diag = pc.L * np.sqrt(3)
+        # height = diag / 3
+
+        xl = np.sqrt(length ** 2 + length ** 2) * 0.5
+        yw = xl
+        zh = width
+
+        vertices = np.array([
+            (0, 0, -zh),
+            (-xl, yw, -zh),
+            (-xl, yw, zh),
+            (0, 0, zh)
+        ])
+        extr = length
+
+    elif np.all(pc.direction == [0., 0., 1.]):
+        length = 2 * pc.r0
+        width = pc.L / 2
+        # diag = pc.L * np.sqrt(3)
+        # height = diag / 3
+
+        xl = np.sqrt(length ** 2 + length ** 2) * 0.5
+        yw = xl
+        zh = width
+
+        vertices = np.array([
+            (0, 0, -zh),
+            (xl, yw, -zh),
+            (0, 2 * yw, -zh),
+            (-xl, yw, -zh)
+        ])
+        extr = 2 * width
+
+    elif np.all(pc.direction == [1., 1., 1.]):
+        length = 2 * pc.r0
+        width = pc.L / 2
+        # diag = pc.L * np.sqrt(3)
+        # height = diag / 3
+
+        xl = -(3 - 2) * pc.L / 3
+        yw = -(3 - 2) * pc.L / 3
+        zh = -(3 - 2) * pc.L / 3
+
+        vertices = np.array([
+            (-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh),
+            (0 + xl, -width + yw, 0 + zh),
+            (width / 3 + xl, -2 * width / 3 + yw, -2 * width / 3 + zh),
+            (0 + xl, 0 + yw, -width + zh),
+            (-2 * width / 3 + xl, width / 3 + yw, -2 * width / 3 + zh),
+            (-width + xl, 0 + yw, 0 + zh)
+        ])
+        extr = np.sqrt(3) * pc.L
+        
+    else:
+        raise Exception(f"Direction { pc.direction } is not implemented")
+    
+    #   Cell
+    circuit = ng_occ.Wire([ 
+        ng_occ.Segment(ng_occ.Pnt(*v1), ng_occ.Pnt(*v2)) 
+        for v1, v2 in zip(vertices, np.roll(vertices, -1, axis = 0)) 
+    ])
+    inletface = ng_occ.Face(circuit)
+    inletface.name = "inlet"
+
+    vecFlow = utils.normal(inletface)
+    pc.cell = inletface.Extrude(extr, ng_occ.Vec(*vecFlow))
+    
+    #   Boundary faces
+    symetryId = 0
+
+    for face in pc.cell.faces:
+        fNorm = utils.normal(face)
+        fAngle = utils.angle(vecFlow, fNorm)
+        
+        if fAngle == 0 or fAngle == np.pi:
+            if np.all(utils.pos(face.center) == utils.pos(inletface.center)):
+                face.name = "inlet"
+
+            else:
+                face.name = "outlet"
+
+        else:
+            face.name = f"symetry{ symetryId }"
+            symetryId += 1
+
+    #   Main shape
+    pc.shape = pc.cell - pc.lattice
+
+    assert len(pc.shape.solids) == 1, "Expected single solid shape"
+
+    pc.shape = pc.shape.solids[0]
+
+    #   Boundary faces (walls)
+    for face in pc.shape.faces:
+        if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
+            face.name = "wall"
+
+    return pc

anisotropy/shaping/shape.py

@@ -1,80 +1,95 @@
 # -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from netgen.occ import *
+from __future__ import annotations
-import numpy
+from numpy import ndarray
-from numpy import linalg
-import os
 from os import PathLike
-from pathlib import Path
 
+import numpy as np
+import netgen.occ as ng_occ
+import pathlib
 
-class ShapeError(Exception):
+from . import utils
-    pass
 
 
 class Shape(object):
     def __init__(self):
+        """A Shape object contains OCC shape.
+        """
         self.groups = {}
         self.shape = None
 
+    @property
+    def geometry(self) -> ng_occ.OCCGeometry:
+        """Shape as OCCGeometry object.
+        """
+        return ng_occ.OCCGeometry(self.shape)
+    
+    @property
+    def type(self) -> ng_occ.TopAbs_ShapeEnum:
+        """Type of the shape. (shortcut)
+        """
+        return self.shape.type
+    
+    @property
+    def volume(self) -> float:
+        """Volume of the shape. (shortcut)
+        """
+        return self.shape.volume
+
+    @property
+    def center(self) -> ndarray:
+        """Center of the shape.
+        """
+        return np.array(utils.pos(self.shape.center))
+
     def write(self, filename: PathLike):
-        """Export a shape.
+        """Export a shape to the file.
-        
         Supported formats: step.
         
         :param filename:
-            Name of the file to store the given shape in.
+            Path of the file.
-        
-        :return:
-            Output, error messages and returncode
         """
-        out, err, returncode = "", "", 0
+        path = pathlib.Path(filename).resolve()
-        path = Path(filename).resolve()
         ext = path.suffix[1: ]
         
-        try:
+        if ext == "step":
-            if ext == "step":
+            self.shape.WriteStep(str(path))
-                self.shape.WriteStep(path)
-            
-            else:
-                raise NotImplementedError(f"{ ext } is not supported")
-            
-        except NotImplementedError as e:
-            err = e
-            returncode = 1
         
-        except Exception as e:
+        else:
-            err = e
+            raise NotImplementedError(f"Shape format '{ ext }' is not supported")
-            returncode = 1
+            
-        
-        return out, err, returncode
-
     def read(self, filename: PathLike):
-        path = Path(filename).resolve()
+        """Import a shape from the file.
+        Supported formats: step, iges, brep.
+
+        :param filename:
+            Path of the file.
+        """        
+        path = pathlib.Path(filename).resolve()
         ext = path.suffix[1: ]
 
         if ext in ["step", "iges", "brep"]:
-            self.shape = OCCGeometry(path).shape
+            self.shape = ng_occ.OCCGeometry(str(path)).shape
 
         else:
             raise NotImplementedError(f"Shape format '{ext}' is not supported")
         
         return self
 
-    def patches(self, group: bool = False, shiftIndex: bool = False, prefix: str = None):
+    def patches(
+        self, 
+        group: bool = False, 
+        shiftIndex: bool = False, 
+        prefix: str = None
+    ) -> list | dict:
         """Get patches indices with their names.
 
         :param group:
             Group indices together with the same patches names.
-
         :param shiftIndex:
             Start numerating with one instead of zero.
-
         :param prefix:
             Add string prefix to the index.
-
         :return:
             List if group = False else dictionary.
         """
@@ -105,24 +120,3 @@ class Shape(object):
                 patches_.append((item, face.name))
 
         return patches_ 
-
-    def normal(self, face: FACE) -> numpy.array:
-        """
-        :return:
-            Normal vector to face.
-        """
-        _, u, v = face.surf.D1(0, 0)
-
-        return numpy.cross([u.x, u.y, u.z], [v.x, v.y, v.z]) 
-
-
-    def angle(self, vec1: numpy.array, vec2: numpy.array) -> float:
-        """
-        :return:
-            Angle between two vectors in radians.
-        """
-        inner = numpy.inner(vec1, vec2)
-        norms = linalg.norm(vec1) * linalg.norm(vec2)
-        cos = inner / norms
-
-        return numpy.arccos(numpy.clip(cos, -1.0, 1.0))

anisotropy/shaping/simple.py

@@ -1,170 +0,0 @@
-# -*- coding: utf-8 -*-
-# This file is part of anisotropy.
-# License: GNU GPL version 3, see the file "LICENSE" for details.
-
-from netgen.occ import *
-import numpy
-from numpy import pi, sqrt
-
-from .occExtended import *
-from . import Periodic
-from . import ShapeError
-
-class Simple(Periodic):
-    def __init__(
-        self,
-        direction: list = None,
-        **kwargs
-    ):
-        Periodic.__init__(
-            self,
-            alpha = kwargs.get("alpha", 0.01),
-            r0 = kwargs.get("r0", 1),
-            filletsEnabled = kwargs.get("filletsEnabled", True)
-        )
-
-        #   Parameters
-        self.direction = direction
-        self.alphaMin = 0.01
-        self.alphaMax = 0.28
-
-        #   Objects
-        self.lattice = None
-        self.cell = None
-        self.shape = None
-
-
-    @property
-    def L(self):
-        return 2 * self.r0
-
-
-    def build(self):
-        #   
-        zero = Pnt(0, 0, 0)
-
-        #   Lattice
-        spheres = numpy.array([], dtype = object)
-
-        for zn in range(3):
-            z = zn * self.L
-
-            for yn in range(3):
-                y = yn * self.L
-
-                for xn in range(3):
-                    x = xn * self.L
-
-                    spheres = numpy.append(spheres, Sphere(Pnt(x, y, z), self.radius))
-
-        lattice = spheres.sum()
-        lattice = lattice.Scale(zero, self.maximize)
-
-        if self.filletsEnabled:
-            lattice = lattice.MakeFillet(lattice.edges, self.fillets * self.maximize)
-
-        self.lattice = lattice.Scale(zero, self.minimize)
-        
-        #   Inlet face
-        if (self.direction == numpy.array([1., 0., 0.])).prod():
-            length = self.L * numpy.sqrt(2)
-            width = self.L * numpy.sqrt(2)
-            height = self.L
-
-            xl = numpy.sqrt(length ** 2 * 0.5)
-            yw = xl
-            zh = height
-
-            vertices = numpy.array([
-                (xl, 0, 0),
-                (0, yw, 0),
-                (0, yw, zh),
-                (xl, 0, zh)
-            ])
-            extr = width
-
-        elif (self.direction == numpy.array([0., 0., 1.])).prod():
-            length = self.L * numpy.sqrt(2)
-            width = self.L * numpy.sqrt(2)
-            height = self.L
-
-            xl = numpy.sqrt(length ** 2 * 0.5)
-            yw = xl
-            zh = height
-
-            vertices = numpy.array([
-                (0, yw, 0),
-                (xl, 0, 0),
-                (2 * xl, yw, 0),
-                (xl, 2 * yw, 0)
-            ])
-            extr = height
-
-        elif (self.direction == numpy.array([1., 1., 1.])).prod():
-            length = self.L * numpy.sqrt(2)
-            width = self.L * numpy.sqrt(2)
-            height = self.L
-
-            xl = -self.L - self.L / 6 
-            yw = -self.L - self.L / 6
-            zh = -self.L / 6
-
-            vertices = numpy.array([
-                (self.L + xl, self.L + yw, self.L + zh),
-                (5 * self.L / 3 + xl, 2 * self.L / 3 + yw, 2 * self.L / 3 + zh),
-                (2 * self.L + xl, self.L + yw, 0 + zh),
-                (5 * self.L / 3 + xl, 5 * self.L / 3 + yw, -self.L / 3 + zh),
-                (self.L + xl, 2 * self.L + yw, 0 + zh),
-                (2 * self.L / 3 + xl, 5 * self.L / 3 + yw, 2 * self.L / 3 + zh)
-            ])
-            extr = self.L * sqrt(3) 
-
-        else:
-            raise Exception(f"Direction { self.direction } is not implemented")
-        
-        #
-        
-
-        #   Cell
-        circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ])
-        inletface = Face(circuit)
-        inletface.name = "inlet"
-
-        vecFlow = self.normal(inletface)
-        # ISSUE: don't use face.Extrude(length), only face.Extrude(length, vector) 
-        self.cell = inletface.Extrude(extr, Vec(*vecFlow))
-        
-        #   Boundaries
-        symetryId = 0
-
-        for face in self.cell.faces:
-            fNorm = self.normal(face)
-            fAngle = self.angle(vecFlow, fNorm)
-            
-            if fAngle == 0 or fAngle == numpy.pi:
-                if (face.center.pos() == inletface.center.pos()).prod():
-                    face.name = "inlet"
-                    
-                else:
-                    face.name = "outlet"
-
-            else:
-                face.name = f"symetry{ symetryId }"
-                symetryId += 1
-
-        
-        #   Main shape
-        self.shape = self.cell - self.lattice
-
-        if not len(self.shape.solids) == 1:
-            raise ShapeError("Expected single solid shape")
-
-        else:
-            self.shape = self.shape.solids[0]
-
-        #   Boundaries (walls)
-        for face in self.shape.faces:
-            if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
-                face.name = "wall"
-
-

anisotropy/shaping/utils.py

@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+
+from numpy import ndarray
+
+import numpy as np
+import netgen.occ as ng_occ
+
+
+def pos(point: ng_occ.gp_Pnt) -> ndarray:
+    """Extract coordinates from point.
+
+    :param point:
+        OCC point object
+    :return:
+        Array of coordinates.
+    """
+    return np.array([ point.x, point.y, point.z ])
+
+
+def normal(face: ng_occ.Face) -> ndarray: 
+    """Calculate normal vector from face.
+
+    :param face:
+        OCC face object.
+    :return:
+        Normal vector from face.
+    """
+    _, u, v = face.surf.D1(0, 0)
+
+    return np.cross([u.x, u.y, u.z], [v.x, v.y, v.z]) 
+
+
+def angle(vec1: ndarray, vec2: ndarray) -> float:
+    """Angle between two vectors.
+
+    :param vec1:
+        Array of points that represents first vector.
+    :param vec2:
+        Array of points that represents second vector.
+    :return:
+        Angle between two vectors in radians.
+    """
+    inner = np.inner(vec1, vec2)
+    norms = np.linalg.norm(vec1) * np.linalg.norm(vec2)
+    cos = inner / norms
+
+    return np.arccos(np.clip(cos, -1.0, 1.0))