L-Nafaryus/anisotropy
L-Nafaryus/anisotropy
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Mod: new tests

Browse Source 
Mod: working core runner
This commit is contained in:
L-Nafaryus 2021-11-19 23:35:12 +05:00
parent c56504f120
commit 8f4c664a74
No known key found for this signature in database
GPG Key ID: C76D8DCD2727DBB7
24 changed files with 282 additions and 2238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
anisotropy/core/__init__.py
View File

@ -2,4 +2,5 @@
# This file is part of anisotropy. # This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details. # License: GNU GPL version 3, see the file "LICENSE" for details.
from .config import Config, DefaultConfig
from .runner import UltimateRunner

4
anisotropy/core/config.py
View File

@ -96,14 +96,14 @@ class DefaultConfig(Config):
self.options = deepcopy(self.content["options"]) self.options = deepcopy(self.content["options"])
labels = ["simple", "bodyCentered", "faceCentered"] labels = ["simple", "bodyCentered", "faceCentered"]
thetas = array([[0.01, 0.28], [0.01, 0.17], [0.01, 0.13]], dtype = float) thetas = [[0.01, 0.28], [0.01, 0.17], [0.01, 0.13]]
for label, theta in zip(labels, thetas): for label, theta in zip(labels, thetas):
self.content["structures"].append({ self.content["structures"].append({
"label": label, "label": label,
"theta": theta, "theta": theta,
"thetaStep": 0.01, "thetaStep": 0.01,
"directions": array([[1, 0, 0], [0, 0, 1], [1, 1, 1]], dtype = float), "directions": [[1, 0, 0], [0, 0, 1], [1, 1, 1]],
"r0": 1, "r0": 1,
"filletsEnabled": True "filletsEnabled": True
}) })

405
anisotropy/core/database.py
View File

@ -1,405 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import os
import logging
from copy import deepcopy
from anisotropy import env
from anisotropy.core.utils import setupLogger
from anisotropy.core.models import (
db, JOIN,
Structure,
Mesh, SubMesh, MeshResult,
Flow, FlowApproximation, FlowResult
)
from peewee import OperationalError
logger = logging.getLogger(env["logger_name"])
#setupLogger(logger, logging.INFO, env["LOG"])
# NOTE: Temporary solution for multiprocessing mode when database is locked
def tryUntilDone(func):
def inner(*args, **kwargs):
done = False
while not done:
try:
ret = func(*args, **kwargs)
done = True
except OperationalError as e:
logger.error(e)
return ret
return inner
class Database(object):
def __init__(self, name: str, filepath: str):
self.name = name
self.filepath = filepath
self.__db = db
def setup(self):
os.makedirs(self.filepath, exist_ok = True)
fullpath = os.path.join(self.filepath, "{}.db".format(self.name))
self.__db.init(fullpath)
if not os.path.exists(fullpath):
self.__db.create_tables([
Structure,
Mesh,
SubMesh,
MeshResult,
Flow,
FlowApproximation,
FlowResult
])
def isempty(self) -> bool:
"""Checks DB for main table existence (Structure)
:return: True if exists
:rtype: bool
"""
query = Structure.select()
return not query.exists()
def load(self, structure_type: str, structure_direction: list, structure_theta: float) -> dict:
structureQuery = (
Structure
.select()
.where(
Structure.type == structure_type,
Structure.direction == str(structure_direction),
Structure.theta == structure_theta
)
)
params = {}
with self.__db.atomic():
if structureQuery.exists():
params["structure"] = structureQuery.dicts().get()
meshQuery = structureQuery.get().meshes
if meshQuery.exists():
params["mesh"] = meshQuery.dicts().get()
submeshQuery = meshQuery.get().submeshes
if submeshQuery.exists():
params["submesh"] = [ entry for entry in submeshQuery.dicts() ]
meshresultQuery = meshQuery.get().meshresults
if meshresultQuery.exists():
params["meshresult"] = meshresultQuery.dicts().get()
flowQuery = structureQuery.get().flows
if flowQuery.exists():
params["flow"] = flowQuery.dicts().get()
flowapproximationQuery = flowQuery.get().flowapproximations
if flowapproximationQuery.exists():
params["flowapproximation"] = flowapproximationQuery.dicts().get()
flowresultsQuery = flowQuery.get().flowresults
if flowresultsQuery.exists():
params["flowresult"] = flowresultsQuery.dicts().get()
else:
logger.error("Missed Structure table")
return params
def loadGeneral(self, type: str = None, direction: list = None, theta: float = None) -> list:
query = (
Structure
.select()
.order_by(Structure.type, Structure.direction, Structure.theta)
)
response = []
if type:
query = query.where(Structure.type == type)
if direction:
query = query.where(Structure.direction == str(direction))
if theta:
query = query.where(Structure.theta == theta)
for entry in query.dicts():
response.append({ "structure": entry })
return response
def update(self, params: dict):
if not params:
logger.error("Trying to update db from empty parameters")
return
query = (
Structure
.select(Structure, Mesh, Flow)
.join(
Mesh,
JOIN.INNER,
on = (Mesh.structure_id == Structure.structure_id)
)
.join(
Flow,
JOIN.INNER,
on = (Flow.structure_id == Structure.structure_id)
)
.where(
Structure.type == params["structure"]["type"],
Structure.direction == str(params["structure"]["direction"]),
Structure.theta == params["structure"]["theta"]
)
)
structureID = tryUntilDone(self._updateStructure)(params.get("structure", {}), query)
meshID = tryUntilDone(self._updateMesh)(params.get("mesh", {}), query, structureID)
for submeshParams in params.get("submesh", []):
tryUntilDone(self._updateSubMesh)(submeshParams, query, meshID)
tryUntilDone(self._updateMeshResult)(params.get("meshresult", {}), query, meshID)
flowID = tryUntilDone(self._updateFlow)(params.get("flow", {}), query, structureID)
tryUntilDone(self._updateFlowApproximation)(params.get("flowapproximation", {}), query, flowID)
tryUntilDone(self._updateFlowResult)(params.get("flowresult", {}), query, flowID)
def search(self, args: list):
result = {}
query = (
Structure
.select(Structure, Mesh, SubMesh, MeshResult, Flow, FlowApproximation, FlowResult)
.join(
Mesh,
JOIN.INNER,
on = (Mesh.structure_id == Structure.structure_id)
)
.join(
SubMesh,
JOIN.INNER,
on = (SubMesh.mesh_id == Mesh.mesh_id)
)
.join(
MeshResult,
JOIN.INNER,
on = (MeshResult.mesh_id == Mesh.mesh_id)
)
.join(
Flow,
JOIN.INNER,
on = (Flow.structure_id == Structure.structure_id)
)
.join(
FlowApproximation,
JOIN.INNER,
on = (FlowApproximation.flow_id == Flow.flow_id)
)
.join(
FlowResult,
JOIN.INNER,
on = (FlowResult.flow_id == Flow.flow_id)
)
)
for arg in args:
query = query.where(arg)
with self.__db.atomic():
if not self.isempty():
result = [ entry for entry in query.dicts() ]
else:
logger.error("Missed Structure table")
return result
def _updateStructure(self, src: dict, queryMain) -> int:
raw = deepcopy(src)
with self.__db.atomic():
if not queryMain.exists():
tabID = Structure.create(**raw)
else:
req = queryMain.dicts().get()
tabID = req["structure_id"]
query = (
Structure.update(**raw)
.where(
Structure.type == req["type"],
Structure.direction == str(req["direction"]),
Structure.theta == req["theta"]
)
)
query.execute()
return tabID
def _updateMesh(self, src: dict, queryMain, structureID) -> int:
raw = deepcopy(src)
with self.__db.atomic():
if not queryMain.exists():
tabID = Mesh.create(
structure_id = structureID,
**raw
)
else:
req = queryMain.dicts().get()
tabID = req["mesh_id"]
query = (
Mesh.update(**raw)
.where(
Mesh.structure_id == structureID
)
)
query.execute()
return tabID
def _updateSubMesh(self, src: dict, queryMain, meshID):
#if not src:
# return
raw = deepcopy(src)
with self.__db.atomic():
if not SubMesh.select().where(SubMesh.mesh_id == meshID).exists():
tabID = SubMesh.create(
mesh_id = meshID,
**raw
)
logger.debug(f"[ DB ] Created SubMesh entry { tabID }")
else:
query = (
SubMesh.update(**raw)
.where(
SubMesh.mesh_id == meshID,
SubMesh.name == src["name"]
)
)
query.execute()
def _updateMeshResult(self, src: dict, queryMain, meshID):
#if not src:
# return
raw = deepcopy(src)
with self.__db.atomic():
if not MeshResult.select().where(MeshResult.mesh_id == meshID).exists():
tabID = MeshResult.create(
mesh_id = meshID,
**raw
)
logger.debug(f"[ DB ] Created MeshResult entry { tabID }")
else:
query = (
MeshResult.update(**raw)
.where(
MeshResult.mesh_id == meshID
)
)
query.execute()
def _updateFlow(self, src: dict, queryMain, structureID):
raw = deepcopy(src)
with self.__db.atomic():
if not queryMain.exists():
tabID = Flow.create(
structure_id = structureID,
**raw
)
else:
req = queryMain.dicts().get()
tabID = req["flow_id"]
query = (
Flow.update(**raw)
.where(
Flow.structure_id == structureID
)
)
query.execute()
return tabID
def _updateFlowApproximation(self, src: dict, queryMain, flowID):
#if not src:
# return
raw = deepcopy(src)
with self.__db.atomic():
if not FlowApproximation.select().where(FlowApproximation.flow_id == flowID).exists():
tabID = FlowApproximation.create(
flow_id = flowID,
**raw
)
logger.debug(f"[ DB ] Created FlowApproximation entry { tabID }")
else:
query = (
FlowApproximation.update(**raw)
.where(
FlowApproximation.flow_id == flowID
)
)
query.execute()
def _updateFlowResult(self, src: dict, queryMain, flowID):
#if not src:
# return
raw = deepcopy(src)
with self.__db.atomic():
if not FlowResult.select().where(FlowResult.flow_id == flowID).exists():
tabID = FlowResult.create(
flow_id = flowID,
**raw
)
logger.debug(f"[ DB ] Created FlowResult entry { tabID }")
else:
query = (
FlowResult.update(**raw)
.where(
FlowResult.flow_id == flowID
)
)
query.execute()

