Mod: occ geometry

This commit is contained in:
L-Nafaryus 2021-11-04 18:26:13 +05:00
parent 06dcb26391
commit 6d90fa1c33
No known key found for this signature in database
GPG Key ID: C76D8DCD2727DBB7
15 changed files with 722 additions and 7 deletions

1
.gitignore vendored
View File

@ -4,6 +4,7 @@ logs/
storage/
temp/
env/
*output/
*.gz
*.xz
*.fls

View File

@ -2,8 +2,4 @@
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from netgen.occ import *
class Shape(object):
def __init__(self):
pass

View File

@ -0,0 +1,5 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.

View File

@ -2,6 +2,6 @@
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from anisotropy.samples.simple import Simple
from anisotropy.samples.bodyCentered import BodyCentered
from anisotropy.samples.faceCentered import FaceCentered
#from anisotropy.samples.simple import Simple
#from anisotropy.samples.bodyCentered import BodyCentered
#from anisotropy.samples.faceCentered import FaceCentered

View File

@ -0,0 +1,9 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .shape import Shape
from .periodic import Periodic
from .simple import Simple
from .faceCentered import FaceCentered
from .bodyCentered import BodyCentered

View File

@ -0,0 +1,151 @@
# -*- 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 . import Periodic
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 in [[1, 0, 0], [0, 0, 1]]:
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
if self.direction == [1, 0, 0]:
vertices = numpy.array([
(xl, 0, 0),
(0, yw, 0),
(0, yw, zh),
(xl, 0, zh)
])
extr = diag
elif self.direction == [0, 0, 1]:
vertices = numpy.array([
(0, yw, 0),
(xl, 0, 0),
(2 * xl, yw, 0),
(xl, 2 * yw, 0)
])
extr = height
elif self.direction == [1, 1, 1]:
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)
self.cell = inletface.Extrude(extr)
# Boundaries
symetryId = 0
for face in self.cell.faces:
fNorm = self.normal(face)
fAngle = self.angle(vecFlow, fNorm)
if fAngle == 0 and not face.center == inletface.center:
face.name = "outlet"
else:
face.name = f"symetry{ symetryId }"
symetryId += 1
# Main shape
self.shape = self.cell - self.lattice
for face in self.shape.faces:
if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
face.name = "wall"

View File

@ -0,0 +1,160 @@
# -*- 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 . import Periodic
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 in [[1, 0, 0], [0, 0, 1]]:
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
if self.direction == [1, 0, 0]:
vertices = numpy.array([
(0, 0, -zh),
(-xl, yw, -zh),
(-xl, yw, zh),
(0, 0, zh)
])
extr = length
elif self.direction == [0, 0, 1]:
vertices = numpy.array([
(0, 0, -zh),
(xl, yw, -zh),
(0, 2 * yw, -zh),
(-xl, yw, -zh)
])
extr = 2 * width
elif self.direction == [1, 1, 1]:
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)
self.cell = inletface.Extrude(extr)
# Boundaries
symetryId = 0
for face in self.cell.faces:
fNorm = self.normal(face)
fAngle = self.angle(vecFlow, fNorm)
if fAngle == 0 and not face.center == inletface.center:
face.name = "outlet"
else:
face.name = f"symetry{ symetryId }"
symetryId += 1
# Main shape
self.shape = self.cell - self.lattice
for face in self.shape.faces:
if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
face.name = "wall"

View File

