From 00f7279c6dfd7574a5139d39276344b14803f20e Mon Sep 17 00:00:00 2001 From: L-Nafaryus Date: Sat, 22 Jan 2022 22:08:47 +0500 Subject: [PATCH] Mod: shaping is stable now Mod: shaping improved documentation Remove: not necessary modules or moved --- anisotropy/meshing/__init__.py | 7 +- anisotropy/meshing/mesh.py | 4 +- anisotropy/shaping/__init__.py | 24 +- anisotropy/shaping/bodyCentered.py | 173 ---------- anisotropy/shaping/faceCentered.py | 182 ----------- anisotropy/shaping/occExtended.py | 41 --- anisotropy/shaping/periodic.py | 76 ++--- anisotropy/shaping/primitives.py | 499 +++++++++++++++++++++++++++++ anisotropy/shaping/shape.py | 114 ++++--- anisotropy/shaping/simple.py | 170 ---------- anisotropy/shaping/utils.py | 47 +++ 11 files changed, 647 insertions(+), 690 deletions(-) delete mode 100644 anisotropy/shaping/bodyCentered.py delete mode 100644 anisotropy/shaping/faceCentered.py delete mode 100644 anisotropy/shaping/occExtended.py create mode 100644 anisotropy/shaping/primitives.py delete mode 100644 anisotropy/shaping/simple.py create mode 100644 anisotropy/shaping/utils.py diff --git a/anisotropy/meshing/__init__.py b/anisotropy/meshing/__init__.py index 45c93ea..982fad5 100644 --- a/anisotropy/meshing/__init__.py +++ b/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" ] diff --git a/anisotropy/meshing/mesh.py b/anisotropy/meshing/mesh.py index e24c50d..77b3b4e 100644 --- a/anisotropy/meshing/mesh.py +++ b/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: diff --git a/anisotropy/shaping/__init__.py b/anisotropy/shaping/__init__.py index bf99910..68afd17 100644 --- a/anisotropy/shaping/__init__.py +++ b/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 .bodyCentered import BodyCentered + +from . import primitives + + +__all__ = [ + "utils", + "primitives", + "Shape", + "Periodic" +] diff --git a/anisotropy/shaping/bodyCentered.py b/anisotropy/shaping/bodyCentered.py deleted file mode 100644 index a14a056..0000000 --- a/anisotropy/shaping/bodyCentered.py +++ /dev/null @@ -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" diff --git a/anisotropy/shaping/faceCentered.py b/anisotropy/shaping/faceCentered.py deleted file mode 100644 index 7b06776..0000000 --- a/anisotropy/shaping/faceCentered.py +++ /dev/null @@ -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" - - diff --git a/anisotropy/shaping/occExtended.py b/anisotropy/shaping/occExtended.py deleted file mode 100644 index d2cccb6..0000000 --- a/anisotropy/shaping/occExtended.py +++ /dev/null @@ -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) - diff --git a/anisotropy/shaping/periodic.py b/anisotropy/shaping/periodic.py index c1e05ad..dbd2667 100644 --- a/anisotropy/shaping/periodic.py +++ b/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 * -from numpy import pi, sqrt, cos, arccos, fix +import numpy as np from . import Shape @@ -13,85 +10,60 @@ class Periodic(Shape): self, alpha: float = None, r0: float = 1, - #L: float = None, - #radius: float = None, + L: 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): - """(Override) Parameter depending on the ``r0``. - """ - pass - - @property - def radius(self): + def radius(self) -> float: """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 - # ISSUE: MakeFilletAll : Fillet can't be computed on the given shape with the given radius. - # Temporary solution: degrade the precision (minRound <= analytical). + """Fillets radius beetween spheres. + + :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) diff --git a/anisotropy/shaping/primitives.py b/anisotropy/shaping/primitives.py new file mode 100644 index 0000000..bac3381 --- /dev/null +++ b/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 diff --git a/anisotropy/shaping/shape.py b/anisotropy/shaping/shape.py index 69ee2a2..346c47b 100644 --- a/anisotropy/shaping/shape.py +++ b/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 * -import numpy -from numpy import linalg -import os +from __future__ import annotations +from numpy import ndarray from os import PathLike -from pathlib import Path +import numpy as np +import netgen.occ as ng_occ +import pathlib -class ShapeError(Exception): - pass +from . import utils 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. - - :return: - Output, error messages and returncode + Path of the file. """ - out, err, returncode = "", "", 0 - path = Path(filename).resolve() + path = pathlib.Path(filename).resolve() ext = path.suffix[1: ] - try: - if ext == "step": - self.shape.WriteStep(path) - - else: - raise NotImplementedError(f"{ ext } is not supported") - - except NotImplementedError as e: - err = e - returncode = 1 + if ext == "step": + self.shape.WriteStep(str(path)) - except Exception as e: - err = e - returncode = 1 - - return out, err, returncode - + else: + raise NotImplementedError(f"Shape format '{ ext }' is not supported") + 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)) diff --git a/anisotropy/shaping/simple.py b/anisotropy/shaping/simple.py deleted file mode 100644 index b5b86de..0000000 --- a/anisotropy/shaping/simple.py +++ /dev/null @@ -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" - - diff --git a/anisotropy/shaping/utils.py b/anisotropy/shaping/utils.py new file mode 100644 index 0000000..72679eb --- /dev/null +++ b/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))