185
anisotropy/core/models.py
View File

@ -1,185 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from peewee import (
SqliteDatabase, JOIN,
Model, Field,
AutoField, ForeignKeyField,
TextField, FloatField,
IntegerField, BooleanField,
TimeField, DateTimeField
)
import json
db = SqliteDatabase(
None,
pragmas = { "foreign_keys": 1 },
field_types = { "list": "text" }
)
class BaseModel(Model):
class Meta:
database = db
class ListField(Field):
field_type = "list"
def db_value(self, value):
return str(value)
def python_value(self, value):
pval = []
for entry in value[1 : -1].split(","):
try:
pval.append(float(entry))
except:
pval.append(entry.strip().replace("'", ""))
return pval
class JSONField(TextField):
def db_value(self, value):
return json.dumps(value)
def python_value(self, value):
if value is not None:
return json.loads(value)
class Execution(Model):
execution_id = AutoField()
date = DateTimeField()
executionTime = TimeField()
class Meta:
database = db
db_table = "executions"
class Structure(BaseModel):
structure_id = AutoField()
type = TextField()
direction = ListField()
theta = FloatField()
r0 = FloatField(null = True)
L = FloatField(null = True)
radius = FloatField(null = True)
filletsEnabled = BooleanField(null = True)
fillets = FloatField(null = True)
#path = TextField()
class Mesh(BaseModel):
mesh_id = AutoField()
structure_id = ForeignKeyField(Structure, backref = "meshes")
maxSize = FloatField(null = True)
minSize = FloatField(null = True)
fineness = IntegerField(null = True)
growthRate = FloatField(null = True)
nbSegPerEdge = FloatField(null = True)
nbSegPerRadius = FloatField(null = True)
chordalErrorEnabled = BooleanField(null = True)
chordalError = FloatField(null = True)
secondOrder = BooleanField(null = True)
optimize = BooleanField(null = True)
quadAllowed = BooleanField(null = True)
useSurfaceCurvature = BooleanField(null = True)
fuseEdges = BooleanField(null = True)
checkChartBoundary = BooleanField(null = True)
viscousLayers = BooleanField(null = True)
thickness = FloatField(null = True)
numberOfLayers = IntegerField(null = True)
stretchFactor = FloatField(null = True)
isFacesToIgnore = BooleanField(null = True)
facesToIgnore = ListField(null = True)
#faces = []
extrusionMethod = TextField(null = True)
class SubMesh(BaseModel):
submesh_id = AutoField()
mesh_id = ForeignKeyField(Mesh, backref = "submeshes")
name = TextField()
maxSize = FloatField(null = True)
minSize = FloatField(null = True)
fineness = IntegerField(null = True)
growthRate = FloatField(null = True)
nbSegPerEdge = FloatField(null = True)
nbSegPerRadius = FloatField(null = True)
chordalErrorEnabled = BooleanField(null = True)
chordalError = FloatField(null = True)
secondOrder = BooleanField(null = True)
optimize = BooleanField(null = True)
quadAllowed = BooleanField(null = True)
useSurfaceCurvature = BooleanField(null = True)
fuseEdges = BooleanField(null = True)
checkChartBoundary = BooleanField(null = True)
class MeshResult(BaseModel):
meshresult_id = AutoField()
mesh_id = ForeignKeyField(Mesh, backref = "meshresults")
surfaceArea = FloatField(null = True)
volume = FloatField(null = True)
volumeCell = FloatField(null = True)
elements = IntegerField(null = True)
edges = IntegerField(null = True)
faces = IntegerField(null = True)
volumes = IntegerField(null = True)
tetrahedrons = IntegerField(null = True)
prisms = IntegerField(null = True)
pyramids = IntegerField(null = True)
meshStatus = TextField(null = True, default = "Idle")
meshCalculationTime = TimeField(null = True)
class Flow(BaseModel):
flow_id = AutoField()
structure_id = ForeignKeyField(Structure, backref = "flows")
scale = ListField(null = True)
pressure = JSONField(null = True)
velocity = JSONField(null = True)
transportProperties = JSONField(null = True)
class FlowApproximation(BaseModel):
flow_approximation_id = AutoField()
flow_id = ForeignKeyField(Flow, backref = "flowapproximations")
pressure = JSONField(null = True)
velocity = JSONField(null = True)
transportProperties = JSONField(null = True)
class FlowResult(BaseModel):
flowresult_id = AutoField()
flow_id = ForeignKeyField(Flow, backref = "flowresults")
flowRate = FloatField(null = True)
porosity = FloatField(null = True)
permeability = FloatField(null = True)
flowStatus = TextField(null = True, default = "Idle")
flowCalculationTime = TimeField(null = True)

44
anisotropy/core/runner.py
View File