@ -0,0 +1,97 @@
# -*- 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
from . import Shape
class Periodic(Shape):
def __init__(
self,
alpha: float = None,
r0: float = 1,
#L: float = None,
#radius: float = None,
filletsEnabled: bool = True,
#fillets: float = None,
gamma = None,
**kwargs
):
"""Constructor method.
:param alpha:
Spheres overlap parameter.
:param r0:
Initial spheres radius.
:param filletsEnabled:
Enable fillets beetween spheres.
"""
Shape.__init__(self)
# Parameters
self.alpha = alpha
self.r0 = r0
self.theta = 0.5 * pi
self.gamma = gamma or pi - 2 * (0.5 * self.theta)
self.filletsEnabled = filletsEnabled
self.filletsScale = 0.8
# 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):
"""Spheres 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).
rTol = 3
minRound = 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)

View File

@ -0,0 +1,66 @@
# -*- 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
class Shape(object):
def __init__(self):
self.groups = {}
self.shape = None
def export(self, filename: str):
"""Export a shape.
Supported formats: step.
:param filename:
Name of the file to store the given shape in.
:return:
Output, error messages and returncode
"""
out, err, returncode = "", "", 0
ext = os.path.splitext(filename)[1][1: ]
try:
if ext == "step":
self.shape.WriteStep(filename)
else:
raise NotImplementedError(f"{ ext } is not supported")
except NotImplementedError as e:
err = e
returncode = 1
except Exception as e:
err = e
returncode = 1
return out, err, returncode
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))

View File

@ -0,0 +1,145 @@
# -*- 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 . import Periodic
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 in [[1, 0, 0], [0, 0, 1]]:
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
if self.direction == [1, 0, 0]:
vertices = numpy.array([
(xl, 0, 0),
(0, yw, 0),
(0, yw, zh),
(xl, 0, zh)
])
extr = width
elif self.direction == [0, 0, 1]:
vertices = numpy.array([
(0, yw, 0),
(xl, 0, 0),
(2 * xl, yw, 0),
(xl, 2 * yw, 0)
])
extr = height
elif self.direction == [1, 1, 1]:
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)
self.cell = inletface.Extrude(extr)
# Boundaries
symetryId = 0
for face in self.cell.faces:
fNorm = self.normal(face)
fAngle = self.angle(vecFlow, fNorm)
if fAngle == 0 and not face.center == inletface.center:
face.name = "outlet"
else:
face.name = f"symetry{ symetryId }"
symetryId += 1
# Main shape
self.shape = self.cell - self.lattice
for face in self.shape.faces:
if face.name not in ["inlet", "outlet", *[ f"symetry{ n }" for n in range(6) ]]:
face.name = "wall"

View File

View File

@ -0,0 +1,5 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.

View File

@ -0,0 +1,22 @@
import os
import unittest
unittest.TestLoader.sortTestMethodsUsing = None
class TestFaceCentered(unittest.TestCase):
def setUp(self):
self.outputPath = os.path.join(os.path.abspath("."), "tests/test_shaping_output")
os.makedirs(self.outputPath, exist_ok = True)
def test_faceCentered_lattice(self):
from anisotropy.shaping import FaceCentered
fc = FaceCentered([1, 0, 0], alpha = 0.01, filletsEnabled = True)
fc.build()
fc.lattice.WriteStep(os.path.join(self.outputPath, "fc_lattice.step"))
def tearDown(self):
pass
if __name__ == "__main__":
unittest.main()

58
tests/test_shaping.py Normal file
View File

@ -0,0 +1,58 @@
import os
import unittest
unittest.TestLoader.sortTestMethodsUsing = None
class TestShaping(unittest.TestCase):
def setUp(self):
from anisotropy import shaping
self.shaping = shaping
self.outputPath = os.path.join(os.path.abspath("."), "tests/test_shaping_output")
os.makedirs(self.outputPath, exist_ok = True)
def test_simple(self):
simple100 = self.shaping.Simple(direction = [1, 0, 0], alpha = 0.01)
simple001 = self.shaping.Simple(direction = [0, 0, 1], alpha = 0.01)
simple111 = self.shaping.Simple(direction = [1, 1, 1], alpha = 0.01)
simple100.build()
simple001.build()
simple111.build()
simple100.export(os.path.join(self.outputPath, "simple100.step"))
simple001.export(os.path.join(self.outputPath, "simple001.step"))
simple111.export(os.path.join(self.outputPath, "simple111.step"))
def test_bodyCentered(self):
bodyCentered100 = self.shaping.BodyCentered(direction = [1, 0, 0], alpha = 0.01)
bodyCentered001 = self.shaping.BodyCentered(direction = [0, 0, 1], alpha = 0.01)
bodyCentered111 = self.shaping.BodyCentered(direction = [1, 1, 1], alpha = 0.01)
bodyCentered100.build()
bodyCentered001.build()
bodyCentered111.build()
bodyCentered100.export(os.path.join(self.outputPath, "bodyCentered100.step"))
bodyCentered001.export(os.path.join(self.outputPath, "bodyCentered001.step"))
bodyCentered111.export(os.path.join(self.outputPath, "bodyCentered111.step"))
def test_faceCentered(self):
faceCentered100 = self.shaping.FaceCentered(direction = [1, 0, 0], alpha = 0.01)
faceCentered001 = self.shaping.FaceCentered(direction = [0, 0, 1], alpha = 0.01)
faceCentered111 = self.shaping.FaceCentered(direction = [1, 1, 1], alpha = 0.01)
faceCentered100.build()
faceCentered001.build()
faceCentered111.build()
faceCentered100.export(os.path.join(self.outputPath, "faceCentered100.step"))
faceCentered001.export(os.path.join(self.outputPath, "faceCentered001.step"))
faceCentered111.export(os.path.join(self.outputPath, "faceCentered111.step"))
def tearDown(self):
pass
if __name__ == "__main__":
unittest.main()