@ -28,44 +28,48 @@ class UltimateRunner(object):
self.mesh = None self.mesh = None
self.flow = None self.flow = None
def casePath(self): def casepath(self):
case = self.config.cases[0] params = self.config.cases[0]
return path.join( return path.abspath(path.join(
self.config["build"], self.config["build"],
case["label"], params["label"],
"direction-[{},{},{}]".format(*[ str(d) for d in case["direction"] ]), "direction-[{},{},{}]".format(*[ str(d) for d in params["direction"] ]),
"theta-{}".format(case["theta"]) "theta-{}".format(params["theta"])
) ))
def computeShape(self): def computeShape(self):
case = self.config.cases[0] params = self.config.cases[0]
filename = "shape.step" filename = "shape.step"
match case["label"]: match params["label"]:
case "simple": case "simple":
self.shape = Simple(case["direction"]) self.shape = Simple(params["direction"])
case "bodyCentered": case "bodyCentered":
self.shape = BodyCentered(case["direction"]) self.shape = BodyCentered(params["direction"])
case "faceCentered": case "faceCentered":
self.shape = FaceCentered(case["direction"]) self.shape = FaceCentered(params["direction"])
self.shape.build() self.shape.build()
self.shape.export(path.join(case, filename))
os.makedirs(self.casepath(), exist_ok = True)
self.shape.export(path.join(self.casepath(), filename))
def computeMesh(self): def computeMesh(self):
case = self.config.cases[0] params = self.config.cases[0]
filename = "mesh.mesh" filename = "mesh.mesh"
self.mesh = Mesh(self.shape.shape) self.mesh = Mesh(self.shape.shape)
self.mesh.build() self.mesh.build()
self.mesh.export(path.join(case, filename))
os.makedirs(self.casepath(), exist_ok = True)
self.mesh.export(path.join(self.casepath(), filename))
def computeFlow(self): def computeFlow(self):
case = self.config.cases[0] params = self.config.cases[0]
flow = OnePhaseFlow() flow = OnePhaseFlow(path = self.casepath())
# initial 43 unnamed patches -> # initial 43 unnamed patches ->
# 6 named patches (inlet, outlet, wall, symetry0 - 3/5) -> # 6 named patches (inlet, outlet, wall, symetry0 - 3/5) ->
@ -75,13 +79,15 @@ class UltimateRunner(object):
patches = {} patches = {}
for n, patch in enumerate(self.shape.shape.faces): for n, patch in enumerate(self.shape.shape.faces):
# shifted index
n += 1
name = patch.name name = patch.name
if patches.get(name): if patches.get(name):
patches[name].append(n) patches[name].append(f"patch{ n }")
else: else:
patches[name] = [n] patches[name] = [ f"patch{ n }" ]
for name in patches.keys(): for name in patches.keys():
match name: match name:

4
anisotropy/database/database.py
View File

@ -20,13 +20,13 @@ class Database(object):
def setup(self): def setup(self):
path = os.path.abspath(self.filename) path = os.path.abspath(self.filename)
os.makedirs(path, exist_ok = True) #os.makedirs(path, exist_ok = True)
self.database.init(path) self.database.init(path)
if not os.path.exists(path): if not os.path.exists(path):
self.database.create_tables([Execution])
self.database.create_tables([ self.database.create_tables([
Execution,
Physics, Physics,
Shape, Shape,
Mesh, Mesh,

10
anisotropy/database/models.py
View File

@ -26,7 +26,7 @@ class Execution(Model):
class Meta: class Meta:
database = sqliteDB database = sqliteDB
db_table = "executions" table_name = "executions"
class Physics(Model): class Physics(Model):
@ -43,7 +43,7 @@ class Physics(Model):
class Meta: class Meta:
database = sqliteDB database = sqliteDB
db_table = "physics" table_name = "physics"
depends_on = Execution depends_on = Execution
@ -64,7 +64,7 @@ class Shape(Model):
class Meta: class Meta:
database = sqliteDB database = sqliteDB
db_table = "shapes" table_name = "shapes"
depends_on = Execution depends_on = Execution
@ -85,7 +85,7 @@ class Mesh(Model):
class Meta: class Meta:
database = sqliteDB database = sqliteDB
db_table = "meshes" table_name = "meshes"
depends_on = Execution depends_on = Execution
@ -98,7 +98,7 @@ class Flow(Model):
class Meta: class Meta:
database = sqliteDB database = sqliteDB
db_table = "flows" table_name = "flows"
depends_on = Execution depends_on = Execution

2
anisotropy/openfoam/__init__.py
View File

@ -3,7 +3,7 @@
# License: GNU GPL version 3, see the file "LICENSE" for details. # License: GNU GPL version 3, see the file "LICENSE" for details.
from .meshConversion import ideasUnvToFoam from .meshConversion import ideasUnvToFoam, netgenNeutralToFoam
from .meshManipulation import createPatch, transformPoints, checkMesh, renumberMesh from .meshManipulation import createPatch, transformPoints, checkMesh, renumberMesh
from .miscellaneous import foamDictionary from .miscellaneous import foamDictionary
from .parallelProcessing import decomposePar from .parallelProcessing import decomposePar

5
anisotropy/salomepl/__init__.py
View File

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

216
anisotropy/salomepl/geometry.py
View File

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

304
anisotropy/salomepl/mesh.py
View File

@ -1,304 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import logging
SMESH_IMPORTED = False
try:
import SMESH
from salome.smesh import smeshBuilder
SMESH_IMPORTED = True
except:
pass
class Mesh(object):
def __init__(self, geom):
# Mesh module
if not SMESH_IMPORTED:
raise ImportError("Cannot find the salome mesh modules.")
else:
self.smesh = smeshBuilder.New()
self.smeshBuilder = smeshBuilder
# General attributes
self.geom = geom
self.mesh = self.smesh.Mesh(self.geom.shape, self.geom.name)
def algo3d(self, algo, type = "tetrahedron"):
smeshAlgo = self.mesh.__dict__.get(type.capitalize())
self.meshAlgorithm3d = algo()
self.meshAlgorithm3d.initialize(smeshAlgo(algo = self.meshAlgorithm3d.key))
self.mesh.AddHypothesis(self.meshAlgorithm3d.hypo)
return self.meshAlgorithm3d
def algo2d(self, algo, type = "triangle"):
smeshAlgo = self.mesh.__dict__.get(type.capitalize())
self.meshAlgorithm2d = algo()
self.meshAlgorithm2d.initialize(smeshAlgo(algo = self.meshAlgorithm2d.key))
self.mesh.AddHypothesis(self.meshAlgorithm2d.hypo)
return self.meshAlgorithm2d
def algo1d(self, algo, type = "segment"):
smeshAlgo = self.mesh.__dict__.get(type.capitalize())
self.meshAlgorithm1d = algo()
self.meshAlgorithm1d.initialize(smeshAlgo(algo = self.meshAlgorithm1d.key))
self.mesh.AddHypothesis(self.meshAlgorithm1d.hypo)
return self.meshAlgorithm1d
def createGroups(self, prefix = None):
prefix = prefix or ""
for group in self.shape.groups:
name = group.GetName()
if name:
name = prefix + name
self.mesh.GroupOnGeom(group, name, SMESH.FACE)
def compute(self):
"""Compute mesh.
"""
isDone = self.mesh.Compute()
out = ""
err = self.mesh.GetComputeErrors()
returncode = int(not isDone)
return out, err, returncode
def stats(self):
return {
"elements": self.mesh.NbElements(),
"edges": self.mesh.NbEdges(),
"faces": self.mesh.NbFaces(),
"volumes": self.mesh.NbVolumes(),
"tetrahedrons": self.mesh.NbTetras(),
"prisms": self.mesh.NbPrisms(),
"pyramids": self.mesh.NbPyramids()
}
def removePyramids(self):
if self.mesh.NbPyramids() > 0:
pyramidCriterion = smesh.GetCriterion(
SMESH.VOLUME,
SMESH.FT_ElemGeomType,
SMESH.FT_Undefined,
SMESH.Geom_PYRAMID
)
pyramidGroup = self.mesh.MakeGroupByCriterion("pyramids", pyramidCriterion)
pyramidVolumes = self.mesh.GetIDSource(pyramidGroup.GetIDs(), SMESH.VOLUME)
self.mesh.SplitVolumesIntoTetra(pyramidVolumes, smesh.Hex_5Tet)
self.mesh.RemoveGroup(pyramidGroup)
self.mesh.RenumberElements()
def export(
filename: str
):
"""Export a mesh.
Supported formats: unv.
:param filename:
Name of the file to store the given mesh in.
:return:
Output, error messages and returncode
"""
out, err, returncode = "", "", 0
ext = os.path.splitext(filename)[1][1: ]
try:
if ext == "unv":
self.mesh.ExportUNV(self.mesh, filename)
else:
raise NotImplementedError(f"{ ext } is not supported")
except NotImplementedError as e:
err = e
returncode = 1
except Exception as e:
err = e.details.text
returncode = 1
return out, err, returncode
class MeshAlgorithm(object):
pass
class AlgorithmHypothesis(object):
pass
class Netgen3D(MeshAlgorithm):
"""
MaxElementVolume
Parameters
ViscousLayers
"""
def __init__(self, **kwargs):
self.key = smeshBuilder.NETGEN_3D
def initialize(self, algo, hypo): #thesises: list):
self.algo = algo
#self.hypo = self.algo.Parameters()
#for hypo in hypothesises:
self.hypo = self.__dict__[hypo.__name__]()
class ViscousLayers(AlgorithmHypothesis):
def __init__(self,
algo,
thickness = 1,
numberOfLayers = 1,
stretchFactor = 0,
faces = [],
isFacesToIgnore = True,
extrMethod = "SURF_OFFSET_SMOOTH",
**kwargs
):
extrusionMethod = dict(
SURF_OFFSET_SMOOTH = smeshBuilder.SURF_OFFSET_SMOOTH,
FACE_OFFSET = smeshBuilder.FACE_OFFSET,
NODE_OFFSET = smeshBuilder.NODE_OFFSET,
).get(extrMethod, smeshBuilder.SURF_OFFSET_SMOOTH)
self.hypo = self.algo.ViscousLayers(
thickness,
numberOfLayers,
stretchFactor,
faces,
isFacesToIgnore,
extrusionMethod
)
class Parameters(AlgorithmHypothesis):
def __init__(self, algo):
self.hypo = self.algo.Parameters()
@property
def minSize(self):
return self.hypo.GetMinSize()
@minSize.setter
def minSize(self, value):
self.hypo.SetMinSize(value)
@property
def maxSize(self):
return self.hypo.GetMaxSize()
@maxSize.setter
def maxSize(self, value):
self.hypo.SetMaxSize(value)
@property
def fineness(self):
return self.hypo.GetFineness()
@fineness.setter
def fineness(self, value):
self.hypo.SetFineness(value)
@property
def growthRate(self):
return self.hypo.GetGrowthRate()
@growthRate.setter
def growthRate(self, value):
self.hypo.SetGrowthRate(value)
@property
def nbSegPerEdge(self):
return self.hypo.GetNbSegPerEdge()
@nbSegPerEdge.setter
def nbSegPerEdge(self, value):
self.hypo.SetNbSegPerEdge(value)
@property
def nbSegPerRadius(self):
return self.hypo.GetNbSegPerRadius()
@nbSegPerRadius.setter
def nbSegPerRadius(self, value):
self.hypo.SetNbSegPerRadius(value)
@property
def chordalErrorEnabled(self):
return self.hypo.GetChordalErrorEnabled()
@chordalErrorEnabled.setter
def chordalErrorEnabled(self, value):
self.hypo.SetChordalErrorEnabled(value)
@property
def chordalError(self):
return self.hypo.GetChordalError()
@chordalError.setter
def chordalError(self, value):
self.hypo.SetChordalError(value)
@property
def secondOrder(self):
return self.hypo.GetSecondOrder()
@secondOrder.setter
def secondOrder(self, value):
self.hypo.SetSecondOrder(value)
@property
def optimize(self):
return self.hypo.GetOptimize()
@optimize.setter
def optimize(self, value):
self.hypo.SetOptimize(value)
@property
def quadAllowed(self):
return self.hypo.GetQuadAllowed()
@quadAllowed.setter
def quadAllowed(self, value):
self.hypo.SetQuadAllowed(value)
@property
def useSurfaceCurvature(self):
return self.hypo.GetUseSurfaceCurvature()
@useSurfaceCurvature.setter
def useSurfaceCurvature(self, value):
self.hypo.SetUseSurfaceCurvature(value)
@property
def fuseEdges(self):
return self.hypo.GetFuseEdges()
@fuseEdges.setter
def fuseEdges(self, value):
self.hypo.SetFuseEdges(value)
@property
def checkChartBoundary(self):
return self.hypo.GetCheckChartBoundary()
@checkChartBoundary.setter
def GetCheckChartBoundary(self, value):
self.hypo.SetCheckChartBoundary(value)
class MEFISTO(MeshAlgorithm):
pass

112
anisotropy/salomepl/utils.py
View File

@ -1,112 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import subprocess
import logging
import sys, os
import re
class SalomeNotFound(Exception):
pass
class SalomeManager(object):
def __init__(self):
self.__port = None
self.__lastproc = None
def runner(self, cmdargs: list, **kwargs):
timeout = kwargs.pop("timeout") if kwargs.get("timeout") else None
try:
with subprocess.Popen(
cmdargs,
stdout = kwargs.pop("stdout") if kwargs.get("stdout") else subprocess.PIPE,
stderr = kwargs.pop("stdout") if kwargs.get("stderr") else subprocess.PIPE,
**kwargs
) as proc:
self.__lastproc = proc
out, err = proc.communicate(timeout = timeout)
except FileNotFoundError:
raise SalomeNotFound()
return out, err, proc.returncode
def port(self) -> int:
out, err, returncode = self.runner(["salome", "start", "--print-port"], text = True)
if returncode == 0:
reg = re.search("(?!port:)([0-9]+)", out)
if reg:
return int(reg.group())
return 2810
def version(self) -> int:
out, err, returncode = self.runner(["salome", "--version"], text = True)
return out.strip().split(" ")[-1]
def kill(self, port: int = None):
return self.runner(["salome", "kill", str(self.__port or port)])
def killall(self):
return self.runner(["salome", "killall"])
def execute(self, scriptpath: str, *args, root: str = None, logpath: str = None, timeout: int = None, **kwargs):
if not root:
root = os.environ["HOME"]
# ISSUE: salome removes commas from string list
args = list(args)
args.insert(1, root)
salomeargs = "args:"
salomeargs += ",".join([ '"{}"'.format(str(arg)) for arg in args ])
if kwargs:
salomeargs += "," + ",".join([ '{}="{}"'.format(k, v) for k, v in kwargs.items() ])
###
self.__port = self.port()
cmd = [
"salome",
"start",
"-t",
"--shutdown-servers=1",
"--port", str(self.__port),
scriptpath,
salomeargs
]
try:
out, err, returncode = self.runner(cmd, timeout = timeout)
except subprocess.TimeoutExpired:
lastproc = self.__lastproc
self.kill()
out, err = lastproc.communicate()
returncode = lastproc.returncode
if logpath:
os.makedirs(logpath, exist_ok = True)
with open(os.path.join(logpath, "salome.log"), "wb") as io:
io.write(out)
io.write(err)
return str(out, "utf-8"), str(err, "utf-8"), returncode

7
anisotropy/samples/__init__.py
View File

@ -1,7 +0,0 @@
# -*- coding: utf-8 -*-
# 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

251
anisotropy/samples/bodyCentered.py
View File

@ -1,251 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from anisotropy.salomepl.geometry import StructureGeometry
from numpy import pi, sqrt, cos, arccos, fix
import logging
class BodyCentered(StructureGeometry):
@property
def name(self):
"""Shape name.
"""
return "bodyCentered"
@property
def L(self):
return self.r0 * 4 / sqrt(3)
@property
def thetaMin(self):
return 0.01
@property
def thetaMax(self):
return 0.18
@property
def fillets(self):
analytical = self.r0 * (sqrt(2) / sqrt(1 - cos(pi - 2 * arccos(sqrt(2 / 3)))) -
1 / (1 - self.theta))
# 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.filletScale
def build(self):
###
# Pore Cell
##
if self.direction in [[1, 0, 0], [0, 0, 1]]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = 2 * self.r0
width = self.L / 2
diag = self.L * sqrt(2)
height = self.L
point = []
xl = sqrt(diag ** 2 + diag ** 2) * 0.5
yw = xl
zh = height
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
spos1 = (0, 0, 0)
spos2 = (0, 0, 0)
###
# Bounding box
##
if self.direction == [1, 0, 0]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(xl, 0, 0)
sk.addPointsAbsolute(0, yw, 0)
sk.addPointsAbsolute(0, yw, zh)
sk.addPointsAbsolute(xl, 0, zh)
sk.addPointsAbsolute(xl, 0, 0)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, diag)
elif self.direction == [0, 0, 1]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(0, yw, 0)
sk.addPointsAbsolute(xl, 0, 0)
sk.addPointsAbsolute(2 * xl, yw, 0)
sk.addPointsAbsolute(xl, 2 * yw, 0)
sk.addPointsAbsolute(0, yw, 0)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, zh)
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
elif self.direction == [1, 1, 1]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = 2 * self.r0
width = self.L / 2
diag = self.L * sqrt(2)
height = diag / 3
point = []
xl, yw, zh = -self.L / 4, -self.L / 4, -self.L / 4
point.append((self.L / 3 + xl, self.L / 3 + yw, 4 * self.L / 3 + zh))
point.append((self.L + xl, 0 + yw, self.L + zh))
point.append((4 * self.L / 3 + xl, self.L / 3 + yw, self.L / 3 + zh))
point.append((self.L + xl, self.L + yw, 0 + zh))
point.append((self.L / 3 + xl, 4 * self.L / 3 + yw, self.L / 3 + zh))
point.append((0 + xl, self.L + yw, self.L + zh))
point.append((self.L / 3 + xl, self.L / 3 + yw, 4 * self.L / 3 + zh))
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
spos1 = (0, 0, 0)
spos2 = (0, 0, 0)
###
# Bounding box
##
sk = self.geo.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = self.geo.MakeFaceWires([sk.wire()], False)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
else:
raise Exception(f"Direction { self.direction } is not implemented")
###
# Lattice
##
ox = self.geo.MakeVectorDXDYDZ(1, 0, 0)
oy = self.geo.MakeVectorDXDYDZ(0, 1, 0)
oz = self.geo.MakeVectorDXDYDZ(0, 0, 1)
xy = self.geo.MakeVectorDXDYDZ(1, 1, 0)
xmy = self.geo.MakeVectorDXDYDZ(1, -1, 0)
grain = self.geo.MakeSpherePntR(self.geo.MakeVertex(*spos1), self.radius)
lattice1 = self.geo.MakeMultiTranslation2D(grain, ox, self.L, xn, oy, self.L, yn)
lattice1 = self.geo.MakeMultiTranslation1D(lattice1, oz, self.L, zn)
lattice2 = self.geo.MakeTranslation(lattice1, -0.5 * self.L, -0.5 * self.L, -0.5 * self.L)
grains = self.geo.ExtractShapes(lattice1, self.geo.ShapeType["SOLID"], True)
grains += self.geo.ExtractShapes(lattice2, self.geo.ShapeType["SOLID"], True)
grains = self.geo.MakeFuseList(grains, False, False)
grains = self.geo.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
grains = self.geo.MakeFilletAll(grains, self.fillets * scale)
self.shapeLattice = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
###
# Shape
##
self.shape = self.geo.MakeCutList(poreCell, [grains])
self.shape = self.geo.MakeScaleTransform(self.shape, oo, 1 / scale, theName = self.name)
isValid, _ = self.isValid()
if not isValid:
self.heal()
###
# Groups
#
# inlet, outlet, simetry(N), strips(optional), wall
##
self.groups = []
groupAll = self.createGroupAll(self.shape)
self.groups.append(self.createGroup(self.shape, inletface, "inlet", [grains], 1 / scale))
self.groups.append(self.createGroup(self.shape, outletface, "outlet", [grains], 1 / scale))
for n, face in enumerate(symetryface):
self.groups.append(self.createGroup(self.shape, face, f"symetry{ n }", [grains], 1 / scale))
if self.filletsEnabled:
shapeShell = self.geo.ExtractShapes(self.shape, self.geo.ShapeType["SHELL"], True)[0]
self.groups.append(self.createGroup(self.shape, shapeShell, "strips", self.groups + [grainsOrigin]))
self.groups.append(self.geo.CutListOfGroups([groupAll], self.groups, theName = "wall"))
class BodyCenteredMesh(Mesh):
def build(self):
algo2d = self.mesh.Triangle(algo = self.smeshBuilder.NETGEN_1D2D)
hypo2d = algo2d.Parameters()
hypo2d.SetMaxSize(0.1)
hypo2d.SetMinSize(0.001)
hypo2d.SetFineness(5)
hypo2d.SetGrowthRate(0.3)
hypo2d.SetNbSegPerEdge(2)
hypo2d.SetNbSegPerRadius(3)
hypo2d.SetChordalErrorEnabled(True)
hypo2d.SetChordalError(0.05)
hypo2d.SetOptimize(True)
hypo2d.SetUseSurfaceCurvature(True)
algo3d = self.mesh.Tetrahedron(algo = self.smeshBuilder.NETGEN_3D)
#hypo3d = algo3d.Parameters()
#faces = [ group for group in self.geom.groups if group.GetName() in ["inlet", "outlet"] ]
#hypo3dVL = algo3d.ViscousLayers(...)

255
anisotropy/samples/faceCentered.py
View File

@ -1,255 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from anisotropy.salomepl.geometry import StructureGeometry
from anisotropy.salomepl.mesh import Mesh
from numpy import pi, sqrt, fix
import logging
class FaceCentered(StructureGeometry):
@property
def name(self):
"""Shape name.
"""
return "faceCentered"
@property
def L(self):
return self.r0 * 4 / sqrt(2)
@property
def thetaMin(self):
return 0.01
@property
def thetaMax(self):
return 0.13
@property
def fillets(self):
analytical = self.r0 * (2 * sqrt(3) / 3 - 1 / (1 - self.theta))
# 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.filletScale
def build(self):
###
# Pore Cell
##
if self.direction in [[1, 0, 0], [0, 0, 1]]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = 2 * self.r0
width = self.L / 2
diag = self.L * sqrt(3)
height = diag / 3
point = []
xl = sqrt(length ** 2 + length ** 2) * 0.5
yw = xl
zh = width
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
spos1 = (-width * (xn - 1), 0, -width * (zn - 2))
spos2 = (-width * xn, 0, -width * (zn - 1))
###
# Bounding box
##
if self.direction == [1, 0, 0]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(0, 0, -zh)
sk.addPointsAbsolute(-xl, yw, -zh)
sk.addPointsAbsolute(-xl, yw, zh)
sk.addPointsAbsolute(0, 0, zh)
sk.addPointsAbsolute(0, 0, -zh)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, length)
elif self.direction == [0, 0, 1]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(0, 0, -zh)
sk.addPointsAbsolute(xl, yw, -zh)
sk.addPointsAbsolute(0, 2 * yw, -zh)
sk.addPointsAbsolute(-xl, yw, -zh)
sk.addPointsAbsolute(0, 0, -zh)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, 2 * zh)
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
elif self.direction == [1, 1, 1]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = 2 * self.r0
width = self.L / 2
diag = self.L * sqrt(3)
height = diag / 3
point = []
xl, yw, zh = -(xn - 2) * self.L / 3, -(yn - 2) * self.L / 3, -(zn - 2) * self.L / 3
point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh))
point.append((0 + xl, -width + yw, 0 + zh))
point.append((width / 3 + xl, -2 * width / 3 + yw, -2 * width / 3 + zh))
point.append((0 + xl, 0 + yw, -width + zh))
point.append((-2 * width / 3 + xl, width / 3 + yw, -2 * width / 3 + zh))
point.append((-width + xl, 0 + yw, 0 + zh))
point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh))
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
spos1 = (-width * (xn - 1), 0, -width * (zn - 2))
spos2 = (-width * xn, 0, -width * (zn - 1))
###
# Bounding box
##
sk = self.geo.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = self.geo.MakeFaceWires([sk.wire()], False)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
else:
raise Exception(f"Direction { self.direction } is not implemented")
###
# Grains
##
ox = self.geo.MakeVectorDXDYDZ(1, 0, 0)
oy = self.geo.MakeVectorDXDYDZ(0, 1, 0)
oz = self.geo.MakeVectorDXDYDZ(0, 0, 1)
xy = self.geo.MakeVectorDXDYDZ(1, 1, 0)
xmy = self.geo.MakeVectorDXDYDZ(1, -1, 0)
grain = self.geo.MakeSpherePntR(self.geo.MakeVertex(*spos1), self.radius)
lattice1 = self.geo.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn)
lattice1 = self.geo.MakeMultiTranslation1D(lattice1, oz, self.L, zn - 1)
grain = self.geo.MakeSpherePntR(self.geo.MakeVertex(*spos2), self.radius)
lattice2 = self.geo.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
lattice2 = self.geo.MakeMultiTranslation1D(lattice2, oz, self.L, zn)
grains = self.geo.ExtractShapes(lattice1, self.geo.ShapeType["SOLID"], True)
grains += self.geo.ExtractShapes(lattice2, self.geo.ShapeType["SOLID"], True)
grains = self.geo.MakeFuseList(grains, False, False)
grains = self.geo.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
grains = self.geo.MakeFilletAll(grains, self.fillets * scale)
self.shapeLattice = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
###
# Shape
##
#poreCell = self.geo.LimitTolerance(poreCell, 1e-12)
#grains = self.geo.LimitTolerance(grains, 1e-12)
self.shape = self.geo.MakeCutList(poreCell, [grains])
self.shape = self.geo.MakeScaleTransform(self.shape, oo, 1 / scale, theName = self.name)
isValid, msg = self.isValid()
if not isValid:
logging.warning(msg)
self.heal()
###
# Groups
#
# inlet, outlet, simetry(N), strips(optional), wall
##
self.groups = []
groupAll = self.createGroupAll(self.shape)
self.groups.append(self.createGroup(self.shape, inletface, "inlet", [grains], 1 / scale))
self.groups.append(self.createGroup(self.shape, outletface, "outlet", [grains], 1 / scale))
for n, face in enumerate(symetryface):
self.groups.append(self.createGroup(self.shape, face, f"symetry{ n }", [grains], 1 / scale))
if self.filletsEnabled:
shapeShell = self.geo.ExtractShapes(self.shape, self.geo.ShapeType["SHELL"], True)[0]
self.groups.append(self.createGroup(self.shape, shapeShell, "strips", self.groups + [grainsOrigin]))
self.groups.append(self.geo.CutListOfGroups([groupAll], self.groups, theName = "wall"))
class FaceCenteredMesh(Mesh):
def build(self):
algo2d = self.mesh.Triangle(algo = self.smeshBuilder.NETGEN_1D2D)
hypo2d = algo2d.Parameters()
hypo2d.SetMaxSize(0.1)
hypo2d.SetMinSize(0.001)
hypo2d.SetFineness(5)
hypo2d.SetGrowthRate(0.3)
hypo2d.SetNbSegPerEdge(2)
hypo2d.SetNbSegPerRadius(3)
hypo2d.SetChordalErrorEnabled(True)
hypo2d.SetChordalError(0.05)
hypo2d.SetOptimize(True)
hypo2d.SetUseSurfaceCurvature(True)
algo3d = self.mesh.Tetrahedron(algo = self.smeshBuilder.NETGEN_3D)
#hypo3d = algo3d.Parameters()
#faces = [ group for group in self.geom.groups if group.GetName() in ["inlet", "outlet"] ]
#hypo3dVL = algo3d.ViscousLayers(...)

405
anisotropy/samples/simple.py
View File

@ -1,405 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from numpy import pi, sqrt, fix
import logging
from anisotropy.salomepl.geometry import StructureGeometry
from anisotropy.salomepl.mesh import Mesh
from anisotropy.openfoam.presets import (
ControlDict, FvSchemes, FvSolution,
TransportProperties, TurbulenceProperties, CreatePatchDict,
P, U
)
from anisotropy.openfoam.foamcase import FoamCase
import anisotropy.openfoam as openfoam
class _Geometry(StructureGeometry):
@property
def name(self):
"""Shape name.
"""
return "simple"
@property
def L(self):
return 2 * self.r0
@property
def thetaMin(self):
return 0.01
@property
def thetaMax(self):
return 0.28
@property
def fillets(self):
analytical = self.r0 * (sqrt(2) - 1 / (1 - self.theta))
# 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.filletScale
def build(self):
###
# Pore Cell
##
if self.direction in [[1, 0, 0], [0, 0, 1]]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = self.L * sqrt(2)
width = self.L * sqrt(2)
height = self.L
xl = sqrt(length ** 2 * 0.5)
yw = xl
zh = height
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
###
# Bounding box
##
if self.direction == [1, 0, 0]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(xl, 0, 0)
sk.addPointsAbsolute(0, yw, 0)
sk.addPointsAbsolute(0, yw, zh)
sk.addPointsAbsolute(xl, 0, zh)
sk.addPointsAbsolute(xl, 0, 0)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, width)
elif self.direction == [0, 0, 1]:
sk = self.geo.Sketcher3D()
sk.addPointsAbsolute(0, yw, 0)
sk.addPointsAbsolute(xl, 0, 0)
sk.addPointsAbsolute(2 * xl, yw, 0)
sk.addPointsAbsolute(xl, 2 * yw, 0)
sk.addPointsAbsolute(0, yw, 0)
inletface = self.geo.MakeFaceWires([sk.wire()], True)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, height)
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
elif self.direction == [1, 1, 1]:
###
# Parameters
##
xn, yn, zn = 3, 3, 3
length = self.L * sqrt(2)
width = self.L * sqrt(2)
height = self.L
point = []
xl, yw, zh = -self.L - self.L / 6, -self.L - self.L / 6, -self.L / 6
point.append((self.L + xl, self.L + yw, self.L + zh))
point.append((5 * self.L / 3 + xl, 2 * self.L / 3 + yw, 2 * self.L / 3 + zh))
point.append((2 * self.L + xl, self.L + yw, 0 + zh))
point.append((5 * self.L / 3 + xl, 5 * self.L / 3 + yw, -self.L / 3 + zh))
point.append((self.L + xl, 2 * self.L + yw, 0 + zh))
point.append((2 * self.L / 3 + xl, 5 * self.L / 3 + yw, 2 * self.L / 3 + zh))
point.append((self.L + xl, self.L + yw, self.L + zh))
scale = 100
oo = self.geo.MakeVertex(0, 0, 0)
###
# Bounding box
##
sk = self.geo.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = self.geo.MakeFaceWires([sk.wire()], False)
vecflow = self.geo.GetNormal(inletface)
poreCell = self.geo.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
self.shapeCell = poreCell
inletface = self.geo.MakeScaleTransform(inletface, oo, scale)
poreCell = self.geo.MakeScaleTransform(poreCell, oo, scale)
faces = self.geo.ExtractShapes(poreCell, self.geo.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = self.geo.GetNormal(face)
angle = round(self.geo.GetAngle(norm, vecflow), 0)
if (angle == 0 or angle == 180) and not face == inletface:
outletface = face
else:
symetryface.append(face)
else:
raise Exception(f"Direction { self.direction } is not implemented")
###
# Grains
##
ox = self.geo.MakeVectorDXDYDZ(1, 0, 0)
oy = self.geo.MakeVectorDXDYDZ(0, 1, 0)
oz = self.geo.MakeVectorDXDYDZ(0, 0, 1)
grain = self.geo.MakeSphereR(self.radius)
lattice = self.geo.MakeMultiTranslation2D(grain, ox, self.L, xn, oy, self.L, yn)
lattice = self.geo.MakeMultiTranslation1D(lattice, oz, self.L, zn)
grains = self.geo.ExtractShapes(lattice, self.geo.ShapeType["SOLID"], True)
grains = self.geo.MakeFuseList(grains, False, False)
grains = self.geo.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
grains = self.geo.MakeFilletAll(grains, self.fillets * scale)
self.shapeLattice = self.geo.MakeScaleTransform(grains, oo, 1 / scale)
###
# Shape
##
self.shape = self.geo.MakeCutList(poreCell, [grains])
self.shape = self.geo.MakeScaleTransform(self.shape, oo, 1 / scale, theName = self.name)
isValid, msg = self.isValid()
if not isValid:
logging.warning(msg)
self.heal()
###
# Groups
#
# inlet, outlet, simetry(N), strips(optional), wall
##
self.groups = []
groupAll = self.createGroupAll(self.shape)
self.groups.append(self.createGroup(self.shape, inletface, "inlet", [grains], 1 / scale))
self.groups.append(self.createGroup(self.shape, outletface, "outlet", [grains], 1 / scale))
for n, face in enumerate(symetryface):
self.groups.append(self.createGroup(self.shape, face, f"symetry{ n }", [grains], 1 / scale))
if self.filletsEnabled:
shapeShell = self.geo.ExtractShapes(self.shape, self.geo.ShapeType["SHELL"], True)[0]
self.groups.append(self.createGroup(self.shape, shapeShell, "strips", self.groups + [grainsOrigin]))
self.groups.append(self.geo.CutListOfGroups([groupAll], self.groups, theName = "wall"))
class _Mesh(Mesh):
def __init__(self, shape):
Mesh.__init__(self, shape)
algo2d = self.mesh.Triangle(algo = self.smeshBuilder.NETGEN_1D2D)
hypo2d = algo2d.Parameters()
hypo2d.SetMaxSize(0.1)
hypo2d.SetMinSize(0.001)
hypo2d.SetFineness(5)
hypo2d.SetGrowthRate(0.3)
hypo2d.SetNbSegPerEdge(2)
hypo2d.SetNbSegPerRadius(3)
hypo2d.SetChordalErrorEnabled(True)
hypo2d.SetChordalError(0.05)
hypo2d.SetOptimize(True)
hypo2d.SetUseSurfaceCurvature(True)
algo3d = self.mesh.Tetrahedron(algo = self.smeshBuilder.NETGEN_3D)
#hypo3d = algo3d.Parameters()
#faces = [ group for group in self.geom.groups if group.GetName() in ["inlet", "outlet"] ]
#hypo3dVL = algo3d.ViscousLayers(...)
def build(self):
out, err, returncode = self.mesh.compute()
if not returncode:
self.mesh.removePyramids()
self.mesh.createGroups()
class _OnePhaseFlow(object):
def __init__(self):
controlDict = ControlDict()
controlDict.update(
startFrom = "latestTime",
endTime = 5000,
writeInterval = 100,
runTimeModifiable = "true"
)
fvSchemes = FvSchemes()
fvSolution = FvSolution()
fvSolution["solvers"]["U"].update(
nSweeps = 2,
tolerance = 1e-08
)
fvSolution["solvers"]["Phi"] = dict(
solver = "GAMG",
smoother = "DIC",
cacheAgglomeration = "yes",
agglomerator = "faceAreaPair",
nCellsInCoarsestLevel = 10,
mergeLevels = 1,
tolerance = 1e-06,
relTol = 0.01
)
fvSolution["potentialFlow"] = dict(
nNonOrthogonalCorrectors = 20,
PhiRefCell = 0,
PhiRefPoint = 0,
PhiRefValue = 0,
Phi = 0
)
fvSolution["cache"] = { "grad(U)": None }
fvSolution["SIMPLE"].update(
nNonOrthogonalCorrectors = 10,
residualControl = dict(
p = 1e-05,
U = 1e-05
)
)
fvSolution["relaxationFactors"]["equations"]["U"] = 0.5
transportProperties = TransportProperties()
transportProperties.update(
nu = 1e-06
)
turbulenceProperties = TurbulenceProperties()
turbulenceProperties.content = dict(
simulationType = "laminar"
)
createPatchDict = CreatePatchDict()
createPatchDict["patches"] = []
for patch in ["inlet", "outlet", "wall", "strips", *[ f"symetry{ n }" for n in range(6) ]]:
newPatch = dict(
name = patch,
patchInfo = dict(
type = "patch",
inGroups = [patch]
),
constructFrom = "patches",
patches = [ f"smesh_{ patch }" ]
)
match patch:
case "wall" | "strips":
newPatch["patchInfo"].update(
type = "wall",
inGroups = [ "wall" ]
)
case patch if patch.find("symetry") == 0:
newPatch["patchInfo"]["inGroups"] = [ "symetryPlane" ]
createPatchDict["patches"].append(newPatch)
boundaries = [ "inlet", "outlet", "symetryPlane", "wall", "strips" ]
p = P()
p["boundaryField"] = {}
u = U()
u["boundaryField"] = {}
# ISSUE: add proxy from geometry direction to outlet boundaryField.
for boundary in boundaries:
match boundary:
case "inlet":
p["boundaryField"][boundary] = dict(
type = "fixedValue",
value = "uniform 1e-3"
)
u["boundaryField"][boundary] = dict(
type = "fixedValue",
value = "uniform (0 0 -6e-5)" # * direction
)
case "outlet":
p["boundaryField"][boundary] = dict(
type = "fixedValue",
value = "uniform 0"
)
u["boundaryField"][boundary] = dict(
type = "zeroGradient",
)
case _:
p["boundaryField"][boundary] = dict(
type = "zeroGradient"
)
u["boundaryField"][boundary] = dict(
type = "fixedValue",
value = "uniform (0 0 0)"
)
self.solution = FoamCase([
controlDict, fvSchemes, fvSolution, createPatchDict,
transportProperties, turbulenceProperties,
p, u
])
def build(self):
self.solution.write()
openfoam.ideasUnvToFoam("mesh.unv")
openfoam.createPatch()
openfoam.checkMesh()
openfoam.transformPoints((1e-5, 1e-5, 1e-5))
openfoam.renumberMesh()
openfoam.potentialFoam()
self.solution.read()
self.solution.U["boundaryField"]["outlet"] = dict(
type = "pressureInletVelocity",
value = "uniform (0 0 0)" # * direction
)
self.solution.write()
openfoam.simpleFoam()
class Simple(object):
geometry = _Geometry
mesh = _Mesh
onephaseflow = _OnePhaseFlow

16
anisotropy/shaping/bodyCentered.py
View File

@ -70,7 +70,7 @@ class BodyCentered(Periodic):
self.lattice = lattice.Scale(zero, self.minimize) self.lattice = lattice.Scale(zero, self.minimize)
# Inlet face # Inlet face
if self.direction in [[1, 0, 0], [0, 0, 1]]: if (self.direction == numpy.array([1., 0., 0.])).prod():
length = 2 * self.r0 length = 2 * self.r0
width = self.L / 2 width = self.L / 2
diag = self.L * sqrt(2) diag = self.L * sqrt(2)
@ -80,7 +80,6 @@ class BodyCentered(Periodic):
yw = xl yw = xl
zh = height zh = height
if self.direction == [1, 0, 0]:
vertices = numpy.array([ vertices = numpy.array([
(xl, 0, 0), (xl, 0, 0),
(0, yw, 0), (0, yw, 0),
@ -89,7 +88,16 @@ class BodyCentered(Periodic):
]) ])
extr = diag extr = diag
elif self.direction == [0, 0, 1]: 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([ vertices = numpy.array([
(0, yw, 0), (0, yw, 0),
(xl, 0, 0), (xl, 0, 0),
@ -98,7 +106,7 @@ class BodyCentered(Periodic):
]) ])
extr = height extr = height
elif self.direction == [1, 1, 1]: elif (self.direction == numpy.array([1., 1., 1.])).prod():
length = 2 * self.r0 length = 2 * self.r0
width = self.L / 2 width = self.L / 2
diag = self.L * sqrt(2) diag = self.L * sqrt(2)

17
anisotropy/shaping/faceCentered.py
View File

@ -77,7 +77,7 @@ class FaceCentered(Periodic):
self.lattice = lattice.Scale(zero, self.minimize) self.lattice = lattice.Scale(zero, self.minimize)
# Inlet face # Inlet face
if self.direction in [[1, 0, 0], [0, 0, 1]]: if (self.direction == numpy.array([1., 0., 0.])).prod():
length = 2 * self.r0 length = 2 * self.r0
width = self.L / 2 width = self.L / 2
diag = self.L * sqrt(3) diag = self.L * sqrt(3)
@ -87,7 +87,6 @@ class FaceCentered(Periodic):
yw = xl yw = xl
zh = width zh = width
if self.direction == [1, 0, 0]:
vertices = numpy.array([ vertices = numpy.array([
(0, 0, -zh), (0, 0, -zh),
(-xl, yw, -zh), (-xl, yw, -zh),
@ -96,7 +95,16 @@ class FaceCentered(Periodic):
]) ])
extr = length extr = length
elif self.direction == [0, 0, 1]: 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([ vertices = numpy.array([
(0, 0, -zh), (0, 0, -zh),
(xl, yw, -zh), (xl, yw, -zh),
@ -105,7 +113,7 @@ class FaceCentered(Periodic):
]) ])
extr = 2 * width extr = 2 * width
elif self.direction == [1, 1, 1]: elif (self.direction == numpy.array([1., 1., 1.])).prod():
length = 2 * self.r0 length = 2 * self.r0
width = self.L / 2 width = self.L / 2
diag = self.L * sqrt(3) diag = self.L * sqrt(3)
@ -139,6 +147,7 @@ class FaceCentered(Periodic):
# Boundaries # Boundaries
symetryId = 0 symetryId = 0
# ISSUE: inlet and outlet faces have the same name and normal vector after extrusion
for face in self.cell.faces: for face in self.cell.faces:
fNorm = self.normal(face) fNorm = self.normal(face)
fAngle = self.angle(vecFlow, fNorm) fAngle = self.angle(vecFlow, fNorm)

27
anisotropy/shaping/occExtended.py Normal file
View File

@ -0,0 +1,27 @@
# -*- 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])

48
anisotropy/shaping/simple.py
View File

@ -6,8 +6,26 @@ from netgen.occ import *
import numpy import numpy
from numpy import pi, sqrt from numpy import pi, sqrt
from .occExtended import *
from . import Periodic from . import Periodic
def reconstruct(solid):
solidNew = []
for face in solid.faces:
vertices = numpy.array([ v.p.pos() for v in face.vertices ])
# TODO: correct following vertices
vertices = numpy.unique(vertices, axis = 0)
circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ])
faceNew = Face(circuit)
faceNew.name = face.name
solidNew.append(faceNew)
return numpy.array(solidNew).sum()
class Simple(Periodic): class Simple(Periodic):
def __init__( def __init__(
self, self,
@ -64,7 +82,7 @@ class Simple(Periodic):
self.lattice = lattice.Scale(zero, self.minimize) self.lattice = lattice.Scale(zero, self.minimize)
# Inlet face # Inlet face
if self.direction in [[1, 0, 0], [0, 0, 1]]: if (self.direction == numpy.array([1., 0., 0.])).prod():
length = self.L * numpy.sqrt(2) length = self.L * numpy.sqrt(2)
width = self.L * numpy.sqrt(2) width = self.L * numpy.sqrt(2)
height = self.L height = self.L
@ -73,8 +91,6 @@ class Simple(Periodic):
yw = xl yw = xl
zh = height zh = height
# TODO: correct compasion for arrays
if self.direction == [1, 0, 0]:
vertices = numpy.array([ vertices = numpy.array([
(xl, 0, 0), (xl, 0, 0),
(0, yw, 0), (0, yw, 0),
@ -83,7 +99,15 @@ class Simple(Periodic):
]) ])
extr = width extr = width
elif self.direction == [0, 0, 1]: 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([ vertices = numpy.array([
(0, yw, 0), (0, yw, 0),
(xl, 0, 0), (xl, 0, 0),
@ -92,7 +116,7 @@ class Simple(Periodic):
]) ])
extr = height extr = height
elif self.direction == [1, 1, 1]: elif (self.direction == numpy.array([1., 1., 1.])).prod():
length = self.L * numpy.sqrt(2) length = self.L * numpy.sqrt(2)
width = self.L * numpy.sqrt(2) width = self.L * numpy.sqrt(2)
height = self.L height = self.L
@ -114,22 +138,32 @@ class Simple(Periodic):
else: else:
raise Exception(f"Direction { self.direction } is not implemented") raise Exception(f"Direction { self.direction } is not implemented")
#
# Cell # Cell
circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ]) circuit = Wire([ Segment(Pnt(*v1), Pnt(*v2)) for v1, v2 in zip(vertices, numpy.roll(vertices, -1, axis = 0)) ])
inletface = Face(circuit) inletface = Face(circuit)
inletface.name = "inlet" inletface.name = "inlet"
vecFlow = self.normal(inletface) vecFlow = self.normal(inletface)
# ISSUE: netgen.occ.Face.Extrude: the opposite face has the same name and normal vector as an initial face.
self.cell = inletface.Extrude(extr) self.cell = inletface.Extrude(extr)
self.cell = reconstruct(self.cell)
print([ f.name for f in self.cell.faces ])
# Boundaries # Boundaries
symetryId = 0 symetryId = 0
for face in self.cell.faces: for face in self.cell.faces:
fNorm = self.normal(face) fNorm = self.normal(face)
fAngle = self.angle(vecFlow, fNorm) fAngle = self.angle(vecFlow, fNorm)
print((face.center.pos() == inletface.center.pos()))
print(fAngle)
if fAngle == 0:
if (face.center.pos() == inletface.center.pos()).prod():
face.name = "inlet"
if fAngle == 0 and not face.center == inletface.center: else:
face.name = "outlet" face.name = "outlet"
else: else:

6
anisotropy/solving/onephase.py
View File

@ -11,8 +11,8 @@ from anisotropy.openfoam.presets import (
from anisotropy.openfoam.foamcase import FoamCase from anisotropy.openfoam.foamcase import FoamCase
class OnePhaseFlow(FoamCase): class OnePhaseFlow(FoamCase):
def __init__(self): def __init__(self, path: str = None):
FoamCase.__init__(self) FoamCase.__init__(self, path = path)
controlDict = ControlDict() controlDict = ControlDict()
controlDict.update( controlDict.update(
@ -127,7 +127,7 @@ class OnePhaseFlow(FoamCase):
self.read() self.read()
self.solution.U["boundaryField"]["outlet"] = dict( self.U["boundaryField"]["outlet"] = dict(
type = "pressureInletVelocity", type = "pressureInletVelocity",
value = "uniform (0 0 0)" # * direction value = "uniform (0 0 0)" # * direction
) )

65
tests/test_core.py Normal file
View File

@ -0,0 +1,65 @@
import os
from os import path
import unittest
unittest.TestLoader.sortTestMethodsUsing = None
class TestCore(unittest.TestCase):
def setUp(self):
try:
import netgen
NETGEN_MODULE = True
except ImportError:
NETGEN_MODULE = False
if not NETGEN_MODULE:
self.skipTest("Missing Netgen.")
else:
from anisotropy import core
self.core = core
self.outputPath = os.path.join(os.path.abspath("."), "tests/test_core_output")
os.makedirs(self.outputPath, exist_ok = True)
def test_config(self):
from copy import deepcopy
config = self.core.DefaultConfig()
contentOld = deepcopy(config.content)
filepath = os.path.join(self.outputPath, "test_config.toml")
config.dump(filepath)
config = self.core.Config()
config.load(filepath)
self.assertEqual(contentOld, config.content)
def test_runner(self):
os.chdir(self.outputPath)
pathOld = os.path.abspath(".")
config = self.core.DefaultConfig()
config.expand()
config.cases = [ config.cases[0] ]
runner = self.core.UltimateRunner(config = config, exec_id = True)
runner.computeShape()
self.assertTrue(path.isfile(path.join(runner.casepath(), "shape.step")))
runner.computeMesh()
self.assertTrue(path.isfile(path.join(runner.casepath(), "mesh.mesh")))
runner.computeFlow()
#self.assertTrue(path.isfile(path.join(runner.casepath(), "mesh.mesh")))
os.chdir(pathOld)
def tearDown(self):
pass
if __name__ == "__main__":
unittest.main()

39
tests/test_database.py Normal file
View File

@ -0,0 +1,39 @@
import os
import unittest
unittest.TestLoader.sortTestMethodsUsing = None
class TestDatabase(unittest.TestCase):
def setUp(self):
from anisotropy import database
self.database = database
self.outputPath = os.path.join(os.path.abspath("."), "tests/test_database_output")
os.makedirs(self.outputPath, exist_ok = True)
def test_setup(self):
filepath = os.path.join(self.outputPath, "test_database.db")
tables = [
self.database.Execution,
self.database.Physics,
self.database.Shape,
self.database.Mesh,
self.database.Flow
]
db = self.database.Database(filepath)
db.setup()
self.assertTrue(
os.path.exists(filepath) and os.path.isfile(filepath),
"database wasn't created"
)
for table in tables:
self.assertTrue(table.table_exists())
def tearDown(self):
pass
if __name__ == "__main__":
unittest.main()

4
tests/test_shaping.py
View File

@ -7,10 +7,10 @@ class TestShaping(unittest.TestCase):
def setUp(self): def setUp(self):
try: try:
import netgen import netgen
NETGEN_MODULE = False NETGEN_MODULE = True
except ImportError: except ImportError:
NETGEN_MODULE = True NETGEN_MODULE = False
if not NETGEN_MODULE: if not NETGEN_MODULE:
self.skipTest("Missing Netgen.") self.skipTest("Missing Netgen.")