Compare commits

...

81 Commits

Author SHA1 Message Date
f3191aeafb
stage changes 2024-03-14 16:04:24 +05:00
f06eb9aaad
new: test new gui (dpg) 2022-04-21 19:19:52 +05:00
e0c6686d2b
new: data resources 2022-04-21 19:19:33 +05:00
c6136b6e4b
mod: minor changes 2022-04-21 19:19:14 +05:00
e26eaf01eb
Mod: update data analyze figures 2022-03-05 01:08:18 +05:00
da63856466 Mod: data files 2022-03-04 17:19:56 +05:00
fc1b8432cb
Mod: core: improve pipeline
Fix: openfoam: wrong path expressions
2022-02-18 00:28:05 +05:00
c5d2bf1783
Mod: real volume 2022-02-16 12:27:36 +05:00
fea53d0db7
Mod: refactoring gui 2022-02-12 13:58:24 +05:00
13e02b57a4
Fix: shape primitives
Mod: default config
2022-02-12 01:10:15 +05:00
c701f0452b
Fix: digging for errors 2022-02-11 00:01:36 +05:00
fb52e4c6c8 Mod: refactoring ultimate runner 2022-02-01 16:06:52 +05:00
34313205d6
Mod: clean up
Mod: cli is standalone submodule now
Mod: gui layouts collected together
Mod: pyproject cli entry
2022-01-29 17:41:36 +05:00
12ee8ab0b0
Mod: improved file and case usage
Mod: openfoam submodule documentation
2022-01-29 16:44:46 +05:00
33ff54bd3b
Mod: formatting openfoam submodule 2022-01-28 21:16:35 +05:00
6bfbb83b77
Mod: openfoam conda requirement 2022-01-28 12:41:09 +05:00
a4e9a8f8dc Mod: database is now stable.
Improved convenience and documention.
2022-01-24 17:50:57 +05:00
00f7279c6d
Mod: shaping is stable now
Mod: shaping improved documentation
Remove: not necessary modules or moved
2022-01-22 22:08:47 +05:00
ac9e938a50
Mod: improved meshing submodule
New: netgen neutral format conversion
Mod: updated desc
2022-01-22 00:07:22 +05:00
7dce03e97a
Mod: improved mesh and shape classes
Mod: exception handling
2022-01-19 23:02:10 +05:00
d724b45e85 Mod: meshio object 2022-01-14 14:14:56 +05:00
c7111dc9cc - 2022-01-14 09:26:18 +05:00
8fa62bba38
Mod: mesh and shape io
Mod: mesh patches
2022-01-14 00:12:08 +05:00
e501b8b52b
Mod: improved plot for each data
Mod: cwd for gui settings
Mod: smooth stages
2022-01-12 23:17:42 +05:00
479e7d666a
New: mesh visualization
New: data visualization
Mod: python deps
2022-01-11 22:57:19 +05:00
3eb1029cb5 Mod: extra gui dependencies
Fix: parallel database usage
Fix: reduced memory usage
Fix: env paths in gui
2022-01-10 15:38:55 +05:00
b41c730e48
Mod: update packages
Mov: documentation
2021-12-31 20:01:59 +05:00
29a9a7fa7f
New: gui
Mod: global env
2021-12-31 01:18:39 +05:00
59e44dd08d
Fix: postProcess incorrect log path 2021-12-25 13:36:27 +05:00
11b612a99f
Fix: missed viscosity value
Fix: postProcess incorrect log path
Mod: flake8 poetry deps
2021-12-25 13:35:23 +05:00
9d4aac4652
Mod: track new model variables
Mod: return mesh elements
Mod: update netgen conda channel
Fix: some mistakes
2021-12-25 01:02:28 +05:00
fba3e53bb6
Mod: database new fields 2021-12-22 22:50:15 +05:00
49368cc681
Mod: improved openfoam interface
Mod: improved runners exception handlers
Delete: some useless trash
New: post process pipeline block
2021-12-21 20:36:01 +05:00
8e769ec1ce
Mod: working database improvements
Fix: FoamRunner exceptions
Mod: configs for subprocess runners
2021-12-20 20:52:46 +05:00
6601525d22 Mod: database improvements, new methods
Mod: conda netgen version changed to master
2021-12-20 16:16:10 +05:00
c166be7801
Mod: openfoam runner and presets 2021-12-15 22:21:04 +05:00
49a7390473
Mod: tests and conda package channel for netgen 2021-12-14 17:50:00 +05:00
f94333ea0b
Mod: deps and runner 2021-12-13 04:18:42 +05:00
ace2cd0247 Mod: working mesh stage 2021-12-10 19:33:54 +05:00
7454714d70 Mod: working shape stage 2021-12-09 17:24:44 +05:00
c35d5cfe3c Mod: db improved 2021-12-09 13:59:32 +05:00
64a5cf1a6d Mod: improved database
Mod: new cli
Mod: improved ultimate runner
Mod: improved multiprocessing
2021-12-08 19:04:38 +05:00
7d0fe5a4a9 Mod: minor changes 2021-12-07 15:57:27 +05:00
c2ef046f74 Mod: minor changes 2021-12-06 19:19:34 +05:00
9e98df1ffd Mod: graph 2021-12-03 19:43:33 +05:00
17934fe17a Mod: conda environment 2021-12-02 18:14:20 +05:00
604a877b6e
Fix: missed inlet/outlet boundary 2021-12-01 16:58:06 +05:00
8df065d5aa
New: conda environment 2021-11-25 22:08:51 +05:00
e977c39496
Fix: faces after extrusion had shared object 2021-11-23 00:46:56 +05:00
f0ce49a0d8
Mod: trying to fix reconstruct method 2021-11-22 15:59:10 +05:00
71d25cb510
New: oce, OpenFOAM and ThirdParty externals 2021-11-21 16:20:25 +05:00
f64864a29a
New: netgen external 2021-11-21 16:09:24 +05:00
c22aa7f3f2
Mod: trying to fix issue 2021-11-21 15:00:04 +05:00
8f4c664a74
Mod: new tests
Mod: working core runner
2021-11-19 23:35:12 +05:00
c56504f120
Mod: merge cae to pipeline 2021-11-18 22:00:30 +05:00
88f07abf4a
Mod: skip test if thirdparty not found 2021-11-18 12:34:43 +05:00
b08a798c08
Mod: up pytest version 2021-11-18 11:52:35 +05:00
3cf2391399 Updated config.yml 2021-11-18 11:17:06 +05:00
2c27c25389
New: circleci 2021-11-18 11:12:59 +05:00
8920e89928
New: poetry
Mod: shaping test
2021-11-18 00:48:17 +05:00
6d90fa1c33
Mod: occ geometry 2021-11-04 18:26:13 +05:00
06dcb26391
New: netgen to foam mesh conversion 2021-11-02 22:01:43 +05:00
085ec5a0e2
Mod: - 2021-10-26 13:47:10 +05:00
120c2f260c
New: config class 2021-10-25 17:19:15 +05:00
0859b6bdf5
Mod: database models 2021-10-23 17:36:17 +05:00
493094749d
Mod: working flow process for sample 2021-10-22 21:26:36 +05:00
b2343aabb4
Mod: simple, foamcase for approximation 2021-10-22 14:02:50 +05:00
3071a81d1b
Mod: foamfile presets 2021-10-21 16:04:19 +05:00
b45954d153
Mod: improved mesh class
Mod: mesh subclass for each sample
New: foamfile, reader and writer
2021-10-20 21:48:59 +05:00
77647b3370
Mod: filletScale parameter 2021-10-18 17:10:30 +05:00
c57ac1154a
Fix: fillet radius on all structures
Move: StructureGeometry class to geometry.py
2021-10-13 21:47:46 +05:00
1d36353628
Mod: merge 2021-10-12 22:33:51 +05:00
d43c4bd012
Mod: - 2021-10-12 22:31:57 +05:00
0b9a6d3f5c
Fix: faceCentered fillets radius 2021-10-11 16:30:43 +05:00
3b896fb3fc
Mod: - 2021-10-10 21:46:54 +05:00
73e6c98675
Mod: geometry export function 2021-10-07 22:16:58 +05:00
1374ce0000
Mod: improved geometry construction by adding flexibility
New: StructureGeometry class
Fix: fixed invalid shapes (external edges on shapes, etc)
New: issues list
2021-10-06 22:15:05 +05:00
eda9a56a28
New: extra requirements + jupyter 2021-10-06 01:02:41 +05:00
8c2953b3ec
New: failed analytics 2021-10-06 00:58:37 +05:00
2ea14a5770
New: playground notebook 2021-10-04 16:08:29 +05:00
8ed50bc9a3
New: calc tests for comparison b/n mesh w(/o) prismatic layer 2021-10-04 12:16:01 +05:00
127 changed files with 9122 additions and 4249 deletions

28
.circleci/config.yml Normal file
View File

@ -0,0 +1,28 @@
version: 2
workflows:
version: 2
test:
jobs:
- test-3.10
jobs:
test-3.10:
docker:
- image: cimg/python:3.10
steps:
- checkout
- run:
name: install packages
command: |
pip install poetry
poetry config virtualenvs.in-project true
poetry install
- run:
name: run tests
command: |
mkdir test-reports
poetry run pytest --cov=anisotropy --json-report --json-report-file test-reports/pytest.json -v tests
- store_artifacts:
path: test-reports

9
.gitignore vendored
View File

@ -4,6 +4,7 @@ logs/
storage/
temp/
env/
*output/
*.gz
*.xz
*.fls
@ -18,3 +19,11 @@ env/
*.toc
*.out
*.egg-info
.ipynb_checkpoints/
.coverage
.venv
.pytest_cache
dist/
doc/anisotropy/*
!doc/anisotropy/.gitkeep
/.idea

15
.vscode/launch.json vendored Normal file
View File

@ -0,0 +1,15 @@
{
"version": "0.2.0",
"configurations": [
{
"name": "anisotropy compute",
"type": "python",
"request": "launch",
"program": "${workspaceFolder}/anisotropy/core/cli.py",
"args": ["compute", "-v"],
"console": "integratedTerminal",
"cwd": "/tmp/anisotropy",
"preLaunchTask": "prepare"
}
]
}

12
.vscode/settings.json vendored Normal file
View File

@ -0,0 +1,12 @@
{
"python.pythonPath": "${env:HOME}/.local/share/mambaforge/envs/anisotropy/bin/python",
"python.testing.pytestArgs": [
"tests"
],
"python.testing.unittestEnabled": false,
"python.testing.pytestEnabled": true,
"python.linting.enabled": true,
"python.linting.pylintEnabled": false,
"python.linting.flake8Enabled": true,
"python.linting.flake8Args": ["--ignore=E402,E251,E501,E201,E202,W293,W291,W504,E203", "--verbose"]
}

9
.vscode/tasks.json vendored Normal file
View File

@ -0,0 +1,9 @@
{
"version": "2.0.0",
"tasks": [{
"label": "prepare",
"command": "mkdir",
"args": ["-p", "/tmp/anisotropy"],
"type": "shell",
}]
}

2
ISSUES.rst Normal file
View File

@ -0,0 +1,2 @@
List of known issues
====================

View File

@ -1,3 +1,6 @@
.. image:: https://circleci.com/gh/L-Nafaryus/anisotropy/tree/devel.svg?style=shield&circle-token=423bc964a997ded671ebd4ceacc25f9967acdffa
:target: https://circleci.com/gh/L-Nafaryus/anisotropy/tree/devel
anisotropy
==========
@ -31,7 +34,14 @@ Dependencies
.. Installation
.. include:: INSTALL.rst
.. code-block:: bash
$ git clone https://github.com/L-Nafaryus/anisotropy.git
$ cd anisotropy
$ conda env create
$ conda activate anisotropy
$ poetry install
$ anisotropy --help
Getting Started
===============

View File

@ -1,40 +1,42 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
"""anisotropy
*anisotropy* is a ``Python`` package that is the result of science-research work
on the anisotropy of permeability in the periodic porous media.
A project uses own wrappers around external applications
for constructing a shapes and meshes (``Salome``) and computing a flow (``OpenFOAM``).
"""
__license__ = "GPL3"
__version__ = "1.2.0"
__author__ = __maintainer__ = "George Kusayko"
__email__ = "gkusayko@gmail.com"
__repository__ = "https://github.com/L-Nafaryus/anisotropy"
###
# Environment
##
import os
from os import path, environ
env = dict(
ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
)
env.update(
BUILD = os.path.join(env["ROOT"], "build"),
LOG = os.path.join(env["ROOT"], "logs"),
CONFIG = os.path.join(env["ROOT"], "anisotropy/config/default.toml"),
DOCS = os.path.join(env["ROOT"], "docs")
)
env.update(
logger_name = "anisotropy",
db_name = "anisotropy",
db_path = env["BUILD"],
salome_timeout = 15 * 60,
openfoam_template = os.path.join(env["ROOT"], "anisotropy/openfoam/template")
)
del os
PROJECT = path.abspath(path.dirname(__file__))
TMP = "/tmp/anisotropy"
env = {
"PROJECT": path.abspath(path.dirname(__file__)),
"TMP": TMP,
"CWD": TMP,
"BUILD_DIR": "build",
"CONF_FILE": "anisotropy.toml",
"LOG_FILE": "anisotropy.log",
"DB_FILE": "anisotropy.db"
}
def loadEnv():
prefix = "AP_"
for k, v in env.items():
environ[f"{ prefix }{ k }"] = v
del path, PROJECT, TMP

View File

@ -0,0 +1,15 @@
# -*- coding: utf-8 -*-
from . import utils
from .cli import anisotropy_cli
__all__ = [
"utils",
"anisotropy_cli"
]
if __name__ == "__main__":
anisotropy_cli()

181
anisotropy/cli/cli.py Normal file
View File

@ -0,0 +1,181 @@
# -*- coding: utf-8 -*-
import click
import pathlib
import os
import logging
from . import utils
@click.group()
@click.version_option()
def anisotropy_cli():
pass
@anisotropy_cli.command(
help = "Initialize new anisotropy project."
)
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
@click.option(
"-v", "--verbose", "verbose",
count = True,
help = "Increase verbose level"
)
def init(path, verbose):
from anisotropy.core import config as core_config
from anisotropy.core import utils as core_utils
core_utils.setupLogger(utils.verbose_level(verbose))
logger = logging.getLogger(__name__)
config = core_config.default_config()
filepath = os.path.abspath(os.path.join(path, "anisotropy.toml"))
logger.info(f"Saving file at { filepath }")
config.dump(filepath)
@anisotropy_cli.command()
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
@click.option(
"-C", "--conf", "configFile"
)
@click.option(
"-N", "--nprocs", "nprocs",
type = click.INT,
# default = 1,
help = "Count of parallel processes"
)
@click.option(
"-s", "--stage", "stage",
type = click.Choice(["all", "shape", "mesh", "flow", "postProcess"]),
# default = "all",
help = "Current computation stage"
)
@click.option(
"-f", "--force", "overwrite",
is_flag = True,
default = None,
help = "Overwrite existing entries"
)
@click.option(
"-p", "--params", "params",
metavar = "key=value",
multiple = True,
cls = utils.KeyValueOption,
help = "Overwrite existing parameter (except control variables)"
)
@click.option(
"-v", "--verbose", "verbose",
count = True,
help = "Increase verbose level"
)
@click.option(
"--exec-id", "execution"
)
@click.option(
"--pid", "pid",
help = "Specify pid file path"
)
@click.option(
"--logfile", "logfile",
help = "Specify log file path"
)
def compute(path, configFile, nprocs, stage, overwrite, params, verbose, execution, pid, logfile):
import anisotropy
from anisotropy.core import UltimateRunner
from anisotropy.core import config as core_config
from anisotropy.core import utils as core_utils
anisotropy.loadEnv()
if path:
path = pathlib.Path(path).resolve()
os.makedirs(path, exist_ok = True)
os.chdir(path)
os.environ["AP_CWD"] = str(path)
core_utils.setupLogger(utils.verbose_level(verbose), logfile)
logger = logging.getLogger(__name__)
config = core_config.default_config()
if configFile:
configFile = pathlib.Path(configFile).resolve()
logger.info(f"Loading file from { configFile }")
try:
config.load(configFile)
except FileNotFoundError:
config.dump(configFile)
else:
logger.info("Using default configuration")
args = {
"nprocs": nprocs,
"stage": stage,
"overwrite": overwrite
}
for k, v in args.items():
if v is not None:
config.update(**{ k: v })
if pid:
pid = pathlib.Path(pid).resolve()
with open(pid, "w") as io:
io.write(str(os.getpid()))
runner = UltimateRunner(config = config, exec_id = execution)
runner.prepare_queue()
runner.start()
if pid:
os.remove(pid)
logger.info("Computation done. Exiting ...")
@anisotropy_cli.command()
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
@click.option(
"-v", "--verbose", "verbose",
count = True,
help = "Increase verbose level"
)
def gui(path, verbose):
import anisotropy
from anisotropy.core import utils as core_utils
from anisotropy import gui
anisotropy.loadEnv()
path = pathlib.Path(path).resolve()
path.mkdir(parents = True, exist_ok = True)
os.chdir(path)
os.environ["ANISOTROPY_CWD"] = str(path)
core_utils.setupLogger(utils.verbose_level(verbose))
# logger = logging.getLogger(__name__)
gui.run(debug = True)

74
anisotropy/cli/utils.py Normal file
View File

@ -0,0 +1,74 @@
# -*- coding: utf-8 -*-
import click
import ast
import logging
class LiteralOption(click.Option):
def type_cast_value(self, ctx, value):
try:
return ast.literal_eval(value)
except Exception:
raise click.BadParameter(f"{ value } (Type error)")
class KeyValueOption(click.Option):
def _convert(self, ctx, value):
if not value:
return {}
if value.find("=") == -1:
raise click.BadParameter(f"{ value } (Missed '=')")
params = value.split("=")
if not len(params) == 2:
raise click.BadParameter(f"{ value } (Syntax error)")
key, val = params[0].strip(), params[1].strip()
if val[0].isalpha():
val = f"'{ val }'"
try:
return { key: ast.literal_eval(val) }
except Exception:
raise click.BadParameter(f"{ value } (Type error)")
def type_cast_value(self, ctx, value):
if isinstance(value, list):
return [ self._convert(ctx, val) for val in value ]
else:
return self._convert(ctx, value)
class CliListOption(click.Option):
def _convert(self, ctx, value):
if not value:
return []
output = [ val for val in value.split(",") ]
if "" in output:
raise click.BadParameter(f"{ value } (Trying to pass empty item)")
return output
def type_cast_value(self, ctx, value):
if isinstance(value, list):
return [ self._convert(ctx, val) for val in value ]
else:
return self._convert(ctx, value)
def verbose_level(level: int):
return {
0: logging.ERROR,
1: logging.INFO,
2: logging.DEBUG
}.get(level, logging.ERROR)

View File

@ -1,3 +0,0 @@
# Example of anisotropy bashrc entries (modify it)
export PATH="${HOME}/programs/salome/SALOME-9.7.0-MPI:${PATH}"
source "${HOME}/programs/OpenFOAM/OpenFOAM-v2012/etc/bashrc"

View File

@ -1,281 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
[logger]
name = "anisotropy"
format = "[ %(levelname)s ] %(message)s"
[base]
simple = true
bodyCentered = true
faceCentered = true
[[structures]]
[structures.structure]
type = "simple"
# auto # from theta: list # theta: float
theta = [0.01, 0.28, 0.01] # [min, max, step]
# auto # from directions:list # direction: list
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
# r0 =
# L =
# radius =
filletsEnabled = true
# auto # fillets: float
[structures.mesh]
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.5
nbSegPerEdge = 2
nbSegPerRadius = 1
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
viscousLayers = true
thickness = [0.01, 0.005] # [min, max] # step is controlled by theta count
numberOfLayers = 1
stretchFactor = 1
isFacesToIgnore = true
facesToIgnore = ["inlet", "outlet"]
# auto # faces: list
extrusionMethod = "SURF_OFFSET_SMOOTH"
[[structures.submesh]]
name = "strips"
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.2
nbSegPerEdge = 2
nbSegPerRadius = 3
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
[structures.flowapproximation]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
# multiplication velocity value with flow direction vector
velocity.boundaryField.inlet = { type = "fixedValue", value = 6e-5 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }
[structures.flow]
scale = [ 1e-5, 1e-5, 1e-5 ]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
velocity.boundaryField.inlet = { type = "pressureInletVelocity", value = 0.0 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }
[[structures]]
[structures.structure]
type = "bodyCentered"
# auto # from theta: list # theta: float
theta = [0.01, 0.17, 0.01] # [min, max, step]
# auto # from directions:list # direction: list
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
# r0 =
# L =
# radius =
filletsEnabled = true
# auto # fillets: float
[structures.mesh]
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.5
nbSegPerEdge = 2
nbSegPerRadius = 1
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
viscousLayers = true
thickness = [0.005, 0.0005] # [min, max] # step is controlled by theta count
numberOfLayers = 1
stretchFactor = 1
isFacesToIgnore = true
facesToIgnore = ["inlet", "outlet"]
# auto # faces: list
extrusionMethod = "SURF_OFFSET_SMOOTH"
[[structures.submesh]]
name = "strips"
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.2
nbSegPerEdge = 2
nbSegPerRadius = 3
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
[structures.flowapproximation]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
# multiplication velocity value with direction vector
velocity.boundaryField.inlet = { type = "fixedValue", value = 6e-5 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }
[structures.flow]
scale = [ 1e-5, 1e-5, 1e-5 ]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
velocity.boundaryField.inlet = { type = "pressureInletVelocity", value = 0.0 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }
[[structures]]
[structures.structure]
type = "faceCentered"
# auto # from theta: list # theta: float
theta = [0.01, 0.12, 0.01] # [min, max, step]
# auto # from directions:list # direction: list
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
# r0 =
# L =
# radius =
filletsEnabled = true
# auto # fillets: float
[structures.mesh]
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.5
nbSegPerEdge = 2
nbSegPerRadius = 1
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
viscousLayers = true
thickness = [0.001, 0.0005] # [min, max] # step is controlled by theta count
numberOfLayers = 1
stretchFactor = 1
isFacesToIgnore = true
facesToIgnore = ["inlet", "outlet"]
# auto # faces: list
extrusionMethod = "SURF_OFFSET_SMOOTH"
[[structures.submesh]]
name = "strips"
maxSize = 0.5
minSize = 0.05
fineness = 5
growthRate = 0.2
nbSegPerEdge = 2
nbSegPerRadius = 3
chordalErrorEnabled = true
chordalError = 0.25
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
[structures.flowapproximation]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
# multiplication velocity value with direction vector
velocity.boundaryField.inlet = { type = "fixedValue", value = 6e-5 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }
[structures.flow]
scale = [ 1e-5, 1e-5, 1e-5 ]
transportProperties.nu = 1e-6
pressure.boundaryField.inlet = { type = "fixedValue", value = 1e-3 }
pressure.boundaryField.outlet = { type = "fixedValue", value = 0.0 }
velocity.boundaryField.inlet = { type = "pressureInletVelocity", value = 0.0 }
velocity.boundaryField.outlet = { type = "zeroGradient", value = "None" }

View File

@ -1,5 +1,17 @@
# -*- 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 config
from . import postprocess
from .parallel import ParallelRunner
from .runner import UltimateRunner
__all__ = [
"utils",
"config",
"postprocess",
"ParallelRunner",
"UltimateRunner"
]

View File

@ -1,579 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import click
import ast
import os, sys, shutil
import logging
class LiteralOption(click.Option):
def type_cast_value(self, ctx, value):
try:
return ast.literal_eval(value)
except:
raise click.BadParameter(f"{ value } (Type error)")
class KeyValueOption(click.Option):
def _convert(self, ctx, value):
if not value:
return {}
if value.find("=") == -1:
raise click.BadParameter(f"{ value } (Missed '=')")
params = value.split("=")
if not len(params) == 2:
raise click.BadParameter(f"{ value } (Syntax error)")
key, val = params[0].strip(), params[1].strip()
if val[0].isalpha():
val = f"'{ val }'"
try:
return { key: ast.literal_eval(val) }
except:
raise click.BadParameter(f"{ value } (Type error)")
def type_cast_value(self, ctx, value):
if isinstance(value, list):
return [ self._convert(ctx, val) for val in value ]
else:
return self._convert(ctx, value)
class CliListOption(click.Option):
def _convert(self, ctx, value):
if not value:
return []
output = [ val for val in value.split(",") ]
if "" in output:
raise click.BadParameter(f"{ value } (Trying to pass empty item)")
return output
def type_cast_value(self, ctx, value):
if isinstance(value, list):
return [ self._convert(ctx, val) for val in value ]
else:
return self._convert(ctx, value)
def version():
msg = "Missed package anisotropy"
try:
from anisotropy.core.main import Anisotropy
msg = Anisotropy.version()
except ImportError:
pass
return msg
@click.group()
@click.version_option(version = "", message = version())
def anisotropy():
pass
@anisotropy.command(
help = "Initialize new anisotropy project."
)
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
def init(path):
from anisotropy import env
from anisotropy.core.main import Database
if not os.path.exists(path) or not os.path.isdir(path):
click.echo(f"Cannot find directory { path }")
return
wds = [ "build", "logs" ]
for wd in wds:
os.makedirs(os.path.join(path, wd), exist_ok = True)
shutil.copy(env["CONFIG"], os.path.join(path, "anisotropy.toml"), follow_symlinks = True)
db = Database(env["db_name"], path)
db.setup()
click.echo(f"Initialized anisotropy project in { path }")
@anisotropy.command(
help = "Load parameters from configuration file and update database."
)
@click.option(
"-f", "--force", "force",
is_flag = True,
default = False,
help = "Overwrite existing entries"
)
@click.option(
"-p", "--params", "params",
metavar = "key=value",
multiple = True,
cls = KeyValueOption,
help = "Specify control parameters to update (type, direction, theta)"
)
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
def update(force, params, path):
from anisotropy import env
from anisotropy.core.main import Anisotropy, Database
env.update(
LOG = os.path.join(path, "logs"),
BUILD = os.path.join(path, "build"),
CONFIG = os.path.join(path, "anisotropy.toml"),
db_path = path
)
args = dict()
for param in params:
args.update(param)
model = Anisotropy()
database = Database(env["db_name"], env["db_path"])
click.echo("Configuring database ...")
database.setup()
if database.isempty() or update:
paramsAll = model.loadFromScratch(env["CONFIG"])
if args.get("type"):
paramsAll = [ entry for entry in paramsAll if args["type"] == entry["structure"]["type"] ]
if args.get("direction"):
paramsAll = [ entry for entry in paramsAll if args["direction"] == entry["structure"]["direction"] ]
if args.get("theta"):
paramsAll = [ entry for entry in paramsAll if args["theta"] == entry["structure"]["theta"] ]
for entry in paramsAll:
database.update(entry)
click.echo("{} entries was updated.".format(len(paramsAll)))
else:
click.echo("Database was not modified.")
@anisotropy.command(
help = """Compute cases by chain (mesh -> flow)
Control parameters: type, direction, theta (each parameter affects on a queue)
"""
)
@click.option(
"-s", "--stage", "stage",
type = click.Choice(["all", "mesh", "flow", "postProcessing"]),
default = "all",
help = "Current computation stage"
)
@click.option(
"-N", "--nprocs", "nprocs",
type = click.INT,
default = 1,
help = "Count of parallel processes"
)
@click.option(
"-f", "--force", "force",
is_flag = True,
default = False,
help = "Overwrite existing entries"
)
@click.option(
"-p", "--params", "params",
metavar = "key=value",
multiple = True,
cls = KeyValueOption,
help = "Overwrite existing parameter (except control variables)"
)
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
def compute(stage, nprocs, force, params, path):
from anisotropy import env
from anisotropy.core.main import Anisotropy, Database, logger
from anisotropy.core.utils import setupLogger, parallel
env.update(
LOG = os.path.join(path, "logs"),
BUILD = os.path.join(path, "build"),
CONFIG = os.path.join(path, "anisotropy.toml"),
db_path = path
)
setupLogger(logger, logging.INFO, env["LOG"])
args = dict()
for param in params:
args.update(param)
###
logger.info("Writing pid ...")
pidpath = os.path.join(path, "anisotropy.pid")
with open(pidpath, "w") as io:
io.write(str(os.getpid()))
###
# Preparations
##
database = Database(env["db_name"], env["db_path"])
logger.info("Loading database ...")
database.setup()
params = database.loadGeneral(
args.get("type"),
args.get("direction"),
args.get("theta")
)
queueargs = []
for p in params:
s = p["structure"]
queueargs.append((s["type"], s["direction"], s["theta"]))
###
# Wrap function
##
def computeCase(type, direction, theta):
case = Anisotropy()
case.db = database
case.load(type, direction, theta)
case.evalParams()
case.update()
logger.info(f"Case: type = { type }, direction = { direction }, theta = { theta }")
logger.info(f"Stage mode: { stage }")
if stage in ["mesh", "all"]:
case.load(type, direction, theta)
if not case.params["meshresult"]["meshStatus"] == "Done" or force:
logger.info("Current stage: mesh")
out, err, returncode = case.computeMesh(path)
if out: logger.info(out)
if err: logger.error(err)
if returncode:
logger.error("Mesh computation failed. Skipping flow computation ...")
return
else:
logger.info("Mesh exists. Skipping ...")
if stage in ["flow", "all"]:
case.load(type, direction, theta)
if not case.params["flowresult"]["flowStatus"] == "Done" or force:
logger.info("Current stage: flow")
out, err, returncode = case.computeFlow(path)
if out: logger.info(out)
if err: logger.error(err)
if returncode:
logger.error("Flow computation failed.")
return
else:
logger.info("Flow exists. Skipping ...")
if stage in ["postProcessing", "all"]:
case.load(type, direction, theta)
if case.params["meshresult"]["meshStatus"] == "Done":
logger.info("Current stage: mesh postProcessing")
case.porosity()
else:
logger.warning("Cannot compute mesh post processing values.")
if case.params["flowresult"]["flowStatus"] == "Done":
logger.info("Current stage: flow postProcessing")
case.flowRate()
else:
logger.warning("Cannot compute flow post processing values.")
###
# Run
##
if nprocs == 1:
for pos, qarg in enumerate(queueargs):
computeCase(*qarg)
else:
parallel(nprocs, queueargs, computeCase)
if os.path.exists(pidpath):
logger.info("Removing pid ...")
os.remove(pidpath)
logger.info("Computation done.")
@anisotropy.command(
help = "Kill process by pid file"
)
@click.option(
"-P", "--path", "path",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
@click.argument("pidfile")
def kill(path, pidfile):
from anisotropy.salomepl.utils import SalomeManager
try:
with open(os.path.join(path, pidfile), "r") as io:
pid = io.read()
os.kill(int(pid), 9)
except FileNotFoundError:
click.echo(f"Unknown file { pidfile }")
except ProcessLookupError:
click.echo(f"Cannot find process with pid { pid }")
# TODO: killall method kills all salome instances. Not a good way
SalomeManager().killall()
@anisotropy.command(
help = "! Internal command"
)
@click.argument("root")
@click.argument("type")
@click.argument("direction")
@click.argument("theta")
@click.argument("path")
def computemesh(root, type, direction, theta, path):
# ISSUE: can't hide command from help, 'hidden = True' doesn't work
# [Salome Environment]
###
# Args
##
direction = [ float(num) for num in direction[1:-1].split(" ") if num ]
theta = float(theta)
###
# Modules
##
import os, sys
pyversion = "{}.{}".format(3, 9) #(*sys.version_info[:2])
sys.path.extend([
root,
os.path.join(root, "env/lib/python{}/site-packages".format(pyversion)),
os.path.abspath(".local/lib/python{}/site-packages".format(pyversion))
])
from anisotropy import env
from anisotropy.core.main import Anisotropy, Database
import salome
###
model = Anisotropy()
model.db = Database(env["db_name"], path)
model.load(type, direction, theta)
salome.salome_init()
model.genmesh(path)
salome.salome_close()
@anisotropy.command()
@click.option(
"-p", "--params", "params",
metavar = "key=value",
multiple = True,
cls = KeyValueOption,
help = "Select by control parameters (type, direction, theta)"
)
@click.option(
"-P", "--path", "path",
metavar = "PATH",
default = os.getcwd(),
help = "Specify directory to use (instead of cwd)"
)
@click.option(
"--list", "printlist",
is_flag = True,
help = "Print a list of avaliable fields."
)
@click.option(
"--export",
metavar = "PATH",
help = "Export output."
)
@click.option(
"--fields", "fields",
metavar = "f1,f2,...",
multiple = True,
cls = CliListOption,
help = "Select fields to use."
)
@click.argument(
"output",
required = False,
type = click.Choice(["cli", "plot"]),
default = "cli"
)
def show(params, path, printlist, export, fields, output):
from anisotropy import env
from anisotropy.core.database import Database, Structure
from pandas import DataFrame, Series
import matplotlib.pyplot as plt
env.update(
LOG = os.path.join(path, "logs"),
BUILD = os.path.join(path, "build"),
CONFIG = os.path.join(path, "anisotropy.toml"),
db_path = path
)
args = dict()
for param in params:
args.update(param)
###
db = Database(env["db_name"], env["db_path"])
db.setup()
searchargs = []
if args.get("type"):
searchargs.append(Structure.type == args["type"])
if args.get("direction"):
searchargs.append(Structure.direction == str(args["direction"]))
if args.get("theta"):
searchargs.append(Structure.theta == args["theta"])
result = db.search(searchargs)
result.sort(key = lambda src: f"{ src['type'] }{ src['direction'] }{ src['theta'] }")
df = DataFrame(result)
df_keys = [ key for key in df.keys() ]
if printlist:
click.echo("Avaliable fields for query:")
click.echo("\t{}".format("\n\t".join(df_keys)))
return
if not result:
click.echo("Empty result.")
return
tables = []
if fields:
for fieldslist in fields:
for field in fieldslist:
if field not in df_keys:
click.echo(f"Unknown field '{ field }'. Try to use '--list' flag to see all avaliable fields.")
return
tables.append(df[fieldslist])
else:
tables.append(df)
if output == "plot":
fig, ax = plt.subplots(nrows = 1, ncols = 1)
for table in tables:
table.plot(table.keys()[0], table.keys()[1], ax = ax, style = "o")
plt.legend()
plt.grid()
if export:
supported = ["csv", "jpg"]
filepath, ext = os.path.splitext(export)
ext = ext.replace(".", "")
if ext not in supported:
click.echo(f"Unknown extension '{ ext }'.")
return
if ext == "csv":
if len(tables) == 1:
tables[0].to_csv(export, sep = ";")
else:
for n, table in enumerate(tables):
table.to_csv("{}.{}.{}".format(filepath, n, ext), sep = ";")
elif ext == "jpg":
plt.savefig(export)
else:
if output == "cli":
res = "\n\n".join([ table.to_string() for table in tables ])
click.echo(res)
elif output == "plot":
plt.show()
###
# CLI entry
##
if __name__ == "__main__":
#try:
anisotropy()
#except KeyboardInterrupt:
# click.echo("Interrupted!")
#finally:
# from anisotropy.salomepl.utils import SalomeManager
# click.echo("Exiting ...")
# if os.path.exists("anisotropy.pid"):
# os.remove("anisotropy.pid")
# SalomeManager().killall()
# sys.exit(0)

143
anisotropy/core/config.py Normal file
View File

@ -0,0 +1,143 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import os
import toml
import copy
import numpy as np
class Config(object):
def __init__(self):
self.content = {}
self.options = {}
self.params = {}
self.cases = None
def __getitem__(self, key):
return self.options[key]
def __setitem__(self, key, value):
self.options[key] = value
def __delitem__(self, key):
del self.options[key]
def update(self, **kwargs):
self.options.update(**kwargs)
def __len__(self):
return len(self.options)
def __iter__(self):
for key in self.options:
yield key
def load(self, filename: str):
path = os.path.abspath(filename)
if not os.path.exists(path):
raise FileNotFoundError(path)
self.content = toml.load(path)
self.options = copy.deepcopy(self.content["options"])
self.content.pop("options")
def dump(self, filename: str):
path = os.path.abspath(filename)
ext = os.path.splitext(path)[1]
if not ext == ".toml":
path += ".toml"
os.makedirs(os.path.split(path)[0], exist_ok = True)
self.content.update(options = self.options)
with open(path, "w") as io:
toml.dump(self.content, io, encoder = toml.TomlNumpyEncoder())
def minimize(self):
self.content = None
self.cases = None
def purge(self):
self.minimize()
self.params = None
def copy(self):
return copy.deepcopy(self)
def expand(self):
self.cases = []
# Expand structures for each direction and each alpha
for structure in self.content["structures"]:
structure = copy.deepcopy(structure)
alphaA = np.round(np.arange(
structure["alpha"][0],
structure["alpha"][1] + structure["alphaStep"],
structure["alphaStep"]
), 9)
directionA = np.array(structure["directions"], dtype = float)
structure.pop("alpha")
structure.pop("directions")
for direction in directionA:
for alpha in alphaA:
self.cases.append({
"alpha": alpha,
"direction": direction,
**structure
})
def chooseParams(self, idn: int):
if len(self.cases) > 0:
self.params = self.cases[idn]
else:
raise IndexError("list index out of range in cause of zero length of 'cases'")
def default_config() -> Config:
config = Config()
config.options = {
"nprocs": 1,
"stage": "all",
"overwrite": False,
"database": "anisotropy.db",
"build": "build",
"logs": "logs",
"shapefile": "shape.step",
"meshfile": "mesh.mesh"
}
config.content = {
"structures": []
}
labels = ["simple", "bodyCentered", "faceCentered"]
alphas = [[0.01, 0.28], [0.01, 0.17], [0.01, 0.13]]
for label, alpha in zip(labels, alphas):
config.content["structures"].append({
"label": label,
"alpha": alpha,
"alphaStep": 0.01,
"directions": [[1., 0., 0.], [0., 0., 1.], [1., 1., 1.]],
"r0": 1,
"filletsEnabled": True,
"pressureInlet": 1,
"pressureOutlet": 0,
"pressureInternal": 0,
"velocityInlet": [0., 0., 0.],
"velocityOutlet": None,
"velocityInternal": [0., 0., 0.],
"density": 1000,
"viscosity": 1e-3,
"scale": [1e-5, 1e-5, 1e-5]
})
return config

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()

View File

@ -1,580 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import os, sys
import time
from datetime import timedelta, datetime
import shutil
import logging
from copy import deepcopy
from math import sqrt
import toml
from anisotropy import (
__version__, env,
openfoam
)
from anisotropy.core.utils import setupLogger, Timer
from anisotropy.core.database import Database
from anisotropy import salomepl
import anisotropy.salomepl.utils
import anisotropy.salomepl.geometry
import anisotropy.salomepl.mesh
from anisotropy.samples import Simple, FaceCentered, BodyCentered
logger = logging.getLogger(env["logger_name"])
#setupLogger(logger, logging.INFO, env["LOG"])
#peeweeLogger = logging.getLogger("peewee")
#peeweeLogger.setLevel(logging.INFO)
class Anisotropy(object):
"""Ultimate class that organizes whole working process"""
def __init__(self):
"""Constructor method"""
self.env = env
self.db = None #Database(self.env["db_name"], self.env["db_path"])
self.params = []
def load(self, structure_type: str, structure_direction: list, structure_theta: float):
"""Shortcut for `Database.setup` and `Database.load`.
See :class:`anisotropy.core.database.Database` for more details.
"""
self.db.setup()
self.params = self.db.load(structure_type, structure_direction, structure_theta)
def update(self, params: dict = None):
"""Shortcut for `Database.setup` and `Database.update`.
See :class:`anisotropy.core.database.Database` for more details.
"""
self.db.setup()
self.db.update(self.params if not params else params)
@staticmethod
def version():
"""Returns versions of all used main programs
:return:
Versions joined by next line symbol
"""
versions = {
"anisotropy": __version__,
"Python": sys.version.split(" ")[0],
"Salome": "[missed]",
"OpenFOAM": "[missed]"
}
try:
versions["Salome"] = salomepl.utils.SalomeManager().version()
versions["OpenFOAM"] = openfoam.version()
except Exception:
pass
return "\n".join([ f"{ k }: { v }" for k, v in versions.items() ])
def loadFromScratch(self, configpath: str = None) -> list:
"""Loads parameters from configuration file and expands special values
:return:
List of dicts with parameters
"""
config = configpath or self.env["CONFIG"]
if not os.path.exists(config):
logger.error("Missed configuration file")
return
else:
logger.info(f"Configuration file: { config }")
buf = toml.load(config).get("structures")
paramsAll = []
for entry in buf:
# Shortcuts
_theta = entry["structure"]["theta"]
thetaMin = int(_theta[0] / _theta[2])
thetaMax = int(_theta[1] / _theta[2]) + 1
thetaList = list(
map(lambda n: n * _theta[2], range(thetaMin, thetaMax))
)
_thickness = entry["mesh"]["thickness"]
count = len(thetaList)
thicknessList = list(
map(lambda n: _thickness[0] + n * (_thickness[1] - _thickness[0]) / (count - 1), range(0, count))
)
for direction in entry["structure"]["directions"]:
for n, theta in enumerate(thetaList):
mesh = deepcopy(entry["mesh"])
mesh["thickness"] = thicknessList[n]
entryNew = {
"structure": dict(
type = entry["structure"]["type"],
theta = theta,
direction = [ float(num) for num in direction ],
filletsEnabled = entry["structure"]["filletsEnabled"]
),
"mesh": mesh,
"submesh": deepcopy(entry["submesh"]),
"meshresult": dict(),
"flow": deepcopy(entry["flow"]),
"flowapproximation": deepcopy(entry["flowapproximation"]),
"flowresult": dict(),
}
# For `type = fixedValue` only
_velocity = entryNew["flowapproximation"]["velocity"]["boundaryField"]["inlet"]["value"]
entryNew["flowapproximation"]["velocity"]["boundaryField"]["inlet"]["value"] = [
val * _velocity for val in entryNew["structure"]["direction"]
]
_velocity = entryNew["flow"]["velocity"]["boundaryField"]["inlet"]["value"]
entryNew["flow"]["velocity"]["boundaryField"]["inlet"]["value"] = [
val * _velocity for val in entryNew["structure"]["direction"]
]
paramsAll.append(entryNew)
return paramsAll
def evalParams(self):
"""Evals specific geometry(structure) parameters"""
structure = self.params.get("structure")
if not structure:
logger.error("Trying to eval empty parameters")
return
if structure["type"] == "simple":
thetaMin = 0.01
thetaMax = 0.28
r0 = 1
L = 2 * r0
radius = r0 / (1 - structure["theta"])
C1, C2 = 0.8, 0.5
Cf = C1 + (C2 - C1) / (thetaMax - thetaMin) * (structure["theta"] - thetaMin)
delta = 0.2
fillets = delta - Cf * (radius - r0)
elif structure["type"] == "faceCentered":
thetaMin = 0.01
thetaMax = 0.13
L = 1.0
r0 = L * sqrt(2) / 4
radius = r0 / (1 - structure["theta"])
C1, C2 = 0.3, 0.2
Cf = C1 + (C2 - C1) / (thetaMax - thetaMin) * (structure["theta"] - thetaMin)
delta = 0.012
fillets = delta - Cf * (radius - r0)
elif structure["type"] == "bodyCentered":
thetaMin = 0.01
thetaMax = 0.18
L = 1.0
r0 = L * sqrt(3) / 4
radius = r0 / (1 - structure["theta"])
C1, C2 = 0.3, 0.2
Cf = C1 + (C2 - C1) / (thetaMax - thetaMin) * (structure["theta"] - thetaMin)
delta = 0.02
fillets = delta - Cf * (radius - r0)
self.params["structure"].update(
L = L,
r0 = r0,
radius = radius,
fillets = fillets
)
def getCasePath(self, path: str = None) -> str:
"""Constructs case path from control parameters
:return: Absolute path to case
:rtype: str
"""
structure = self.params.get("structure")
if not structure:
logger.error("Trying to use empty parameters")
return
if path:
path = os.path.join(path, "build")
else:
path = self.env["BUILD"]
return os.path.join(
path,
structure["type"],
"direction-{}".format(str(structure['direction']).replace(" ", "")),
f"theta-{ structure['theta'] }"
)
def computeMesh(self, path):
"""Computes a mesh on shape via Salome
:return: Process output, error messages and returncode
:rtype: tuple(str, str, int)
"""
p = self.params["structure"]
scriptpath = os.path.join(self.env["ROOT"], "anisotropy/core/cli.py")
salomeargs = [
"computemesh",
p["type"],
p["direction"],
p["theta"],
path
]
manager = salomepl.utils.SalomeManager()
casepath = self.getCasePath(path)
self.params["meshresult"]["meshStatus"] = "Computing"
self.update()
timer = Timer()
out, err, returncode = manager.execute(
scriptpath,
*salomeargs,
timeout = self.env["salome_timeout"],
root = self.env["ROOT"],
logpath = casepath
)
self.load(p["type"], p["direction"], p["theta"])
if not returncode:
self.params["meshresult"].update(
meshStatus = "Done",
meshCalculationTime = timer.elapsed()
)
else:
self.params["meshresult"].update(
meshStatus = "Failed"
)
self.update()
return out, err, returncode
def genmesh(self, path):
"""Computes a mesh on shape
Warning: Working only inside Salome Environment
"""
setupLogger(logger, logging.INFO, self.env["LOG"])
p = self.params
###
# Shape
##
logger.info("Constructing shape ...")
geompy = salomepl.geometry.getGeom()
structure = dict(
simple = Simple,
bodyCentered = BodyCentered,
faceCentered = FaceCentered
)[p["structure"]["type"]]
shapeGeometry = structure(**p["structure"])
shape, groups = shapeGeometry.build()
[length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06)
###
# Mesh
##
logger.info("Prepairing mesh ...")
mp = p["mesh"]
lengths = [
geompy.BasicProperties(edge)[0] for edge in geompy.SubShapeAll(shape, geompy.ShapeType["EDGE"])
]
meanSize = sum(lengths) / len(lengths)
mp["maxSize"] = meanSize
mp["minSize"] = meanSize * 1e-1
mp["chordalError"] = mp["maxSize"] / 2
faces = []
for group in groups:
if group.GetName() in mp["facesToIgnore"]:
faces.append(group)
mesh = salomepl.mesh.Mesh(shape)
mesh.Tetrahedron(**mp)
if mp["viscousLayers"]:
mesh.ViscousLayers(**mp, faces = faces)
# Submesh
smp = p["submesh"]
for submesh in smp:
for group in groups:
if submesh["name"] == group.GetName():
subshape = group
submesh["maxSize"] = meanSize * 1e-1
submesh["minSize"] = meanSize * 1e-3
submesh["chordalError"] = submesh["minSize"] * 1e+1
mesh.Triangle(subshape, **submesh)
self.update()
logger.info("Computing mesh ...")
out, err, returncode = mesh.compute()
###
# Results
##
#p["meshresult"] = dict()
if not returncode:
mesh.removePyramids()
mesh.assignGroups()
casePath = self.getCasePath(path)
os.makedirs(casePath, exist_ok = True)
logger.info("Exporting mesh ...")
returncode, err = mesh.exportUNV(os.path.join(casePath, "mesh.unv"))
if returncode:
logger.error(err)
meshStats = mesh.stats()
p["meshresult"].update(
surfaceArea = surfaceArea,
volume = volume,
volumeCell = shapeGeometry.volumeCell,
**meshStats
)
self.update()
else:
logger.error(err)
p["meshresult"].update(
surfaceArea = surfaceArea,
volume = volume,
volumeCell = shapeGeometry.volumeCell
)
self.update()
def computeFlow(self, path):
"""Computes a flow on mesh via OpenFOAM
:return:
Process output, error messages and returncode
"""
###
# Case preparation
##
foamCase = [ "0", "constant", "system" ]
#self.params["flowresult"] = dict()
self.params["flowresult"]["flowStatus"] = "Computing"
self.update()
timer = Timer()
flow = self.params["flow"]
flowapproximation = self.params["flowapproximation"]
# ISSUE: ideasUnvToFoam cannot import mesh with '-case' flag so 'os.chdir' for that
casePath = self.getCasePath(path)
if not os.path.exists(casePath):
err = f"Cannot find case path { casePath }"
self.params["flowresult"]["flowStatus"] = "Failed"
self.update()
return "", err, 1
os.chdir(casePath)
openfoam.foamClean()
for d in foamCase:
shutil.copytree(
os.path.join(self.env["openfoam_template"], d),
os.path.join(casePath, d)
)
###
# Mesh manipulations
##
if not os.path.exists("mesh.unv"):
os.chdir(path or self.env["ROOT"])
err = f"Missed 'mesh.unv'"
self.params["flowresult"]["flowStatus"] = "Failed"
self.update()
return "", err, 1
out, err, returncode = openfoam.ideasUnvToFoam("mesh.unv")
if returncode:
os.chdir(path or self.env["ROOT"])
self.params["flowresult"]["flowStatus"] = "Failed"
self.update()
return out, err, returncode
openfoam.createPatch(dictfile = "system/createPatchDict")
openfoam.foamDictionary(
"constant/polyMesh/boundary",
"entry0.defaultFaces.type",
"wall"
)
openfoam.foamDictionary(
"constant/polyMesh/boundary",
"entry0.defaultFaces.inGroups",
"1 (wall)"
)
out, err, returncode = openfoam.checkMesh()
if out: logger.warning(out)
openfoam.transformPoints(flow["scale"])
###
# Decomposition and initial approximation
#
# NOTE: Temporary without decomposition
##
openfoam.foamDictionary(
"constant/transportProperties",
"nu",
str(flow["transportProperties"]["nu"])
)
# openfoam.decomposePar()
openfoam.renumberMesh()
pressureBF = flowapproximation["pressure"]["boundaryField"]
velocityBF = flowapproximation["velocity"]["boundaryField"]
openfoam.foamDictionary(
"0/p",
"boundaryField.inlet.value",
openfoam.uniform(pressureBF["inlet"]["value"])
)
openfoam.foamDictionary(
"0/p",
"boundaryField.outlet.value",
openfoam.uniform(pressureBF["outlet"]["value"])
)
openfoam.foamDictionary(
"0/U",
"boundaryField.inlet.value",
openfoam.uniform(velocityBF["inlet"]["value"])
)
openfoam.potentialFoam()
###
# Main computation
##
pressureBF = flow["pressure"]["boundaryField"]
velocityBF = flow["velocity"]["boundaryField"]
openfoam.foamDictionary(
"0/U",
"boundaryField.inlet.type",
velocityBF["inlet"]["type"]
)
openfoam.foamDictionary(
"0/U",
"boundaryField.inlet.value",
openfoam.uniform(velocityBF["inlet"]["value"])
)
#for n in range(os.cpu_count()):
# openfoam.foamDictionary(
# f"processor{n}/0/U",
# "boundaryField.inlet.type",
# velocityBF.inlet.type
# )
# openfoam.foamDictionary(
# f"processor{n}/0/U",
# "boundaryField.inlet.value",
# openfoam.uniform(velocityBF.inlet.value[direction])
# )
out, err, returncode = openfoam.simpleFoam()
if not returncode:
self.params["flowresult"]["flowCalculationTime"] = timer.elapsed()
self.params["flowresult"]["flowStatus"] = "Done"
else:
self.params["flowresult"]["flowStatus"] = "Failed"
self.update()
os.chdir(path or self.env["ROOT"])
return out, str(err, "utf-8"), returncode
def flowRate(self):
casePath = self.getCasePath()
foamPostProcessing = "postProcessing/flowRatePatch(name=outlet)/0/surfaceFieldValue.dat"
path = os.path.join(casePath, foamPostProcessing)
if not os.path.exists(path):
logger.warning(f"Unable to compute flow rate. Missed { path }")
return
with open(path, "r") as io:
lastLine = io.readlines()[-1]
flowRate = float(lastLine.replace(" ", "").replace("\n", "").split("\t")[1])
self.params["flowresult"]["flowRate"] = flowRate
self.update()
return flowRate
def porosity(self):
mr = self.params["meshresult"]
fr = self.params["flowresult"]
fr["porosity"] = mr["volume"] / mr["volumeCell"]
self.update()
return fr["porosity"]

View File

@ -1,174 +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
)
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 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)

View File

@ -0,0 +1,64 @@
# -*- coding: utf-8 -*-
import multiprocessing as mlp
import dill
class ParallelRunner(object):
def __init__(self, nprocs: int = 1, daemon: bool = True):
self.nprocs = nprocs
self.daemon = daemon
self.processes = []
self.queueInput = mlp.Queue(maxsize = 1)
self.queueOutput = mlp.Queue()
self.__pos = -1
self.output = []
def append(self, command, args = [], kwargs = {}):
self.__pos += 1
self.queueInput.put(dill.dumps((self.__pos, command, args, kwargs)))
def extend(self, commands: list, args: list = [], kwargs: list = []):
for command, cargs, ckwargs in zip(commands, args, kwargs):
self.append(command, cargs, ckwargs)
@staticmethod
def queueRelease(queueInput, queueOutput):
while True:
pos, command, args, kwargs = dill.loads(queueInput.get())
if pos is None or command is None:
break
output = command(*args, **kwargs)
queueOutput.put((pos, output))
def start(self):
for n in range(self.nprocs):
self.processes.append(mlp.Process(
target = self.queueRelease,
args = (self.queueInput, self.queueOutput),
name = f"worker-{ n + 1 }"
))
for proc in self.processes:
proc.daemon = self.daemon
proc.start()
def wait(self):
for _ in range(self.nprocs):
self.append(None)
self.output = [ [] for _ in range(self.queueOutput.qsize()) ]
for _ in range(self.queueOutput.qsize()):
pos, output = self.queueOutput.get()
self.output[pos] = output
for proc in self.processes:
proc.join()
self.__pos = -1

View File

@ -0,0 +1,97 @@
# -*- coding: utf-8 -*-
import pathlib
import numpy as np
import anisotropy.openfoam as of
def flowRate(patch: str, path: str = None):
func = f"patchFlowRate(patch={ patch })"
path = pathlib.Path(path or "").resolve()
outfile = path / "postProcessing" / func / "0/surfaceFieldValue.dat"
of.commands.postProcess(func, cwd = path, logpath = path / "patchFlowRate.log")
surfaceFieldValue = of.conversion.read_dat(outfile)["sum(phi)"][-1]
return surfaceFieldValue
def permeability(viscosity = None, flowRate = None, length = None, areaCross = None, pressureInlet = None, pressureOutlet = None, **kwargs):
viscosity = kwargs.get("viscosity", viscosity)
flowRate = kwargs.get("flowRate", flowRate)
length = kwargs.get("length", length)
areaCross = kwargs.get("areaCross", areaCross)
pressureInlet = kwargs.get("pressureInlet", pressureInlet)
pressureOutlet = kwargs.get("pressureOutlet", pressureOutlet)
return viscosity * length * flowRate / (areaCross * (pressureInlet - pressureOutlet))
def mean_nan(arr):
temp = arr.copy()
if np.isnan(temp[0]):
temp[0] = temp[1]
for n, item in enumerate(temp):
if np.all(np.isnan(item)):
vals = temp[n - 1 : n + 2]
if np.sum(~np.isnan(vals)) <= 1:
vals = temp[n - 2 : n + 3]
temp[n] = vals[~np.isnan(vals)].mean()
return temp
def filter_group(arr, nan = True, qhigh = True, quantile = 0.97):
temp = arr.copy()
check = True
quan = np.quantile(temp[~np.isnan(temp)], quantile)
limit = 1000
while check:
if nan and np.any(np.isnan(temp)):
temp = mean_nan(temp)
check = True
elif qhigh and np.any(quan < temp):
temp[quan < temp] = np.nan
check = True
else:
check = False
if limit <= 0:
break
else:
limit -= 1
return temp
def pad_nan(x, y):
return np.pad(y, (0, x.size - y.size), 'constant', constant_values = np.nan)
"""
def multiplot_y(
x,
y,
labels = [],
lines = [],
axes_labels = ["", ""],
legend = True,
grid = True,
figsize = (12, 6)
):
fig, ax = plt.subplots(nrows = 1, ncols = 1, figsize = figsize)
for n, y_arr in enumerate(y):
ax.plot(x, y_arr, )
"""

388
anisotropy/core/runner.py Normal file
View File

@ -0,0 +1,388 @@
# -*- coding: utf-8 -*-
from os import PathLike
from os import environ
from datetime import datetime
import pathlib
import logging
import numpy as np
from anisotropy.database import Database, tables
from anisotropy import shaping
from anisotropy import meshing
from anisotropy import solving
from anisotropy import openfoam
from . import config as core_config
from . import postprocess
from . import utils
from . import ParallelRunner
logger = logging.getLogger(__name__)
class UltimateRunner(object):
def __init__(self, config = None, exec_id: int = None, _type: str = "master"):
self.type = _type
# Configuration file
self.config = config or core_config.default_config()
# Database preparation
self.database = Database(self.config["database"])
self.exec_id = None
if exec_id is not None:
if self.database.getExecution(int(exec_id)):
self.exec_id = int(exec_id)
else:
logger.warning(f"Execution id '{ exec_id }' not found. Creating new.")
if self.exec_id is None:
with self.database:
self.exec_id = tables.Execution.create(date = datetime.now()).exec_id
if self.type == "master":
logger.info(f"Current execution id: { self.exec_id }")
# Parameters
self.queue = []
def _query_database(self, tableName: str, to_dict: bool = False):
tableName = tableName.lower()
get_row = {
"shape": self.database.getShape,
"mesh": self.database.getMesh,
"flowonephase": self.database.getFlowOnephase
}[tableName]
# query
args = (
self.config.params["label"],
self.config.params["direction"],
self.config.params["alpha"],
self.exec_id,
to_dict
)
return get_row(*args)
def prepare_database(self):
# create a row in each table for the current case
# shape
with self.database:
shape = self.database.getShape(execution = self.exec_id, **self.config.params)
if not shape:
shape = tables.Shape(exec_id = self.exec_id, **self.config.params)
self.database.csave(shape)
# mesh
with self.database:
mesh = self.database.getMesh(execution = self.exec_id, **self.config.params)
if not mesh:
mesh = tables.Mesh(shape_id = shape.shape_id)
self.database.csave(mesh)
# onephase flow
with self.database:
flow = self.database.getFlowOnephase(execution = self.exec_id, **self.config.params)
if not flow:
flow = tables.FlowOnephase(mesh_id = mesh.mesh_id, **self.config.params)
self.database.csave(flow)
def prepare_queue(self):
self.config.expand()
logger.info(f"Preparing queue: { len(self.config.cases) }")
for idn, case in enumerate(self.config.cases):
config = self.config.copy()
config.chooseParams(idn)
config.minimize()
kwargs = {
"config": config,
"exec_id": self.exec_id
}
self.queue.append(kwargs)
def start(self, queue: list = None, nprocs: int = None):
nprocs = nprocs or self.config["nprocs"]
logger.info(f"Starting subprocesses: { nprocs }")
parallel = ParallelRunner(nprocs = nprocs)
parallel.start()
for kwargs in self.queue:
parallel.append(self.subrunner, kwargs = kwargs)
parallel.wait()
@property
def casepath(self) -> PathLike:
params = self.config.params
path = pathlib.Path(environ["AP_CWD"], environ["AP_BUILD_DIR"])
path /= "execution-{}".format(self.exec_id)
path /= "{}-[{},{},{}]-{}".format(
params["label"],
*[ str(d) for d in params["direction"] ],
params["alpha"]
)
return path.resolve()
@property
def shapefile(self) -> PathLike:
return self.casepath / self.config["shapefile"]
@property
def meshfile(self) -> PathLike:
return self.casepath / self.config["meshfile"]
def compute_shape(self):
params = self.config.params
shapeParams = self._query_database("shape")
logger.info("Computing shape for {} with direction = {} and alpha = {}".format(
params["label"], params["direction"], params["alpha"]
))
timer = utils.Timer()
# check physical existence
if self.shapefile.exists() and self.config["overwrite"] is not True:
logger.info("Shape exists, skipping ...")
return
#
generate_shape = {
"simple": shaping.primitives.simple,
"bodyCentered": shaping.primitives.body_centered,
"faceCentered": shaping.primitives.face_centered
}[shapeParams.label]
try:
# generate shape
shape = generate_shape(
shapeParams.alpha,
shapeParams.direction,
r0 = shapeParams.r0,
filletsEnabled = shapeParams.filletsEnabled
)
# export
self.casepath.mkdir(parents = True, exist_ok = True)
shape.write(self.shapefile)
except Exception as e:
logger.error(e, exc_info = True)
shapeParams.shapeStatus = "failed"
else:
shapeParams.shapeStatus = "done"
shapeParams.volume = shape.shape.volume * np.prod(params["scale"])
shapeParams.volumeCell = shape.cell.volume * np.prod(params["scale"])
shapeParams.radius = shape.radius
shapeParams.fillets = shape.fillets
shapeParams.porosity = shapeParams.volume / shapeParams.volumeCell
shapeParams.areaOutlet = np.sum([
face.mass for face in shape.shape.faces if face.name == "outlet"
]) * params["scale"][0] ** 2
shapeParams.length = shape.length * params["scale"][0]
shapeParams.areaCellOutlet = np.sum([
face.mass for face in shape.cell.faces if face.name == "outlet"
]) * params["scale"][0] ** 2
# commit parameters
shapeParams.shapeExecutionTime = timer.elapsed()
self.database.csave(shapeParams)
def compute_mesh(self):
params = self.config.params
meshParams = self._query_database("mesh")
logger.info("Computing mesh for {} with direction = {} and alpha = {}".format(
params["label"], params["direction"], params["alpha"]
))
timer = utils.Timer()
# check physical existence
if self.meshfile.exists() and self.config["overwrite"] is not True:
logger.info("Mesh exists, skipping ...")
return
elif not self.shapefile.exists():
logger.error("Cannot find shape file to build a mesh.")
return
# Shape
shape = None
try:
# load shape
shape = shaping.Shape().read(self.shapefile)
# generate mesh
parameters = meshing.MeshingParameters()
mesh = meshing.Mesh(shape)
mesh.generate(parameters)
# export
self.casepath.mkdir(parents = True, exist_ok = True)
mesh.write(self.meshfile)
mesh.write(self.casepath / "mesh.msh")
except Exception as e:
logger.error(e, exc_info = True)
meshParams.meshStatus = "failed"
else:
meshParams.meshStatus = "done"
meshParams.faces = len(mesh.faces)
meshParams.volumes = len(mesh.volumes)
# commit parameters
meshParams.meshExecutionTime = timer.elapsed()
self.database.csave(meshParams)
def compute_flow(self):
params = self.config.params
flowParams = self._query_database("flowonephase")
logger.info("Computing flow for {} with direction = {} and alpha = {}".format(
params["label"], params["direction"], params["alpha"]
))
timer = utils.Timer()
# check physical existence
if openfoam.FoamCase(path = self.casepath).is_converged() and not self.config["overwrite"]:
logger.info("Solution exists. Skipping ...")
return
elif not self.shapefile.exists():
logger.error("Cannot find shape file to compute a solution.")
return
elif not self.meshfile.exists():
logger.error("Cannot find mesh file to compute a solution.")
return
# precalculate some parameters
flowParams.viscosityKinematic = flowParams.viscosity / flowParams.density
self.database.csave(flowParams)
#
flowParamsDict = self._query_database("flowonephase", to_dict = True)
try:
# load shape
shape = shaping.Shape().read(self.shapefile)
#
flow = solving.FlowOnePhase(
path = self.casepath,
direction = params["direction"],
patches = shape.patches(group = True, shiftIndex = True, prefix = "patch"),
scale = params["scale"],
**flowParamsDict
)
# generate solution
flow.generate()
except Exception as e:
logger.error(e, exc_info = True)
flowParams.flowStatus = "failed"
else:
flowParams.flowStatus = "done"
# commit parameters
flowParams.flowExecutionTime = timer.elapsed()
self.database.csave(flowParams)
def compute_postprocess(self):
params = self.config.params
flowParams = self._query_database("flowonephase")
logger.info("Computing post process for {} with direction = {} and alpha = {}".format(
params["label"], params["direction"], params["alpha"]
))
if flowParams.flowStatus == "done":
if flowParams.flowRate is None and self.config["overwrite"] is not True:
flowParams.flowRate = postprocess.flowRate("outlet", self.casepath)
self.database.csave(flowParams)
def pipeline(self, stage: str = None):
stage = stage or self.config["stage"]
stages = ["shape", "mesh", "flow", "postProcess", "all"]
if stage not in stages:
logger.error(f"Unknown stage '{ stage }'.")
return
for current in stages:
if current == "shape":
params = self._query_database("shape")
# check database entry
if params.shapeStatus == "done" and self.config["overwrite"] is not True:
logger.info("Successful shape entry exists, skipping ...")
continue
self.compute_shape()
if current == "mesh":
params = self._query_database("mesh")
# check database entry
if params.meshStatus == "done" and self.config["overwrite"] is not True:
logger.info("Successful mesh entry exists, skipping ...")
continue
self.compute_mesh()
if current == "flow":
params = self._query_database("flowonephase")
# check database entry
if params.flowStatus == "done" and self.config["overwrite"] is not True:
logger.info("Successful flow entry exists, skipping ...")
continue
self.compute_flow()
if current == "postProcess":
params = self._query_database("flowonephase")
# check database entry
# if params.flowStatus == "done" and self.config["overwrite"] is not True:
# logger.info("Successful flow entry exists, skipping ...")
# continue
self.compute_postprocess()
if current == stage or current == "all":
break
@staticmethod
def subrunner(*args, **kwargs):
runner = UltimateRunner(
config = kwargs["config"],
exec_id = kwargs["exec_id"],
_type = "worker"
)
runner.prepare_database()
runner.pipeline()

View File

@ -1,286 +1,113 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import logging
from multiprocessing import Queue, Process, cpu_count
import socket
import copy
import time
from types import FunctionType
import os
import contextlib
class CustomFormatter(logging.Formatter):
def _getFormat(self, level: int):
def __init__(self, colors: bool = True):
self.colors = colors
def applyColors(self, level: int, info: str, msg: str):
grey = "\x1b[38;21m"
yellow = "\x1b[33;21m"
red = "\x1b[31;21m"
bold_red = "\x1b[31;1m"
reset = "\x1b[0m"
format = "[ %(asctime)s ] [ %(processName)s ] [ %(levelname)s ] %(message)s"
formats = {
logging.DEBUG: grey + format + reset,
logging.INFO: grey + format + reset,
logging.WARNING: yellow + format + reset,
logging.ERROR: red + format + reset,
logging.CRITICAL: bold_red + format + reset
logging.DEBUG: grey + info + reset + msg,
logging.INFO: grey + info + reset + msg,
logging.WARNING: yellow + info + reset + msg,
logging.ERROR: red + info + reset + msg,
logging.CRITICAL: bold_red + info + reset + msg
}
return formats.get(level)
def format(self, record):
log_fmt = self._getFormat(record.levelno)
time_fmt = "%H:%M:%S %d-%m-%y"
info = "[%(levelname)s %(processName)s %(asctime)s %(funcName)s]" # [ %(processName)s ]
msg = " %(message)s"
if self.colors:
log_fmt = self.applyColors(record.levelno, info, msg)
else:
log_fmt = info + msg
time_fmt = "%d-%m-%y %H:%M:%S"
formatter = logging.Formatter(log_fmt, time_fmt)
return formatter.format(record)
def setupLogger(logger, level: int, filepath: str = None):
"""Applies settings to logger
def setupLogger(level: int, filepath: str = None):
"""Applies settings to root logger
:param logger:
Instance of :class:`logging.Logger`
:param level:
Logging level (logging.INFO, logging.WARNING, ..)
:param filepath:
Path to directory
"""
logger.handlers = []
logger.setLevel(level)
logging.addLevelName(logging.INFO, "II")
logging.addLevelName(logging.WARNING, "WW")
logging.addLevelName(logging.ERROR, "EE")
logging.addLevelName(logging.CRITICAL, "CC")
streamhandler = logging.StreamHandler()
streamhandler.setLevel(level)
streamhandler.setFormatter(CustomFormatter())
logger.addHandler(streamhandler)
if filepath:
if not os.path.exists(filepath):
os.makedirs(filepath, exist_ok = True)
filehandler = logging.FileHandler(
os.path.join(filepath, "{}.log".format(logger.name))
)
filehandler.setLevel(level)
filehandler.setFormatter(CustomFormatter())
logger.addHandler(filehandler)
class struct:
def __init__(self, *args, **kwargs):
if len(args) > 0:
if type(args[0]) == dict:
for (k, v) in args[0].items():
if type(v) == dict:
setattr(self, k, struct(v))
else:
setattr(self, k, v)
else:
self.__dict__.update(kwargs)
def __iter__(self):
for k in self.__dict__:
if type(getattr(self, k)) == struct:
yield k, dict(getattr(self, k))
else:
yield k, getattr(self, k)
def __str__(self):
members = []
for key in self.__dict__.keys():
members.append(f"{ key } = ")
if type(self.__dict__[key]) == str:
members[len(members) - 1] += f"\"{ self.__dict__[key] }\""
else:
members[len(members) - 1] += f"{ self.__dict__[key] }"
return f"struct({', '.join(members)})"
def __repr__(self):
return str(self)
def deepupdate(target, src):
for k, v in src.items():
if isinstance(v, dict):
if not k in target:
target[k] = copy.deepcopy(v)
else:
deepupdate(target[k], v)
else:
target[k] = copy.copy(v)
def collapse(source, key = None, level = 0, sep = "_"):
if isinstance(source, dict) and source:
level = level + 1
res = {}
for k, v in source.items():
ret, lvl = collapse(v, k, level)
for kk,vv in ret.items():
if level == lvl:
newkey = k
else:
newkey = "{}{}{}".format(k, sep, kk)
res.update({ newkey: vv })
if level == 1:
return res
else:
return res, level
else:
return { key: source }, level
def expand(source, sep = "_"):
res = {}
for k, v in source.items():
if k.find(sep) == -1:
res.update({ k: v })
else:
keys = k.split(sep)
cur = res
for n, kk in enumerate(keys):
if not len(keys) == n + 1:
if not cur.get(kk):
cur.update({ kk: {} })
cur = cur[kk]
else:
cur[kk] = v
return res
#if os.path.exists(env["CONFIG"]):
# config = toml.load(env["CONFIG"])
# for restricted in ["ROOT", "BUILD", "LOG", "CONFIG"]:
# if config.get(restricted):
# config.pop(restricted)
# TODO: not working if custom config empty and etc
# for m, structure in enumerate(config["structures"]):
# for n, estructure in enumerate(env["structures"]):
# if estructure["name"] == structure["name"]:
# deepupdate(env["structures"][n], config["structures"][m])
# config.pop("structures")
# deepupdate(env, config)
def timer(func: FunctionType) -> (tuple, float):
"""(Decorator) Returns output of inner function and execution time
logging.root.addHandler(streamhandler)
:param func: inner function
:type func: FunctionType
logging.root.setLevel(level)
if filepath:
filehandler = logging.FileHandler(filepath)
filehandler.setLevel(logging.INFO)
filehandler.setFormatter(CustomFormatter(colors = False))
:return: output, elapsed time
:rtype: tuple(tuple, float)
"""
def inner(*args, **kwargs):
start = time.monotonic()
ret = func(*args, **kwargs)
elapsed = time.monotonic() - start
return ret, elapsed
return inner
logging.root.addHandler(filehandler)
class Timer(object):
def __init__(self):
self.start = time.monotonic()
self.update()
def update(self):
self.start = time.monotonic()
def elapsed(self):
return time.monotonic() - self.start
def queue(cmd, qin, qout, *args):
while True:
# Get item from the queue
pos, var = qin.get()
class ErrorHandler(contextlib.AbstractContextManager):
def __init__(self):
self.error = ""
self.returncode = 0
self.traceback = None
def __enter__(self):
return self, self.handle
def __exit__(self, exc_type, exc_value, traceback):
if exc_type:
self.error = exc_value.args
self.returncode = 1
self.traceback = traceback
def handle(self, obj):
def inner(*args, **kwargs):
try:
output = obj(*args, **kwargs)
except Exception as e:
self.error = e.args
self.returncode = 1
else:
return output
# Exit point
if pos is None:
break
# Execute command
res = cmd(*var, *args)
# Put results to the queue
qout.put((pos, res))
return
def parallel(np, var, cmd):
varcount = len(var)
processes = []
nprocs = np if np <= cpu_count() else cpu_count()
qin = Queue(1)
qout = Queue()
# Create processes
for n in range(nprocs):
pargs = [cmd, qin, qout]
p = Process(target = queue, args = tuple(pargs))
processes.append(p)
# Start processes
for p in processes:
p.daemon = True
p.start()
# Fill queue
for n in range(varcount):
qin.put((n, var[n]))
for _ in range(nprocs):
qin.put((None, None))
# Get results
results = [[] for n in range(varcount)]
for n in range(varcount):
index, res = qout.get()
results[index] = res
# Wait until each processor has finished
for p in processes:
p.join()
return results
def portIsFree(address, port):
with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
return s.connect_ex((address, port)) == 0
return inner

View File

@ -0,0 +1,13 @@
# -*- coding: utf-8 -*-
from . import utils
from . import tables
from .database import Database
__all__ = [
"utils",
"tables",
"Database"
]

View File

@ -0,0 +1,309 @@
# -*- coding: utf-8 -*-
from __future__ import annotations
from numpy import ndarray
import peewee as pw
import pathlib
import time
import pandas as pd
from . import tables
class Database(pw.SqliteDatabase):
def __init__(self, filename: str = None, *args, **kwargs):
"""A Database object contains SQLite database with convient
properties and methods
"""
self.path = filename
self.pragmas_ = kwargs.get("pragmas", { "foreign_keys": 1, "journal_mode": "wal" })
self.field_types_ = kwargs.get("field_types", { "list": "text" })
self.autoconnect_ = kwargs.get("autoconnect", False)
pw.SqliteDatabase.__init__(
self,
None,
pragmas = self.pragmas_,
field_types = self.field_types_,
autoconnect = self.autoconnect_
)
if self.path:
self.setup()
@property
def tables(self) -> list:
"""Return all tables as list.
"""
return [ tables.__dict__[table] for table in tables.__all__ ]
def setup(self, filename: str = None):
"""Initialize database and create tables.
:param filename:
Path to the file.
:return:
Self.
"""
self.path = pathlib.Path(filename or self.path).resolve()
self.init(
self.path,
pragmas = self.pragmas_
)
tables.database_proxy.initialize(self)
with self:
self.create_tables(self.tables)
return self
def csave(self, table: pw.Model, tries: int = 5000):
"""Try to save data from model to the database ignoring
peewee.OperationalError. Usefull for concurrent processes.
:param table:
Table to save.
:param tries:
Number of tries. Falling to sleep for 1 second if database
is locked.
"""
while tries >= 0:
if self.is_closed():
self.connect()
try:
table.save()
except pw.OperationalError:
tries -= 1
time.sleep(1)
else:
self.close()
break
def getExecution(
self,
idn: int,
to_dict: bool = False,
**kwargs
) -> tables.Execution | None:
"""Get execution entry from database.
:param idn:
Index of the execution.
:return:
If entry is found returns Model instance else None.
"""
query = tables.Execution.select().where(tables.Execution.exec_id == idn)
with self:
if to_dict:
table = query.dicts().get() if query.exists() else None
else:
table = query.get() if query.exists() else None
return table
def getLatest(self) -> tables.Execution | None:
"""Get latest execution entry from database.
:return:
If entry is found returns Model instance else None.
"""
query = tables.Execution.select()
with self:
table = query[-1] if query.exists() else None
return table
def getShape(
self,
label: str = None,
direction: list[float] | ndarray = None,
alpha: float = None,
execution: int = None,
to_dict: bool = False,
**kwargs
) -> tables.Shape | None:
"""Get shape entry from database.
:param label:
Label of the shape.
:param direction:
Array of floats represents direction vector.
:param alpha:
Spheres overlap parameter.
:param execution:
Index of the execution. If None, use latest.
:return:
If entry is found returns Model instance else None.
"""
execution = execution or self.getLatest()
query = (
tables.Shape
.select()
.join(tables.Execution, pw.JOIN.LEFT_OUTER)
.where(
tables.Execution.exec_id == execution,
tables.Shape.label == label,
tables.Shape.direction == direction,
tables.Shape.alpha == alpha
)
)
with self:
if to_dict:
table = query.dicts().get() if query.exists() else None
else:
table = query.get() if query.exists() else None
return table
def getMesh(
self,
label: str = None,
direction: list[float] | ndarray = None,
alpha: float = None,
execution: int = None,
to_dict: bool = False,
**kwargs
) -> tables.Mesh | None:
"""Get mesh entry from database.
:param label:
Label of the shape.
:param direction:
Array of floats represents direction vector.
:param alpha:
Spheres overlap parameter.
:param execution:
Index of the execution. If None, use latest.
:return:
If entry is found returns Model instance else None.
"""
execution = execution or self.getLatest()
query = (
tables.Mesh
.select()
.join(tables.Shape, pw.JOIN.LEFT_OUTER)
.join(tables.Execution, pw.JOIN.LEFT_OUTER)
.where(
tables.Execution.exec_id == execution,
tables.Shape.label == label,
tables.Shape.direction == direction,
tables.Shape.alpha == alpha
)
)
with self:
if to_dict:
table = query.dicts().get() if query.exists() else None
else:
table = query.get() if query.exists() else None
return table
def getFlowOnephase(
self,
label: str = None,
direction: list[float] | ndarray = None,
alpha: float = None,
execution: int = None,
to_dict: bool = False,
**kwargs
) -> tables.Mesh | dict | None:
"""Get one phase flow entry from database.
:param label:
Label of the shape.
:param direction:
Array of floats represents direction vector.
:param alpha:
Spheres overlap parameter.
:param execution:
Index of the execution. If None, use latest.
:param to_dict:
If True, convert result to dict.
:return:
If entry is found returns Model instance or dict else None.
"""
execution = execution or self.getLatest()
query = (
tables.FlowOnephase
.select()
.join(tables.Mesh, pw.JOIN.LEFT_OUTER)
.join(tables.Shape, pw.JOIN.LEFT_OUTER)
.join(tables.Execution, pw.JOIN.LEFT_OUTER)
.where(
tables.Execution.exec_id == execution,
tables.Shape.label == label,
tables.Shape.direction == direction,
tables.Shape.alpha == alpha
)
)
with self:
if to_dict:
table = query.dicts().get() if query.exists() else None
else:
table = query.get() if query.exists() else None
return table
def load(
self,
field: str,
execution: int = None,
label: str = None,
direction: list[float] | ndarray = None
):
execution = execution or self.getLatest()
for table in self.tables:
try:
column = getattr(table, field)
except AttributeError:
pass
else:
break
query = table.select(tables.Shape.alpha, column, tables.Shape.direction, tables.Shape.label)
for current in reversed(self.tables[ :self.tables.index(table)]):
query = query.join(table, pw.JOIN.LEFT_OUTER)
query = (
query.switch(tables.Shape)
.where(tables.Shape.exec_id == execution)
.order_by(tables.Shape.label, tables.Shape.direction, tables.Shape.alpha)
)
resp = None
with self:
if query.exists():
resp = []
for row in query.dicts():
for k in row.keys():
if type(row[k]) == list:
row[k] = str(row[k])
resp.append(row)
resp = pd.DataFrame(resp)
if label is not None:
resp = resp[resp.label == label]
if direction is not None:
resp = resp[resp.direction == direction]
return resp[field]

View File

@ -0,0 +1,106 @@
# -*- coding: utf-8 -*-
import peewee as pw
from . import utils
# proxy, assign database later
database_proxy = pw.Proxy()
class Execution(pw.Model):
exec_id = pw.AutoField()
date = pw.DateTimeField()
executionTime = pw.TimeField(null = True)
class Meta:
database = database_proxy
table_name = "executions"
class Shape(pw.Model):
shape_id = pw.AutoField()
exec_id = pw.ForeignKeyField(Execution, backref = "executions", on_delete = "CASCADE")
shapeStatus = pw.TextField(null = True, default = "idle")
shapeExecutionTime = pw.TimeField(null = True)
label = pw.TextField(null = True)
direction = utils.JSONField(null = True)
alpha = pw.FloatField(null = True)
r0 = pw.FloatField(null = True)
L = pw.FloatField(null = True)
radius = pw.FloatField(null = True)
filletsEnabled = pw.BooleanField(null = True)
fillets = pw.FloatField(null = True)
volumeCell = pw.FloatField(null = True)
volume = pw.FloatField(null = True)
porosity = pw.FloatField(null = True)
areaOutlet = pw.FloatField(null = True)
length = pw.FloatField(null = True)
areaCellOutlet = pw.FloatField(null = True)
class Meta:
database = database_proxy
table_name = "shapes"
# depends_on = Execution
class Mesh(pw.Model):
mesh_id = pw.AutoField()
shape_id = pw.ForeignKeyField(Shape, backref = "shapes", on_delete = "CASCADE")
meshStatus = pw.TextField(null = True, default = "idle")
meshExecutionTime = pw.TimeField(null = True)
elements = pw.IntegerField(null = True)
edges = pw.IntegerField(null = True)
faces = pw.IntegerField(null = True)
volumes = pw.IntegerField(null = True)
tetrahedrons = pw.IntegerField(null = True)
prisms = pw.IntegerField(null = True)
pyramids = pw.IntegerField(null = True)
class Meta:
database = database_proxy
table_name = "meshes"
# depends_on = Execution
class FlowOnephase(pw.Model):
flow_id = pw.AutoField()
mesh_id = pw.ForeignKeyField(Mesh, backref = "meshes", on_delete = "CASCADE")
flowStatus = pw.TextField(null = True, default = "idle")
flowExecutionTime = pw.TimeField(null = True)
pressureInlet = pw.FloatField(null = True)
pressureOutlet = pw.FloatField(null = True)
pressureInternal = pw.FloatField(null = True)
velocityInlet = utils.JSONField(null = True)
velocityOutlet = utils.JSONField(null = True)
velocityInternal = utils.JSONField(null = True)
viscosity = pw.FloatField(null = True)
viscosityKinematic = pw.FloatField(null = True)
density = pw.FloatField(null = True)
flowRate = pw.FloatField(null = True)
permeability = pw.FloatField(null = True)
class Meta:
database = database_proxy
table_name = "flows"
# depends_on = Execution
__all__ = [
"Execution",
"Shape",
"Mesh",
"FlowOnephase"
]

View File

@ -0,0 +1,160 @@
# -*- coding: utf-8 -*-
from numpy import ndarray
import peewee as pw
import json
class ListField(pw.TextField):
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 Exception:
pval.append(entry.strip().replace("'", ""))
return pval
class JSONPath(pw.ColumnBase):
def __init__(self, field, path = None):
super(JSONPath, self).__init__()
self._field = field
self._path = path or ()
@property
def path(self):
return pw.Value('$%s' % ''.join(self._path))
def __getitem__(self, idx):
if isinstance(idx, int):
item = '[%s]' % idx
else:
item = '.%s' % idx
return JSONPath(self._field, self._path + (item,))
def set(self, value, as_json = None):
if as_json or isinstance(value, (list, dict)):
value = pw.fn.json(self._field._json_dumps(value))
return pw.fn.json_set(self._field, self.path, value)
def update(self, value):
return self.set(pw.fn.json_patch(self, self._field._json_dumps(value)))
def remove(self):
return pw.fn.json_remove(self._field, self.path)
def json_type(self):
return pw.fn.json_type(self._field, self.path)
def length(self):
return pw.fn.json_array_length(self._field, self.path)
def children(self):
return pw.fn.json_each(self._field, self.path)
def tree(self):
return pw.fn.json_tree(self._field, self.path)
def __sql__(self, ctx):
return ctx.sql(
pw.fn.json_extract(self._field, self.path)
if self._path else self._field
)
class JSONField(pw.TextField):
field_type = 'TEXT'
unpack = False
def __init__(self, json_dumps = None, json_loads = None, **kwargs):
super(JSONField, self).__init__(**kwargs)
self._json_dumps = json_dumps or json.dumps
self._json_loads = json_loads or json.loads
def python_value(self, value):
if value is not None:
try:
return json.loads(value)
except (TypeError, ValueError):
return value
def db_value(self, value):
if value is not None:
if isinstance(value, ndarray):
value = list(value)
if not isinstance(value, pw.Node):
value = json.dumps(value)
return value
def _e(op):
def inner(self, rhs):
if isinstance(rhs, (list, dict)):
rhs = pw.Value(rhs, converter = self.db_value, unpack = False)
return pw.Expression(self, op, rhs)
return inner
__eq__ = _e(pw.OP.EQ)
__ne__ = _e(pw.OP.NE)
__gt__ = _e(pw.OP.GT)
__ge__ = _e(pw.OP.GTE)
__lt__ = _e(pw.OP.LT)
__le__ = _e(pw.OP.LTE)
__hash__ = pw.Field.__hash__
def __getitem__(self, item):
return JSONPath(self)[item]
def set(self, value, as_json = None):
return JSONPath(self).set(value, as_json)
def update(self, data):
return JSONPath(self).update(data)
def remove(self):
return JSONPath(self).remove()
def json_type(self):
return pw.fn.json_type(self)
def length(self):
return pw.fn.json_array_length(self)
def children(self):
"""
Schema of `json_each` and `json_tree`:
key,
value,
type TEXT (object, array, string, etc),
atom (value for primitive/scalar types, NULL for array and object)
id INTEGER (unique identifier for element)
parent INTEGER (unique identifier of parent element or NULL)
fullkey TEXT (full path describing element)
path TEXT (path to the container of the current element)
json JSON hidden (1st input parameter to function)
root TEXT hidden (2nd input parameter, path at which to start)
"""
return pw.fn.json_each(self)
def tree(self):
return pw.fn.json_tree(self)

View File

@ -0,0 +1,15 @@
# -*- coding: utf-8 -*-
from . import layouts
from .layouts import app
app.layout = layouts.base
def run(*args, **kwargs):
app.run_server(*args, **kwargs)
if __name__ == "__main__":
run(debug = True)

8
anisotropy/gui/app.py Normal file
View File

@ -0,0 +1,8 @@
# -*- coding: utf-8 -*-
import dash
import dash_bootstrap_components as dbc
app = dash.Dash(__name__, external_stylesheets = [ dbc.themes.LUX ])
app.title = "Anisotropy"
app.config["update_title"] = None

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.3 KiB

View File

Before

Width:  |  Height:  |  Size: 58 KiB

After

Width:  |  Height:  |  Size: 58 KiB

View File

@ -0,0 +1,29 @@
# -*- coding: utf-8 -*-
from . import (
runner,
settings,
database,
visualization,
about,
base
)
from .base import app
runner = runner.layout
settings = settings.layout
database = database.layout
visualization = visualization.layout
about = about.layout
base = base.layout
__all__ = [
"runner",
"settings",
"database",
"visualization",
"about",
"base",
"app"
]

View File

@ -0,0 +1,21 @@
# -*- coding: utf-8 -*-
from dash import html
import anisotropy
###
# Layout
##
layout = html.Div([
html.H1(anisotropy.__name__),
html.Hr(),
html.P([
"Author: {}".format(anisotropy.__author__),
html.Br(),
"License: {}".format(anisotropy.__license__),
html.Br(),
"Version: {}".format(anisotropy.__version__)
])
])

View File

@ -0,0 +1,81 @@
# -*- coding: utf-8 -*-
from dash import html
from dash import dcc
from dash.dependencies import Input, Output
import dash_bootstrap_components as dbc
from . import (
runner,
settings,
database,
visualization,
about
)
from ..app import app
from .. import styles
import anisotropy
###
# Layout
##
layout = html.Div([
# Location
dcc.Location(id = "url", refresh = False),
# Sidebar
html.Div([
# Sidebar
html.H2([html.Img(src = "/assets/simple.png", height = "150px")], style = styles.logo),
html.Hr(style = { "color": "#ffffff" }),
dbc.Nav([
dbc.NavLink("Runner", href = "/", active = "exact", style = styles.white),
dbc.NavLink("Settings", href = "/settings", active = "exact", style = styles.white),
dbc.NavLink("Database", href = "/database", active = "exact", style = styles.white),
dbc.NavLink("Visualization", href = "/visualization", active = "exact", style = styles.white),
dbc.NavLink("About", href = "/about", active = "exact", style = styles.white),
], vertical = True, pills = True),
# Misc
html.Hr(style = styles.white),
dbc.Container([
dbc.Row([
dbc.Col("v1.2.0"),
dbc.Col(
html.A(
html.Img(src = "/assets/gh-light.png", height = "20px"),
href = anisotropy.__repository__
)
)
])
], style = styles.misc)
], style = styles.sidebar),
# Content
html.Div(id = "page-content", style = styles.page),
], style = styles.content)
###
# Callbacks
##
@app.callback(
Output("page-content", "children"),
[ Input("url", "pathname") ]
)
def displayPage(pathname):
if pathname == "/settings":
return settings.layout
elif pathname == "/about":
return about.layout
elif pathname == "/database":
return database.layout
elif pathname == "/visualization":
return visualization.layout
else:
return runner.layout

View File

@ -0,0 +1,101 @@
# -*- coding: utf-8 -*-
from dash.dash_table import DataTable
from dash import html
from dash import dcc
import dash_bootstrap_components as dbc
from dash.dependencies import Input, Output, State
import pathlib
from os import environ
from ..app import app
from .. import styles
###
# Layout
##
layout = html.Div([
# Messages and timer
dbc.Alert(
id = "db_status",
duration = 10000,
dismissable = True,
is_open = False,
style = styles.message
),
# dcc.Interval(id = "interval", interval = 1000, n_intervals = 0),
# Query
html.H2("Database"),
html.Hr(),
html.P("Query"),
dcc.Textarea(id = "db_input", style = { "min-width": "100%"}),
html.Br(),
dbc.Button("Query", id = "query", style = styles.minWidth),
# Output
html.Hr(),
html.P("Output"),
DataTable(
id = "db_output",
columns = [],
data = [],
style_table = { "overflow": "scroll"},
style_cell = {
'textAlign': 'left',
'width': '150px',
'minWidth': '180px',
'maxWidth': '180px',
'whiteSpace': 'no-wrap',
'overflow': 'hidden',
'textOverflow': 'ellipsis',
}
),
])
###
# Callbacks
##
@app.callback(
Output("db_output", "columns"),
Output("db_output", "data"),
Output("db_status", "children"),
Output("db_status", "is_open"),
Output("db_status", "color"),
[ Input("query", "n_clicks") ],
[ State("db_input", "value") ],
prevent_initial_call = True
)
def db_query(clicks, db_input):
from anisotropy import database
import peewee as pw
path = pathlib.Path(environ["AP_CWD"], environ["AP_DB_FILE"])
db = database.Database(path)
try:
db.connect()
cursor = db.execute_sql(db_input)
except pw.OperationalError as e:
db.close()
return None, None, str(e), True, "danger"
else:
columns = [ { "name": col[0], "id": col[0] } for col in cursor.description ]
data = []
for row in cursor.fetchall():
item = {}
for key, value in zip(columns, row):
item[key["name"]] = value
data.append(item)
db.close()
return columns, data, "", False, "success"

View File

@ -0,0 +1,180 @@
# -*- coding: utf-8 -*-
from dash.dash_table import DataTable
from dash import html
from dash import dcc
import dash_bootstrap_components as dbc
from dash.dependencies import Input, Output, State
import pathlib
import os
from os import environ
from ..app import app
from .. import styles
from .. import utils
###
# Layout
##
layout = html.Div([
# Messages and timer
dbc.Alert(
id = "status",
duration = 10000,
dismissable = True,
is_open = False,
style = styles.message
),
dcc.Interval(id = "interval", interval = 1000, n_intervals = 0),
# Runner
html.H2("Runner"),
html.Hr(),
html.P("Execution (leave zero for the latest)"),
dcc.Input(id = "execution", type = "number", value = 0, min = 0, style = styles.minWidth),
html.Br(),
dbc.Button("Start", id = "start", color = "success", style = styles.minWidth),
dbc.Button("Stop", id = "stop", color = "danger", disabled = True, style = styles.minWidth),
# Monitor
html.H2("Monitor"),
html.Hr(),
html.P(id = "runner-status"),
DataTable(id = "monitor", columns = [], data = [], style_table = styles.table),
# Log
html.H2("Log"),
html.Hr(),
dbc.Button("Delete", id = "delete", style = styles.minWidth),
dcc.Textarea(id = "logger", disabled = True, style = styles.bigText)
])
###
# Callbacks
##
@app.callback(
Output("start", "active"),
[ Input("start", "n_clicks") ],
[ State("execution", "value") ],
prevent_initial_call = True
)
def runnerStart(clicks, execution):
import subprocess
command = [
"anisotropy",
"compute",
"-v",
"--path", environ["AP_CWD"],
"--conf", environ["AP_CONF_FILE"],
"--pid", "anisotropy.pid",
"--logfile", environ["AP_LOG_FILE"],
]
if execution > 0:
command.extend([ "--exec-id", str(execution) ])
subprocess.run(
command,
start_new_session = True,
)
return True
@app.callback(
Output("stop", "active"),
[ Input("stop", "n_clicks") ],
prevent_initial_call = True
)
def runnerStop(clicks):
import psutil
import signal
pidpath = pathlib.Path(environ["AP_CWD"], "anisotropy.pid")
try:
pid = int(open(pidpath, "r").read())
master = psutil.Process(pid)
except (FileNotFoundError, psutil.NoSuchProcess):
return True
else:
os.killpg(master.pid, signal.SIGTERM)
return True
@app.callback(
Output("monitor", "columns"),
Output("monitor", "data"),
Output("runner-status", "children"),
Output("start", "disabled"),
Output("stop", "disabled"),
Output("delete", "disabled"),
[ Input("interval", "n_intervals") ],
)
def monitorUpdate(intervals):
import psutil
pidpath = pathlib.Path(environ["AP_CWD"], "anisotropy.pid")
processes = []
try:
pid = int(open(pidpath, "r").read())
master = psutil.Process(pid)
except (FileNotFoundError, psutil.NoSuchProcess):
return [], [], "Status: not running", False, True, False
else:
for process in [ master, *master.children() ]:
created = psutil.time.localtime(process.create_time())
processes.append({
"name": process.name(),
"pid": process.pid,
"status": process.status(),
"memory": utils.getSize(process.memory_full_info().uss),
"threads": process.num_threads(),
"created": "{}:{}:{}".format(created.tm_hour, created.tm_min, created.tm_sec)
})
columns = [ { "name": col, "id": col } for col in processes[0].keys() ]
return columns, processes, "Status: running", True, False, True
@app.callback(
Output("logger", "value"),
[ Input("interval", "n_intervals") ]
)
def logUpdate(intervals):
logpath = pathlib.Path(environ["AP_CWD"], "anisotropy.log")
if os.path.exists(logpath):
with open(logpath, "r") as io:
log = io.read()
return log
else:
return "Not found"
@app.callback(
Output("delete", "active"),
[ Input("delete", "n_clicks") ],
prevent_initial_call = True
)
def logDelete(clicks):
logpath = pathlib.Path(environ["AP_CWD"], "anisotropy.log")
if os.path.exists(logpath):
os.remove(logpath)
return True

View File

@ -0,0 +1,141 @@
# -*- coding: utf-8 -*-
from dash import html
from dash import dcc
import dash_bootstrap_components as dbc
from dash.dependencies import Input, Output, State
import pathlib
from os import environ
from ..app import app
from .. import styles
###
# Layout
##
layout = html.Div([
# Messages
dbc.Alert(
id = "status",
duration = 10000,
dismissable = True,
is_open = False,
style = styles.message
),
dbc.Alert(
id = "general-status",
duration = 10000,
dismissable = True,
is_open = False,
style = styles.message
),
# General
html.H2("General"),
html.Hr(),
html.P("Path"),
dcc.Input(id = "cwd", style = { "min-width": "500px" }),
html.Br(),
dbc.Button("Save general", id = "general-save", style = styles.minWidth),
# Options
html.H2("Options"),
html.Hr(),
html.P("Nprocs"),
dcc.Input(id = "nprocs", type = "number", style = styles.minWidth),
html.P("Stage"),
dcc.Dropdown(
id = "stage",
options = [ { "label": k, "value": k } for k in ["all", "shape", "mesh", "flow", "postProcess"] ],
style = styles.minWidth
),
dbc.Button("Save", id = "submit", style = styles.minWidth),
# Cases
html.H2("Cases"),
html.Hr(),
dcc.Textarea(id = "cases", style = styles.bigText),
])
###
# Callbacks
##
@app.callback(
Output("general-status", "children"),
Output("general-status", "is_open"),
Output("general-status", "color"),
[ Input("general-save", "n_clicks") ],
[
State("cwd", "value"),
],
prevent_initial_call = True
)
def generalSave(clicks, cwd):
path = pathlib.Path(cwd)
if not path.is_absolute():
return "Cwd path must be absolute", True, "danger"
environ["AP_CWD"] = str(path)
return "General settings saved", True, "success"
@app.callback(
Output("cwd", "value"),
Output("nprocs", "value"),
Output("stage", "value"),
Output("cases", "value"),
[ Input("url", "pathname") ]
)
def settingsLoad(pathname):
from anisotropy.core import config as core_config
import toml
path = pathlib.Path(environ["AP_CWD"], environ["AP_CONF_FILE"])
config = core_config.default_config()
if path.exists():
config.load(path)
return environ["AP_CWD"], config["nprocs"], config["stage"], toml.dumps(config.content)
@app.callback(
Output("status", "children"),
Output("status", "is_open"),
Output("status", "color"),
[ Input("submit", "n_clicks") ],
[
State("nprocs", "value"),
State("stage", "value"),
State("cases", "value")
],
prevent_initial_call = True
)
def settingsSave(nclick, nprocs, stage, cases):
from anisotropy.core import config as core_config
import toml
path = pathlib.Path(environ["AP_CWD"], environ["AP_CONF_FILE"])
config = core_config.default_config()
if path.exists():
config.load(path)
config.update(
nprocs = nprocs,
stage = stage
)
try:
config.content = toml.loads(cases)
config.dump(path)
except Exception as e:
return str(e), True, "danger"
else:
return f"Saved to { path }", True, "success"

View File

@ -0,0 +1,405 @@
# -*- coding: utf-8 -*-
from dash import html
from dash import dcc
import dash_bootstrap_components as dbc
from dash.dependencies import Input, Output, State
import dash_vtk
import dash_vtk.utils
import vtk
import pathlib
from os import environ
from ..app import app
from .. import styles
class MeshRepresentation(object):
def __init__(self, path):
self.mesh = vtk.vtkXMLUnstructuredGridReader()
self.mesh.SetFileName(path)
self.mesh.Update()
self.clip = False
self.crinkle = False
self.wireframe = True
self.planeNormal = [0, 0, -1]
self.planeOrigin = [0, 0, 1]
def getRepresentation(self):
if not self.clip:
dataset = self.mesh.GetOutput()
else:
plane = vtk.vtkPlane()
plane.SetNormal(self.planeNormal)
plane.SetOrigin(self.planeOrigin)
if self.crinkle:
implfunc = vtk.vtkImplicitBoolean()
implfunc.AddFunction(plane)
implfunc.Modified()
clipper = vtk.vtkExtractGeometry()
clipper.SetInputConnection(self.mesh.GetOutputPort())
clipper.SetImplicitFunction(implfunc)
clipper.SetExtractInside(1)
clipper.SetExtractOnlyBoundaryCells(0)
clipper.SetExtractBoundaryCells(1)
clipper.ExtractInsideOff()
else:
clipper = vtk.vtkTableBasedClipDataSet()
clipper.SetInputData(self.mesh.GetOutput())
clipper.SetClipFunction(plane)
clipper.SetValue(0.0)
clipper.GenerateClippedOutputOn()
clipper.Update()
dataset = clipper.GetOutput()
representation = dash_vtk.GeometryRepresentation(
id = "mesh-representation",
children = [
dash_vtk.Mesh(id = "mesh", state = dash_vtk.utils.to_mesh_state(dataset))
],
property = { "edgeVisibility": self.wireframe }
)
return representation
def databaseColumns():
from anisotropy.database import Database
import re
db = Database()
idcol = re.compile(r"\s*_id")
columns = []
for table in db.tables:
for column in table._meta.columns.keys():
if not idcol.search(column):
columns.append(column)
return columns
###
# Layout
##
controls = html.Div([
html.P("Execution"),
dcc.Input(id = "exec_id", type = "number", value = 0, min = 0, style = styles.minWidth),
html.Br(),
html.P("Structure"),
dcc.Dropdown(
id = "structure",
options = [ { "label": v, "value": v } for v in [ "simple", "bodyCentered", "faceCentered" ] ],
value = "simple",
style = styles.minWidth
),
html.Br(),
html.P("Direction"),
dcc.Dropdown(
id = "direction",
options = [ { "label": str(v), "value": str(v)} for v in [ [1., 0., 0.], [0., 0., 1.], [1., 1., 1.] ] ],
value = str([1., 0., 0.]),
),
html.Br(),
html.P("Alpha"),
dcc.Input(id = "alpha", type = "number", value = 0.01, min = 0.01, step = 0.01, style = styles.minWidth),
html.Br(),
html.P("Clip"),
dbc.Checkbox(id = "clip", value = False),
html.Br(),
html.P("Crinkle"),
dbc.Checkbox(id = "crinkle", value = False),
html.Br(),
html.P("Plane normal"),
dcc.Input(id = "normal-x", type = "number", value = 0, style = { "width": "100px", "margin-bottom": "20px" }),
dcc.Input(id = "normal-y", type = "number", value = 0, style = { "width": "100px", "margin-bottom": "20px" }),
dcc.Input(id = "normal-z", type = "number", value = -1, style = { "width": "100px", "margin-bottom": "20px" }),
html.Br(),
html.P("Plane origin"),
dcc.Input(id = "origin-x", type = "number", value = 0, style = { "width": "100px", "margin-bottom": "20px" }),
dcc.Input(id = "origin-y", type = "number", value = 0, style = { "width": "100px", "margin-bottom": "20px" }),
dcc.Input(id = "origin-z", type = "number", value = 1, style = { "width": "100px", "margin-bottom": "20px" }),
html.Br(),
html.P("Wireframe"),
dbc.Checkbox(id = "wireframe", value = True),
html.Br(),
dbc.Button("Draw", id = "draw"),
])
plotcontrols = html.Div([
html.P("Execution"),
dcc.Input(id = "plot-exec_id", type = "number", value = 0, min = 0, style = styles.minWidth),
html.Br(),
html.P("Structure"),
dcc.Dropdown(
id = "plot-structure",
options = [
{ "label": v, "value": v }
for v in [ "simple", "bodyCentered", "faceCentered" ]
],
value = "simple"
),
html.Br(),
html.P("Direction"),
dcc.Dropdown(
id = "plot-direction",
options = [
{ "label": str(v), "value": str(v)}
for v in [ [1., 0., 0.], [0., 0., 1.], [1., 1., 1.], "all" ]
],
value = str([1., 0., 0.]),
),
html.Br(),
html.P("Data"),
dcc.Dropdown(
id = "plot-data",
options = [ { "label": v, "value": v } for v in databaseColumns() ],
value = "porosity",
),
html.Br(),
dbc.Button("Draw", id = "plot-draw"),
])
layout = html.Div([
html.H2("Plot"),
html.Hr(),
dbc.Container(
fluid = True,
children = [
dbc.Row([
dbc.Col(width = 4, children = plotcontrols, style = { "min-width": "350px" }),
dbc.Col(
width = 8,
children = [
html.Div(id = "plot-output", style = { "width": "100%", "min-width": "800px" })
],
style = { "min-width": "800px" }
),
], style = { "height": "100%"}),
]),
html.Br(),
html.H2("Mesh"),
html.Hr(),
dbc.Container(
fluid = True,
children = [
dbc.Row([
dbc.Col(width = 4, children = controls, style = { "min-width": "350px" }),
dbc.Col(
width = 8,
children = [
html.Div(
id = "vtk-output",
style = { "height": "800px", "width": "100%", "min-width": "800px" }
)
],
style = { "min-width": "800px" }),
], style = { "height": "100%"}),
])
])
@app.callback(
[ Output("plot-output", "children") ],
[ Input("plot-draw", "n_clicks") ],
[
State("plot-exec_id", "value"),
State("plot-structure", "value"),
State("plot-direction", "value"),
State("plot-data", "value"),
]
)
def plotDraw(clicks, execution, structure, direction, data):
import peewee as pw
from anisotropy.database import Database, tables
import json
from pandas import DataFrame
import plotly.express as px
path = pathlib.Path(environ["AP_CWD"], environ["AP_DB_FILE"])
if not path.is_file():
return [ "Database not found" ]
db = Database(path)
if not db.getExecution(execution):
return [ "Execution not found" ]
for model in db.tables:
try:
column = getattr(model, data)
except AttributeError:
pass
else:
break
if direction == "all":
select = (tables.Shape.alpha, column, tables.Shape.direction)
else:
select = (tables.Shape.alpha, column)
"""query = (
tables.Shape
.select(*select)
.join(tables.Execution, JOIN.LEFT_OUTER)
.switch(tables.Shape)
.join(tables.Mesh, JOIN.LEFT_OUTER)
#.switch(tables.Shape)
.join(tables.FlowOnephase, JOIN.LEFT_OUTER)
.switch(tables.Shape)
.where(
tables.Shape.exec_id == execution,
tables.Shape.label == structure,
)
)"""
query = model.select(*select)
idn = db.tables.index(model)
for table in reversed(db.tables[ :idn]):
query = query.join(table, pw.JOIN.LEFT_OUTER)
query = query.switch(tables.Shape)
query = query.where(
tables.Shape.exec_id == execution,
tables.Shape.label == structure,
)
query = query.order_by(tables.Shape.alpha)
if not direction == "all":
query = query.where(tables.Shape.direction == json.loads(direction))
with db:
if query.exists():
table = []
for row in query.dicts():
for k in row.keys():
if type(row[k]) == list:
row[k] = str(row[k])
table.append(row)
else:
table = None
if not table:
return [ "Results not found" ]
import plotly.io as plt_io
plt_io.templates["custom_dark"] = plt_io.templates["plotly_dark"]
plt_io.templates["custom_dark"]['layout']['paper_bgcolor'] = '#30404D'
plt_io.templates["custom_dark"]['layout']['plot_bgcolor'] = '#30404D'
plt_io.templates['custom_dark']['layout']['yaxis']['gridcolor'] = '#4f687d'
plt_io.templates['custom_dark']['layout']['xaxis']['gridcolor'] = '#4f687d'
fig = px.line(
DataFrame(table), x = "alpha", y = data, title = structure, markers = True
)
if direction == "all":
fig = px.line(
DataFrame(table), x = "alpha", y = data, title = structure, markers = True, color = "direction"
)
fig.layout.template = "custom_dark"
fig.update_xaxes(showline = True, linewidth = 1, linecolor = '#4f687d', mirror = True)
fig.update_yaxes(showline = True, linewidth = 1, linecolor = '#4f687d', mirror = True)
plot = dcc.Graph(
figure = fig
)
return [ plot ]
@app.callback(
[ Output("alpha", "min"), Output("alpha", "max") ],
[ Input("structure", "value") ]
)
def alphaLimits(label):
if label == "simple":
return 0.01, 0.29
elif label == "bodyCentered":
return 0.01, 0.18
elif label == "faceCentered":
return 0.01, 0.13
@app.callback(
[ Output("vtk-output", "children") ],
[ Input("draw", "n_clicks") ],
[
State("exec_id", "value"),
State("structure", "value"),
State("direction", "value"),
State("alpha", "value"),
State("clip", "value"),
State("crinkle", "value"),
State("wireframe", "value"),
State("normal-x", "value"),
State("normal-y", "value"),
State("normal-z", "value"),
State("origin-x", "value"),
State("origin-y", "value"),
State("origin-z", "value"),
],
prevent_initial_call = True
)
def meshDraw(clicks, execution, structure, direction, alpha, clip, crinkle, wireframe, normal_x, normal_y, normal_z, origin_x, origin_y, origin_z):
import meshio
path = pathlib.Path(environ["AP_CWD"], environ["AP_BUILD_DIR"])
path /= "execution-{}".format(execution)
path /= "{}-{}-{}".format(
structure,
direction.replace(" ", ""),
alpha
)
basemeshpath = path / "mesh.msh"
meshpath = path / "mesh.vtu"
if not basemeshpath.exists():
return [ "Mesh not found" ]
if not meshpath.exists() or not meshpath.is_file():
meshold = meshio.read(basemeshpath)
meshold.write(meshpath)
meshrepr = MeshRepresentation(meshpath)
meshrepr.clip = clip
meshrepr.crinkle = crinkle
meshrepr.wireframe = wireframe
meshrepr.planeNormal = [ normal_x, normal_y, normal_z ]
meshrepr.planeOrigin = [ origin_x, origin_y, origin_z ]
view = dash_vtk.View([ meshrepr.getRepresentation() ])
return [ view ]

57
anisotropy/gui/styles.py Normal file
View File

@ -0,0 +1,57 @@
# -*- coding: utf-8 -*-
minWidth = {
"min-width": "200px",
"max-width": "200px",
"margin-bottom": "15px"
}
bigText = {
"min-width": "100%",
"min-height": "500px"
}
message = {
"position": "fixed",
"top": "30px",
"right": "30px",
}
page = {
"padding": "30px 50px 0px 50px"
}
table = {
"margin-bottom": "15px",
"min-width": "100%"
}
sidebar = {
"position": "fixed",
"top": 0,
"left": 0,
"bottom": 0,
"width": "16rem",
"padding": "2rem 1rem",
"background-color": "#363636"
}
content = {
"margin-left": "18rem",
"margin-right": "2rem",
"padding": "2rem 1rem",
}
misc = {
"color": "#757575",
"text-align": "center"
}
logo = {
"color": "#ffffff",
"text-align": "center"
}
white = {
"color": "#ffffff"
}

8
anisotropy/gui/utils.py Normal file
View File

@ -0,0 +1,8 @@
# -*- coding: utf-8 -*-
def getSize(bytes):
for unit in ["", "K", "M", "G", "T", "P"]:
if bytes < 1024:
return f"{bytes:.2f} {unit}iB"
bytes /= 1024

View File

@ -0,0 +1,16 @@
# -*- coding: utf-8 -*-
from . import utils
from . import conversion
from . import metrics
from .mesh import Mesh, MeshingParameters
__all__ = [
"utils",
"conversion",
"metrics",
"Mesh",
"MeshingParameters"
]

View File

@ -0,0 +1,145 @@
from numpy import ndarray
from os import PathLike
import numpy as np
from pathlib import Path
from . import utils
def write_neutral(
points: ndarray,
cells: list[utils.CellBlock],
dim: int,
filename: PathLike
):
"""Write mesh to the netgen file (Neutral format [.mesh]).
:param points:
Array of points.
:param cells:
List of cell blocks.
:param dim:
Dimension of the mesh.
:param filename:
Path of the file.
"""
precision = "6f"
floatValue = "{value:>{width}." + precision + "}"
intValue = "{value:>{width}}"
outfile = open(Path(filename).resolve(), "w")
# write points
outfile.write(str(len(points)) + "\n")
for n in range(len(points)):
point = points[n]
outfile.write(floatValue.format(value = point[0], width = 10) + " ")
outfile.write(floatValue.format(value = point[1], width = 9) + " ")
if dim == 3:
outfile.write(floatValue.format(value = point[2], width = 9) + " ")
outfile.write("\n")
# write volume cells
if dim == 3:
count = sum([ len(cell.data) for cell in cells if cell.dim == 3])
outfile.write(str(count) + "\n")
for cellBlock in cells:
if cellBlock.dim == 3:
for cell in cellBlock.data:
outfile.write(intValue.format(value = 1, width = 4))
for pointId in cell:
outfile.write(" ")
# shift id up, in netgen it starts from one
outfile.write(intValue.format(value = pointId + 1, width = 8))
outfile.write("\n")
# write faces
count = sum([ len(cell.data) for cell in cells if cell.dim == 2])
outfile.write(str(count) + "\n")
for cellBlock in cells:
if cellBlock.dim == 2:
for index, cell in zip(cellBlock.indices, cellBlock.data):
outfile.write(intValue.format(value = index, width = 4) + " ")
for pointId in cell:
outfile.write(" ")
# shift id up, in netgen it starts from one
outfile.write(intValue.format(value = pointId + 1, width = 8))
outfile.write("\n")
# write segments
# important?
def read_neutral(filename: PathLike) -> tuple[list, list]:
"""Read mesh from netgen file (Neutral format [.mesh]).
:param filename:
Path of the file.
:return:
List of points and list of cell blocks.
"""
infile = open(Path(filename).resolve(), "r")
# read points from file, starts with single integer
npoints = int(infile.readline())
points = []
for n in range(npoints):
points.append(np.fromstring(infile.readline(), dtype = float, sep = " "))
# dimension
dim = None if len(points) == 0 else points[0].size
# read volume cells
nvolumes = int(infile.readline())
volume_indices = []
volumes = []
for n in range(nvolumes):
data = np.fromstring(infile.readline(), dtype = int, sep = " ")
volume_indices.append(data[0])
# shift node indices down, in netgen it starts from one, need from zero
volumes.append(data[1: ] - 1)
# read surface cells
nsurfaces = int(infile.readline())
surface_indices = []
surfaces = []
for n in range(nsurfaces):
data = np.fromstring(infile.readline(), dtype = int, sep = " ")
surface_indices.append(data[0])
# shift node indices down, in netgen it starts from one, need from zero
surfaces.append(data[1: ] - 1)
# read segment cells
# important?
# write data to object
points = np.asarray(points)
cells = []
if len(volumes) > 0:
cells += utils.collect_cells(volumes, 3, volume_indices)
if len(surfaces) > 0:
cellBlocks = utils.collect_cells(surfaces, 2, surface_indices)
if dim == 3:
for cellBlock in cellBlocks:
cellBlock.is_boundary = True
cells += cellBlocks
return points, cells

242
anisotropy/meshing/mesh.py Normal file
View File

@ -0,0 +1,242 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from __future__ import annotations
from numpy import ndarray
from os import PathLike
import numpy as np
import netgen.occ as ng_occ
import netgen.meshing as ng_meshing
import meshio
import pathlib
from . import utils
from . import metrics
from . import conversion
class NoGeometrySpecified(Exception):
pass
class NotSupportedMeshFormat(Exception):
"""Exception for not supported mesh format IO operations.
"""
pass
class MeshingParameters:
def __init__(self, **kwargs):
"""A MeshingParameters object contains parameters for netgen mesher.
"""
self.preset(2, **kwargs)
def preset(self, key: int, **kwargs):
"""Apply predefined parameters.
:param key:
0: very_coarse
1: coarse
2: moderate
3: fine
4: very_fine
"""
self.__dict__.update(kwargs)
self.maxh = kwargs.get("maxh", 0.1)
self.curvaturesafety = kwargs.get("curvaturesafety", [1, 1.5, 2, 3, 5][key])
self.segmentsperedge = kwargs.get("segmentsperedge", [0.3, 0.5, 1, 2, 3][key])
self.grading = kwargs.get("grading", [0.7, 0.5, 0.3, 0.2, 0.1][key])
self.chartdistfac = kwargs.get("chartdistfac", [0.8, 1, 1.5, 2, 5][key])
self.linelengthfac = kwargs.get("linelengthfac", [0.2, 0.35, 0.5, 1.5, 3][key])
self.closeedgefac = kwargs.get("closeedgefac", [0.5, 1, 2, 3.5, 5][key])
self.minedgelen = kwargs.get("minedgelen", [0.002, 0.02, 0.2, 1.0, 2.0][key])
self.surfmeshcurvfac = kwargs.get("surfmeshcurvfac", [1, 1.5, 2, 3, 5.0][key])
self.optsteps2d = kwargs.get("optsteps2d", 5)
self.optsteps3d = kwargs.get("optsetps3d", 5)
return self
def get(self) -> ng_meshing.MeshingParameters:
return ng_meshing.MeshingParameters(**self.__dict__)
def __repr__(self):
items = [ f"{ key }: { val }" for key, val in self.__dict__.items() ]
return "<MeshingParameters, " + ", ".join(items) + ">"
class Mesh:
def __init__(self, shape = None):
"""A Mesh object contains a mesh
:param shape:
OCC shape.
"""
self.shape = shape
self.mesh = None
self.points = []
self.cells = []
self._dim = None
@property
def dim(self) -> float | None:
"""Mesh dimension.
:return:
float or None if mesh contains no points.
"""
return None if len(self.points) == 0 else self.points[0].size
@property
def geometry(self) -> ng_occ.OCCGeometry | None:
"""Netgen OCCGeometry object.
:return:
OCCGeometry object that can be used for meshing
or None if shape is not defined.
"""
if self.shape is not None:
if hasattr(self.shape, "geometry"):
return self.shape.geometry
else:
return ng_occ.OCCGeometry(self.shape)
return None
def generate(
self,
parameters: ng_meshing.MeshingParameters = None,
refineSteps: int = 0,
refineOptimize: bool = True,
scale: float = 0
):
"""Generate mesh using netgen mesher on occ geometry.
:param parameters:
Meshing parameters object.
:param refineSteps:
Refine mesh after main meshing.
:param refineOptimize:
If True, optimize surface and volume mesh after each refine step.
:param scale:
Scale mesh after all operations.
"""
if not self.geometry:
raise NoGeometrySpecified("Cannot build mesh without geometry")
parameters = parameters or MeshingParameters()
# generate volume mesh
mesh = self.geometry.GenerateMesh(parameters.get())
if refineSteps > 0:
for n in range(refineSteps):
mesh.Refine()
# optimize surface and volume mesh
if refineOptimize:
mesh.OptimizeMesh2d(parameters.get())
mesh.OptimizeVolumeMesh(parameters.get())
if scale > 0:
mesh.Scale(scale)
self.points = utils.extract_netgen_points(mesh)
self.cells = []
if mesh.dim == 3:
volume_cells, volume_indices = utils.extract_netgen_cells(mesh, 3)
self.cells += utils.collect_cells(volume_cells, 3, volume_indices)
surface_cells, surface_indices = utils.extract_netgen_cells(mesh, 2)
cell_blocks = utils.collect_cells(surface_cells, 2, surface_indices)
for cellBlock in cell_blocks:
cellBlock.is_boundary = True
self.cells += cell_blocks
# TODO: dim == 2
return self
def read(self, filename: str):
"""Import a mesh from the file.
:param filename:
Path of the file.
"""
path = pathlib.Path(filename).resolve()
ext = path.suffix[1: ]
# Force netgen neutral format
if ext == "mesh":
self.points, self.cells = conversion.read_neutral(path)
else:
mesh = meshio.read(path)
self.points = mesh.points
self.cells = mesh.cells
return self
def write(self, filename: PathLike):
"""Export mesh to the file.
:param filename:
Path of the file to store the given mesh in.
"""
path = pathlib.Path(filename).resolve()
ext = path.suffix[1: ]
# Force netgen neutral format
if ext == "mesh":
conversion.write_neutral(self.points, self.cells, self.dim, path)
else:
mesh = meshio.Mesh(self.points, self.cells)
mesh.write(path)
@property
def volumes(self) -> list[ndarray]:
"""Volumes.
:return:
List of points.
"""
points = []
for cellBlock in self.cells:
if cellBlock.dim == 3:
points += [ *self.points[cellBlock.data] ]
return points
@property
def volume(self) -> float:
"""Volume of whole mesh.
:return:
Volume.
"""
return np.sum([ metrics.volume(cell) for cell in self.volumes ])
@property
def faces(self) -> list[ndarray]:
"""Boundary faces.
:return:
List of boundary faces.
"""
points = []
for cellBlock in self.cells:
if cellBlock.dim == 2 and cellBlock.is_boundary:
points += [ *self.points[cellBlock.data] ]
return points

View File

@ -0,0 +1,43 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import numpy as np
from numpy import ndarray
from . import utils
def volume(points: ndarray) -> ndarray:
"""Volume of the cell.
:param points:
Array of points that represents volumetric cell.
:return:
Volume.
"""
cellType = utils.cell_type(3, len(points))
if cellType == "tetra":
return 1 / 6 * np.linalg.det(np.append(points.transpose(), np.array([[1, 1, 1, 1]]), axis = 0))
else:
raise Exception(f"Not supported cell type '{ cellType }'")
def area(points: ndarray) -> ndarray:
"""Area of the cell.
:param points:
Array of points that represents surface cell.
:return:
Area.
"""
cellType = utils.cell_type(2, len(points))
vecs = np.roll(points, shift = -1, axis = 0) - points
if cellType == "triangle":
return 0.5 * np.abs(np.cross(*vecs[ :-1]))
else:
raise Exception(f"Not supported cell type '{ cellType }'")

178
anisotropy/meshing/utils.py Normal file
View File

@ -0,0 +1,178 @@
from __future__ import annotations
from numpy import ndarray
import numpy as np
import netgen.meshing as ng_meshing
import meshio
topo_dim = meshio._mesh.topological_dimension
cell_nodes = meshio._common.num_nodes_per_cell
def cell_type(dimension: int, num_nodes: int) -> str | None:
"""Detect a cell topology type.
:param dimension:
Cell dimension.
:param num_nodes:
Number of nodes in the cell.
:return:
Cell type or None if cell type if not supported.
"""
for dim in topo_dim.keys():
for num in cell_nodes.keys():
if (
topo_dim[dim] == dimension and
cell_nodes[num] == num_nodes and
dim == num
):
return dim
return None
class CellBlock:
def __init__(
self,
cellType: str,
data: list | ndarray,
tags: list[str] = [],
indices: list | ndarray = [],
is_boundary: bool = False,
names: list[str] = []
):
"""A CellBlock object contains cells information of the same type.
:param cellType:
Type of the cell.
:param data:
Array of cells.
:param tags:
Some cell tags.
:param indices:
Array of indices of cells. Must be the same length as data length.
Usefull for detecting boundary faces, etc.
:param is_boundary:
Flag that cells are not internal.
:param names:
Array of names of cells. Must be the same length as data length.
Usefull for detecting boundary faces, etc.
"""
self.type = cellType
self.data = data
if cellType.startswith("polyhedron"):
self.dim = 3
else:
self.data = np.asarray(self.data)
self.dim = topo_dim[cellType]
self.tags = tags
self.indices = np.asarray(indices)
self.is_boundary = is_boundary
self.names = names
def __repr__(self):
items = [
"CellBlock",
f"type: { self.type }",
f"num cells: { len(self.data) }",
f"tags: { self.tags }",
]
return "<" + ", ".join(items) + ">"
def __len__(self):
return len(self.data)
def collect_cells(cells: list, dim: int, indices: list | ndarray = None, cellType: str = None) -> list[CellBlock]:
"""Collect cell blocks from a list of cells.
:param cells:
List of cells.
:param dim:
Cell dimension.
:param indices:
List of cell indices. Must be the same length as list of cells.
:return:
List of cell blocks.
"""
cellBlocks = []
if indices is not None:
assert len(cells) == len(indices), "cells and indices arrays must be the same length"
for n, cell in enumerate(cells):
cellType = cellType or cell_type(dim, len(cell))
index = indices[n] if indices else None
is_added = False
# update cell block
for cellBlock in cellBlocks:
if cellBlock.type == cellType:
cellBlock.data += [ cell ]
if index:
cellBlock.indices.append(index)
is_added = True
# or create new
if not is_added:
cellBlock = CellBlock(cellType, [])
cellBlock.data = [ cell ]
if index:
cellBlock.indices = [ index ]
cellBlocks.append(cellBlock)
# convert arrays to numpy
for cellBlock in cellBlocks:
cellBlock.data = np.asarray(cellBlock.data)
cellBlock.indices = np.asarray(cellBlock.indices)
return cellBlocks
def extract_netgen_points(mesh: ng_meshing.Mesh) -> ndarray:
"""Extract points from netgen mesh object.
:param mesh:
Netgen mesh object.
:return:
Array of points.
"""
return np.array([ point.p for point in mesh.Points() ], dtype = float)
def extract_netgen_cells(mesh: ng_meshing.Mesh, dim: int) -> tuple[list, list]:
"""Extract cells from netgen mesh according to cell dimension.
:param mesh:
Netgen mesh object.
:param dim:
Cell dimension.
:return:
List of cells and list of indices.
"""
cells = []
indices = []
elements = {
0: mesh.Elements0D(),
1: mesh.Elements1D(),
2: mesh.Elements2D(),
3: mesh.Elements3D()
}[dim]
if len(elements) > 0:
for cell in elements:
# shift nodes values, they should start from zero
nodes = np.array([ node.nr for node in cell.points ], dtype = int) - 1
cells.append(nodes)
indices.append(cell.index)
return cells, indices

View File

@ -0,0 +1,137 @@
import dearpygui.dearpygui as dpg
def create_geometry(parent):
with dpg.child_window(
tag = "geometry_container",
parent = parent
) as container:
dpg.add_button(label = "Add case", callback = Case.create(container))
with dpg.group(horizontal = True):
dpg.add_text("Cases")
dpg.add_separator()
return container
class Case:
case_list = []
def __init__(self, parent):
self.parent = parent
self.uuid = dpg.generate_uuid()
self.values = {}
Case.case_list.append(self)
def __str__(self):
return f"<Case: { self.uuid }, type: geometry>"
def __repr__(self):
return self.__str__()
@staticmethod
def create(parent):
def create_callback():
case = Case(parent)
with dpg.collapsing_header(label = f"g-case-{ case.uuid }", parent = case.parent) as case._case:
with dpg.child_window(height = 300) as case.child_window:
with dpg.group(horizontal = True):
case._source = dpg.add_combo(["", "from file", "periodic"], label = "Source")
dpg.add_button(label = "Remove", callback = case.remove_case)
dpg.add_button(label = "Edit", callback = case.edit_case)
case.create_stats()
case.create_params_file()
case.create_params_periodic()
return create_callback
def remove_case(self):
dpg.delete_item(self._case)
Case.case_list.remove(self)
def edit_case(self):
source_val = dpg.get_value(self._source)
if source_val == "from file":
dpg.show_item(self._file_params)
elif source_val == "periodic":
dpg.show_item(self._periodic_params)
def create_stats(self):
with dpg.table(header_row = True, parent = self.child_window) as self._case_stats:
dpg.add_table_column(label = "Parameters")
dpg.add_table_column(label = "")
for k, v in self.values.items():
with dpg.table_row():
dpg.add_text(k)
dpg.add_text(v)
def update_stats(self):
dpg.delete_item(self._case_stats)
self.create_stats()
def create_params_file(self):
with dpg.window(show = False, modal = True, height = 600, width = 400) as self._file_params:
with dpg.group(horizontal = True):
self._file_input = dpg.add_input_text(label = "File")
dpg.add_button(label = "..", callback = lambda: dpg.show_item(self._file_dialog))
with dpg.file_dialog(
directory_selector = False,
show = False,
modal = True,
callback = lambda s, d: dpg.set_value(self._file_input, d['file_path_name']),
height = 400,
width = 600
) as self._file_dialog:
dpg.add_file_extension(".geo")
dpg.add_button(label = "Save", width = -1)
def create_params_periodic(self):
def alpha_type_cb(s, d):
type_val = dpg.get_value(self._alpha_type)
if type_val == "float":
dpg.show_item(self._alpha_float)
dpg.hide_item(self._alpha_range)
elif type_val == "range":
dpg.hide_item(self._alpha_float)
dpg.show_item(self._alpha_range)
with dpg.window(show = False, modal = True, height = 600, width = 400) as self._periodic_params:
self._label_str = dpg.add_input_text(label = "label")
with dpg.group(horizontal = True):
self._alpha_type = dpg.add_combo(["float", "range"], default_value = "float", width = 100, callback = alpha_type_cb)
self._alpha_float = dpg.add_input_float(label = "alpha")
self._alpha_range = dpg.add_input_floatx(label = "alpha", size = 3, show = False)
self._r0_float = dpg.add_input_float(label = "r_0")
dpg.add_separator()
dpg.add_button(label = "Save", width = -1, callback = self.save_params)
def save_params(self):
source_val = dpg.get_value(self._source)
if source_val == "from file":
self.values.update(
filepath = dpg.get_value(self._file_input)
)
elif source_val == "periodic":
self.values.update(
label = dpg.get_value(self._label_str),
alpha = dpg.get_value(self._alpha_range) if dpg.get_value(self._alpha_type) == "range" else dpg.get_value(self._alpha_float),
r0 = dpg.get_value(self._r0_float),
)
self.update_stats()
dpg.hide_item(self._periodic_params)

View File

@ -0,0 +1,88 @@
import dearpygui.dearpygui as dpg
import pathlib
import geometry
gui_path = pathlib.Path(__file__).resolve().parent
dpg.create_context()
with dpg.font_registry():
print(gui_path)
default_font = dpg.add_font(gui_path / "fonts/NotoSansMono-Regular-Nerd-Font-Complete.ttf", 20)
with dpg.window(tag = "Primary") as primary:
dpg.bind_font(default_font)
with dpg.menu_bar(tag = "m1"):
with dpg.menu(label = "File"):
dpg.add_menu_item(label = "Preferences", callback = lambda: dpg.configure_item("preferences", show = True))
dpg.add_menu_item(label = "Exit", callback = lambda: dpg.destroy_context())
with dpg.menu(label = "Tools"):
dpg.add_menu_item(label = "Geometry", callback = lambda: dpg.show_item("geometry"))
dpg.add_menu_item(label = "Mesh", callback = lambda: dpg.show_item("mesh"))
dpg.add_menu_item(label = "Solver", callback = lambda: dpg.show_item("solver"))
dpg.add_menu_item(label = "Runner", callback = lambda: dpg.show_item("runner"))
dpg.add_separator()
dpg.add_menu_item(label = "Database", callback = lambda: dpg.show_item("database"))
dpg.add_menu_item(label = "Post-processing", callback = lambda: dpg.show_item("post-processing"))
with dpg.child_window(tag = "Test22", border = False, height = -40):
with dpg.child_window(tag = "Test", border = False):
with dpg.tab_bar(label = "tabbar"):
with dpg.tab(label = "Geometry", tag = "geometry", closable = True, show = False):
geometry.create_geometry("geometry")
with dpg.item_handler_registry(tag="widget handler"):
dpg.add_item_edited_handler(callback = lambda: dpg.delete_item("mesh") if not dpg.is_item_shown("mesh") else None)
with dpg.tab(label = "Mesh", tag = "mesh", closable = True, show = False):
with dpg.child_window(tag = "mesh_container", border = True):
dpg.add_text("")
# bind item handler registry to item
dpg.bind_item_handler_registry("mesh", "widget handler")
with dpg.tab(label = "Solver", tag = "solver", closable = True, show = False):
with dpg.child_window(tag = "solver_container", border = True):
dpg.add_text("")
with dpg.tab(label = "Runner", tag = "runner", closable = True, show = False):
with dpg.child_window(tag = "runner_container", border = True):
dpg.add_text("")
with dpg.tab(label = "Database", tag = "database", closable = True, show = False):
with dpg.child_window(tag = "database_container", border = True):
dpg.add_text("")
with dpg.tab(label = "Post-processing", tag = "post-processing", closable = True, show = False):
with dpg.child_window(tag = "post-processing_container", border = True):
dpg.add_text("")
dpg.add_separator()
with dpg.table(tag = "statusbar", header_row = False):
dpg.add_table_column(label = "")
dpg.add_table_column(label = "")
with dpg.table_row():
dpg.add_button(label = "Debug", callback = lambda: print(dpg.get_value("mesh")))
dpg.add_text("status")
with dpg.window(label = "Preferences", tag = "preferences", modal = True, show = False, width = 400, height = 600):
with dpg.collapsing_header(label = "General"):
with dpg.group(horizontal = True):
dpg.add_text("Font")
dpg.add_text(dpg.get_item_font("Primary"))
dpg.create_viewport(title = "New GUI", width = 1000, height = 800)
dpg.setup_dearpygui()
dpg.show_viewport()
dpg.set_primary_window("Primary", True)
dpg.start_dearpygui()
dpg.destroy_context()

View File

@ -1,37 +1,22 @@
# -*- 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 .file import FoamFile
from .runner import FoamRunner
from .case import FoamCase
from . import presets
from . import commands
from .meshConversion import ideasUnvToFoam
from .meshManipulation import createPatch, transformPoints, checkMesh, renumberMesh
from .miscellaneous import foamDictionary
from .parallelProcessing import decomposePar
from .solvers import potentialFoam, simpleFoam
from .utils import version, foamClean, uniform
__all__ = [
# meshConversion
"ideasUnvToFoam",
# meshManipulation
"createPatch",
"transformPoints",
"checkMesh",
"renumberMesh",
# miscellaneous
"foamDictionary",
# parallelProcessing
"decomposePar",
# solvers
"potentialFoam",
"simpleFoam",
# utils
"version",
"foamClean",
"uniform"
"utils",
"conversion",
"FoamFile",
"FoamRunner",
"FoamCase",
"presets",
"commands"
]

View File

@ -1,51 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import os, sys
import subprocess
import logging
logger = logging.getLogger()
def application(name: str, *args: str, case: str = None, stderr: bool = True, useMPI: bool = False) -> int:
cmd = []
if useMPI:
nprocs = os.cpu_count()
cmd.extend(["mpirun", "-np", str(nprocs), "--oversubscribe"])
cmd.append(name)
if case:
cmd.extend(["-case", case])
if args:
cmd.extend([*args])
logger.info("{}: {}".format(name, [*args]))
logpath = os.path.join(case if case else "", "{}.log".format(name))
with subprocess.Popen(cmd,
stdout = subprocess.PIPE,
stderr = subprocess.PIPE) as p, \
open(logpath, "wb") as logfile:
for line in p.stdout:
#sys.stdout.buffer.write(line)
logfile.write(line)
#for line in p.stderr:
# logfile.write(line)
out, err = p.communicate()
logfile.write(err)
#if err and stderr:
# logger.error("""{}:
# {}""".format(name, str(err, "utf-8")))
return out, err, p.returncode

174
anisotropy/openfoam/case.py Normal file
View File

@ -0,0 +1,174 @@
# -*- coding: utf-8 -*-
from __future__ import annotations
import os
import shutil
import re
import pathlib
import io
from lz.reversal import reverse as lz_reverse
from . import FoamFile
class FoamCase(object):
def __init__(self, files: FoamFile | list[FoamFile] = None, path: str = None):
self.path = path
self._backpath = None
self._files = []
if files is not None:
self.add(files)
def __repr__(self) -> str:
content = [ file.object for file in self._files ]
return "<FoamCase: {}>".format(", ".join(content) or None)
def add(self, files: FoamFile | list[FoamFile]):
if type(files) is not list:
assert type(files) is FoamFile, "passed object is not a FoamFile"
files = [ files ]
for file in files:
assert type(file) is FoamFile, "passed object is not a FoamFile"
assert file.object is not None, "FoamFile object attribute is None"
idn = self.contains(file.object)
if idn is not None:
self._files.pop(idn)
self._files.append(file)
return self
self._files.append(file)
return self
def __add__(self, files: FoamFile | list[FoamFile]):
return self.add(files)
def contains(self, name: str) -> int | None:
for n, file in enumerate(self._files):
if file.object == name:
return n
return None
def write(self, path: str = None):
path = pathlib.Path(path or self.path or "")
for file in self._files:
path_ = path / (
file.location + "/" + file.object
if file.location else file.object
)
file.write(path_.resolve())
def read(self, path: str = None):
path = pathlib.Path(path or self.path or "")
for file in self._files:
path_ = path / (
file.location + "/" + file.object
if file.location else file.object
)
file.read(path_.resolve())
def remove(self, path: str = None):
path = pathlib.Path(path or self.path or "")
for file in self._files:
path_ = path / (
file.location + "/" + file.object
if file.location else file.object
)
file.remove(path_.resolve())
def clean(self, included: list = ["0", "constant", "system"]):
regxs = [
r"^\d+.\d+$",
r"^\d+$",
r"^processor\d+$",
r"^log..+$",
r"^.+.log$",
r"^logs$",
r"^postProcessing$",
r"^polyMesh$"
]
excluded = []
for root, dirs, files in os.walk(os.path.abspath("")):
for _dir in dirs:
excluded += [
os.path.join(root, _dir)
for regx in regxs if re.match(regx, _dir)
]
for file in files:
excluded += [
os.path.join(root, file)
for regx in regxs if re.match(regx, file)
]
for file in excluded:
if os.path.split(file)[1] not in included and os.path.exists(file):
if os.path.isdir(file):
shutil.rmtree(file)
if os.path.isfile(file):
os.remove(file)
def chdir(self, path: str = None):
path = pathlib.Path(path or self.path or "").resolve()
self._backpath = os.getcwd()
os.chdir(path)
def chback(self):
path = pathlib.Path(self._backpath or "").resolve()
os.chdir(path)
def is_converged(self, path: str = None) -> None | bool:
path = pathlib.Path(path or self.path or "").resolve()
controlDict = FoamFile()
if controlDict.exists(path / "system" / "controlDict"):
controlDict.read(path / "system" / "controlDict")
else:
return None
application = controlDict.get("application")
if application is None:
return None
logfile = (
path / f"{ application }.log"
if (path / f"{ application }.log").exists() else (path / f"log.{ application }")
)
status = False
if logfile.exists():
with open(logfile, "r") as infile:
limit = 30
for line in lz_reverse(infile, batch_size = io.DEFAULT_BUFFER_SIZE):
if not line.find("End") < 0:
status = True
if limit <= 0:
status = False
limit -= 1
return status

View File

@ -0,0 +1,211 @@
# -*- coding: utf-8 -*-
from . import FoamRunner
###
# meshConversion
##
def netgenNeutralToFoam(meshfile: str, run: bool = True, **kwargs) -> FoamRunner:
command = "netgenNeutralToFoam"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = [ meshfile ]
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
def ideasUnvToFoam(meshfile: str, run: bool = True, **kwargs) -> FoamRunner:
command = "ideasUnvToFoam"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = [ meshfile ]
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
###
# meshManipulation
##
def createPatch(dictfile: str = None, overwrite: bool = True, run: bool = True, **kwargs) -> FoamRunner:
command = "createPatch"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = []
if dictfile:
args.extend(["-dict", dictfile])
if overwrite:
args.append("-overwrite")
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
def transformPoints(transformations: dict, run: bool = True, **kwargs) -> FoamRunner:
command = "transformPoints"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = []
arg = []
for k, v in transformations.items():
if type(v) == int or type(v) == float:
value = str(v)
elif type(v) == tuple or type(v) == list:
value = "({} {} {})".format(*v)
arg.append("{}={}".format(k, value))
args.append(", ".join(arg))
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
def checkMesh(allGeometry: bool = True, allTopology: bool = True, run: bool = True, **kwargs) -> FoamRunner:
command = "checkMesh"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = []
if allGeometry:
args.append("-allGeometry")
if allTopology:
args.append("-allTopology")
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
def renumberMesh(overwrite: bool = True, run: bool = True, **kwargs) -> FoamRunner:
command = "renumberMesh"
kwargs.update(logpath = kwargs.get("logpath", f"{ command }.log"))
kwargs.update(exit = True)
args = []
if overwrite:
args.append("-overwrite")
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
###
# miscellaneous
##
# def foamDictionary()
###
# parallelProcessing
##
def decomposePar(run: bool = True, **kwargs) -> FoamRunner:
command = "decomposePar"
kwargs.update(logpath = kwargs.get("logpath", f"{command}.log"))
kwargs.update(exit = True)
args = []
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
###
# solvers
##
def potentialFoam(parallel: bool = False, run: bool = True, **kwargs) -> FoamRunner:
command = "potentialFoam"
kwargs.update(logpath = kwargs.get("logpath", f"{command}.log"))
kwargs.update(exit = True)
args = []
if parallel:
args.append("-parallel")
kwargs.update(mpi = True)
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
def simpleFoam(parallel: bool = False, run: bool = True, **kwargs) -> FoamRunner:
command = "simpleFoam"
kwargs.update(logpath = kwargs.get("logpath", f"{command}.log"))
kwargs.update(exit = True)
args = []
if parallel:
args.append("-parallel")
kwargs.update(mpi = True)
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner
###
# postProcessing
##
def postProcess(func: str = None, latestTime: bool = False, run: bool = True, **kwargs) -> FoamRunner:
command = "postProcess"
kwargs.update(logpath=kwargs.get("logpath", f"{command}.log"))
kwargs.update(exit = True)
args = []
if func:
args.extend(["-func", func])
if latestTime:
args.append("-latestTime")
runner = FoamRunner(command, args = args, **kwargs)
if run:
runner.run()
return runner

View File

@ -0,0 +1,120 @@
# -*- coding: utf-8 -*-
import pathlib
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
from PyFoam.Basics.FoamFileGenerator import FoamFileGenerator
import numpy as np
from . import utils
def read_foamfile(filename: str) -> tuple[dict, dict]:
"""Read a FoamFile.
:param filename:
Path to the file.
:return:
Two dictionaries that contains header and content information
of the FoamFile.
"""
path = pathlib.Path(filename).resolve()
ppf = ParsedParameterFile(path)
header = ppf.header or {}
content = ppf.content or {}
if type(content) == dict:
if content.get("internalField"):
content["internalField"] = np.asarray(content["internalField"])
return header, content
def write_foamfile(header: dict, content: dict, filename: str):
"""Write a FoamFile to the file.
:param header:
Header block with the FoamFile metadata.
:param content:
Content block of the FoamFile.
:param filename:
Path to the file.
"""
path = pathlib.Path(filename).resolve()
# preformat
header = (
FoamFileGenerator({}, header = header)
.makeString()[ :-2]
.replace("\n ", "\n" + 4 * " ")
)
content = (
FoamFileGenerator(content)
.makeString()[ :-1]
.replace("\n ", "\n" + 4 * " ")
.replace(" \t// " + 73 * "*" + " //", "")
.replace(" /* empty */ ", "")
)
with open(path, "w") as outfile:
outfile.write(utils.template(header, content) + "\n")
def read_dat(filename: str):
"""Read dat file.
:param filename:
Path to the file.
:return:
Dictionary with arrays. Keys are created according file header
block or numerated with string numbers if header is not found.
"""
path = pathlib.Path(filename).resolve()
header = []
content = []
with open(path, "r") as infile:
for line in infile.readlines():
if line.startswith("#"):
header.append(line)
else:
content.append(line)
columns = []
if header[-1].find(":") < 0:
for column in header[-1].replace("#", "").split("\t"):
columns.append(column.strip())
header.pop(-1)
else:
for column in range(len(content[0].split("\t"))):
columns.append(str(column))
output = {}
for row in header:
key, value = row.replace("#", "").split(":")
try:
value = float(value.strip())
except Exception:
value = value.strip()
output[key.strip()] = value
for column in columns:
output[column] = []
for row in content:
values = row.split("\t")
for column, value in zip(columns, values):
output[column].append(float(value))
for key in output.keys():
output[key] = np.asarray(output[key])
return output

171
anisotropy/openfoam/file.py Normal file
View File

@ -0,0 +1,171 @@
# -*- coding: utf-8 -*-
import pathlib
from . import conversion
class FoamFile(object):
def __init__(
self,
filename: str = None,
_version: float = 2.0,
_format: str = "ascii",
_class: str = "dictionary",
_object: str = None,
_location: str = None
):
"""A FoamFile object.
:param filename:
Can be used as shortcut to set _location and _object,
_location and _object parameters will be ignored.
:param _version:
Version of the file format, current is 2.0.
:param _format:
ASCII or binary representation, currently ascii only
supported.
:param _class:
Class of the file.
:param _object:
Usually contains name of the file.
:param _location:
Path to the parent directory of the file according
to the case root.
"""
if filename:
splitted = filename.split("/")
_object = splitted[-1]
_location = "/".join(splitted[ :-1])
self.header = {
"version": _version,
"format": _format,
"class": _class,
"object": _object
}
self.content = {}
if _location:
self.header["location"] = f'"{ _location }"'
@property
def version(self) -> str:
return self.header.get("version")
@property
def format(self) -> str:
return self.header.get("format")
@property
def class_(self) -> str:
return self.header.get("class")
@property
def object(self) -> str:
return self.header.get("object")
@property
def location(self) -> str:
return (
self.header.get("location").replace('"', "")
if self.header.get("location") else None
)
def __getitem__(self, key):
return self.content[key]
def get(self, key, default = None):
return self.content.get(key, default)
def __setitem__(self, key, value):
self.content[key] = value
def __delitem__(self, key):
del self.content[key]
def update(self, **kwargs):
self.content.update(**kwargs)
def __len__(self):
return len(self.content)
def __iter__(self):
for key in self.content:
yield key
def __repr__(self) -> str:
return "<FoamFile: {}>".format(self.header["object"] or None)
def __add__(self, file):
from . import FoamCase
assert type(file) is FoamFile
return FoamCase([ self, file ])
def exists(self, filename: str = None) -> bool:
filename = filename or (
self.location + "/" + self.object
if self.location else self.object
)
path = pathlib.Path(filename).resolve()
return path.exists() and path.is_file()
def read(self, filename: str = None):
"""Read a FoamFile.
:param filename:
Path to the file. If None, use location and object from header with
the current working directory.
:return:
Self.
"""
filename = filename or (
self.location + "/" + self.object
if self.location else self.object
)
path = pathlib.Path(filename).resolve()
header, content = conversion.read_foamfile(path)
self.header = header
if not header.get("object"):
self.header["object"] = path.parts[-1]
self.content = content
return self
def write(self, filename: str = None):
"""Write a FoamFile to the file.
:param filename:
Path to the file. If None, use location and object from header with
the current working directory.
"""
filename = filename or (
self.location + "/" + self.object
if self.location else self.object
)
path = pathlib.Path(filename).resolve()
path.parent.mkdir(parents = True, exist_ok = True)
conversion.write_foamfile(self.header, self.content, path)
def remove(self, filename: str = None):
"""Remove a FoamFile.
:param filename:
Path to the file. If None, use location and object from header with
the current working directory.
"""
filename = filename or (
self.location + "/" + self.object
if self.location else self.object
)
path = pathlib.Path(filename).resolve()
if path.exists():
pathlib.os.remove(filename)

View File

@ -1,9 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .application import application
def ideasUnvToFoam(mesh: str, case: str = None) -> (str, int):
return application("ideasUnvToFoam", mesh, case = case, stderr = True)

View File

@ -1,42 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .application import application
import re
def createPatch(dictfile: str = None, case: str = None):
args = ["-overwrite"]
if dictfile:
args.extend(["-dict", dictfile])
application("createPatch", *args, case = case, stderr = True)
def transformPoints(scale, case: str = None):
_scale = f"({ scale[0] } { scale[1] } { scale[2] })"
application("transformPoints", "-scale", _scale, case = case, stderr = True)
def checkMesh(case: str = None) -> str:
_, err, returncode = application("checkMesh", "-allGeometry", "-allTopology", case = case, stderr = True)
out = ""
with open("checkMesh.log", "r") as io:
warnings = []
for line in io:
if re.search("\*\*\*", line):
warnings.append(line.replace("***", "").strip())
if warnings:
out = "checkMesh:\n\t{}".format("\n\t".join(warnings))
return out, err, returncode
def renumberMesh(case: str = None):
application("renumberMesh", "-overwrite", useMPI = False, case = case, stderr = True)

View File

@ -1,14 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .application import application
def foamDictionary(filepath: str, entry: str, value: str = None, case: str = None):
args = [filepath, "-entry", entry]
if value:
args.extend(["-set", value])
application("foamDictionary", *args, case = case, stderr = False)

View File

@ -1,9 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .application import application
def decomposePar(case: str = None):
application("decomposePar", case = case, stderr = True)

View File

@ -0,0 +1,223 @@
# -*- coding: utf-8 -*-
from . import FoamFile
def controlDict() -> FoamFile:
ff = FoamFile(filename = "system/controlDict")
ff.content = {
"application": "simpleFoam",
"startFrom": "startTime",
"startTime": 0,
"stopAt": "endTime",
"endTime": 2000,
"deltaT": 1,
"writeControl": "timeStep",
"writeInterval": 100,
"purgeWrite": 0,
"writeFormat": "ascii",
"writePrecision": 6,
"writeCompression": "off",
"timeFormat": "general",
"timePrecision": 6,
"runTimeModifiable": "true"
}
return ff
def fvSolution() -> FoamFile:
ff = FoamFile(filename = "system/fvSolution")
ff.content = {
"solvers": {
"p": {
"solver": "GAMG",
"tolerance": 1e-06,
"relTol": 0.1,
"smoother": "GaussSeidel"
},
"U": {
"solver": "smoothSolver",
"smoother": "symGaussSeidel",
"tolerance": 1e-05,
"relTol": 0.1
}
},
"SIMPLE": {
"nNonOrthogonalCorrectors": 0,
"consistent": "yes",
"residualControl": {
"p": 1e-02,
"U": 1e-03
}
},
"relaxationFactors": {
"fields": {
"p": 0.3
},
"equations": {
"U": 0.7
}
}
}
return ff
def fvSchemes() -> FoamFile:
ff = FoamFile(filename = "system/fvSchemes")
ff.content = {
"ddtSchemes": {
"default": "steadyState"
},
"gradSchemes": {
"default": ("Gauss", "linear")
},
"divSchemes": {
"default": "none",
"div(phi,U)": ("bounded", "Gauss", "linearUpwind", "grad(U)"),
"div((nuEff*dev2(T(grad(U)))))": ("Gauss", "linear"),
"div(nonlinearStress)": ("Gauss", "linear")
},
"laplacianSchemes": {
"default": ("Gauss", "linear", "corrected")
},
"interpolationSchemes": {
"default": "linear"
},
"snGradSchemes": {
"default": "corrected"
}
}
return ff
def transportProperties() -> FoamFile:
ff = FoamFile(filename = "constant/transportProperties")
ff.content = {
"transportModel": "Newtonian",
"nu": 1e-05
}
return ff
def turbulenceProperties() -> FoamFile:
ff = FoamFile(filename = "constant/turbulenceProperties")
ff.content = {
"simulationType": "RAS",
"RAS": {
"RASModel": "kEpsilon",
"turbulence": "on",
"printCoeffs": "on"
}
}
return ff
def p() -> FoamFile:
ff = FoamFile("0/p", _class = "volScalarField")
ff.content = {
"dimensions": "[0 2 -2 0 0 0 0]",
"internalField": "uniform 0",
"boundaryField": {
"inlet": {
"type": "fixedValue",
"value": "uniform 0.001"
},
"outlet": {
"type": "fixedValue",
"value": "uniform 0"
},
"wall": {
"type": "zeroGradient"
}
}
}
return ff
def U() -> FoamFile:
ff = FoamFile("0/U", _class = "volVectorField")
ff.content = {
"dimensions": "[0 1 -1 0 0 0 0]",
"internalField": "uniform (0 0 0)",
"boundaryField": {
"inlet": {
"type": "fixedValue",
"value": "uniform (0 0 -6e-5)"
},
"outlet": {
"type": "zeroGradient",
},
"wall": {
"type": "fixedValue",
"value": "uniform (0 0 0)"
}
}
}
return ff
def createPatchDict() -> FoamFile:
ff = FoamFile("system/createPatchDict")
ff.content = {
"pointSync": False,
"patches": [
{
"name": "inlet",
"patchInfo": {
"type": "patch",
"inGroups": ["inlet"]
},
"constructFrom": "patches",
"patches": ["some_inlet"]
},
{
"name": "output",
"patchInfo": {
"type": "patch",
"inGroups": ["outlet"]
},
"constructFrom": "patches",
"patches": ["some_outlet"]
},
{
"name": "wall",
"patchInfo": {
"type": "wall",
"inGroups": ["wall"]
},
"constructFrom": "patches",
"patches": ["some_wall"]
}
]
}
return ff
def decomposeParDict() -> FoamFile:
ff = FoamFile("system/decomposeParDict")
ff.content = {
"numberOfSubdomains": 4,
"method": "simple",
"coeffs": {
"n": [2, 2, 2]
}
}
return ff

View File

@ -0,0 +1,80 @@
# -*- coding: utf-8 -*-
import os
import subprocess
import logging
logger = logging.getLogger(__name__)
class FoamRunner(object):
def __init__(self, command: str, args: list[str] = None, mpi: bool = False, cwd: str = None, logpath: str = None, exit: bool = False):
self.command = command
self.args = args
self.mpi = mpi
self.cwd = cwd or os.getcwd()
self.logpath = logpath
self.exit = exit
self.output = ""
self.error = ""
self.returncode = 0
def fullcommand(self) -> list[str]:
command = []
if self.mpi:
nprocs = os.cpu_count()
command.extend(["mpirun", "-np", str(nprocs), "--oversubscribe"])
command.append(self.command)
if self.args:
command.extend([ str(arg) for arg in self.args ])
return command
def run(self) -> tuple[str, str, int]:
try:
proc = subprocess.Popen(
self.fullcommand(),
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
encoding = "utf-8",
cwd = self.cwd
)
logger.debug(f"Starting subprocess: { proc.args }")
if self.logpath:
with proc, open(self.logpath, "w") as io:
while True:
output = proc.stdout.read(1)
if output == "" and proc.poll() is not None:
break
if not output == "":
io.write(output)
error = proc.stderr.read()
if not error == "":
self.error = error
io.write(error)
self.returncode = proc.returncode
if self.logpath and self.error:
with open(self.logpath, "a") as io:
io.write(self.error)
except FileNotFoundError as err:
self.error = err.args[1]
self.returncode = 2
logger.error(self.error, exc_info = True)
if not self.returncode == 0 and self.exit:
raise Exception(f"Subprocess failed: { self.error }")
return self.output, self.error, self.returncode

View File

@ -1,34 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from .application import application
import re
def potentialFoam(case: str = None, useMPI: bool = False):
if useMPI:
out, err, returncode = application("potentialFoam", "-parallel", useMPI = True, case = case, stderr = True)
else:
out, err, returncode = application("potentialFoam", case = case, stderr = True)
return out, err, returncode
def simpleFoam(case: str = None, useMPI: bool = False):
if useMPI:
out, err, returncode = application("simpleFoam", "-parallel", useMPI = True, case = case, stderr = True)
else:
out, err, returncode = application("simpleFoam", case = case, stderr = True)
out = ""
with open("simpleFoam.log", "r") as io:
for line in io:
if re.search("solution converged", line):
out = "simpleFoam:\n\t{}".format(line.strip())
return out, err, returncode

View File

@ -1,58 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 -6e-5);
}
outlet
{
type zeroGradient;
}
symetryPlane
{
type fixedValue;
value uniform (0 0 0);
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
strips
{
type fixedValue;
value uniform (0 0 0);
}
defaultFaces
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //

View File

@ -1,55 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type fixedValue;
value uniform 0.001;
}
outlet
{
type fixedValue;
value uniform 0;
}
symetryPlane
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
strips
{
type zeroGradient;
}
defaultFaces
{
type zeroGradient;
}
}
// ************************************************************************* //

View File

@ -1,22 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu 1e-06;
// ************************************************************************* //

View File

@ -1,23 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
// ************************************************************************* //

View File

@ -1,28 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object collapseDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
collapseEdgesCoeffs
{
// Edges shorter than this absolute value will be merged
minimumEdgeLength 2e-7;
// The maximum angle between two edges that share a point attached to
// no other edges
maximumMergeAngle 5;
}
// ************************************************************************* //

View File

@ -1,80 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom latestTime; //startTime;
startTime 0;
stopAt endTime;
endTime 5000;
deltaT 1;
writeControl timeStep;
writeInterval 50;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
flowRatePatch(name=inlet)
{
name inlet;
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl timeStep;
writeInterval 1;
writeFields false;
log false;
regionType patch;
fields ( phi );
operation sum;
}
flowRatePatch(name=outlet)
{
name outlet;
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl timeStep;
writeInterval 1;
writeFields false;
log false;
regionType patch;
fields ( phi );
operation sum;
}
}
// ************************************************************************* //

View File

@ -1,168 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object createPatchDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
pointSync false;
// Patches to create.
patches
(
{
name defaultFaces;
patchInfo
{
type empty;
}
constructFrom patches;
patches (defaultFaces);
}
{
name inlet;
patchInfo
{
type patch;
inGroups (inlet);
}
constructFrom patches;
patches (smesh_inlet);
}
{
name outlet;
patchInfo
{
type patch;
inGroups (outlet);
}
constructFrom patches;
patches (smesh_outlet);
}
{
name wall;
patchInfo
{
type wall;
inGroups (wall);
}
constructFrom patches;
patches (smesh_wall);
}
{
name symetry0;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry0);
}
{
name symetry1;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry1);
}
{
name symetry2;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry2);
}
{
name symetry3;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry3);
}
{
name symetry4;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry4);
}
{
name symetry5;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (smesh_symetry5);
}
{
name strips;
patchInfo
{
type wall;
inGroups (wall);
}
constructFrom patches;
patches (smesh_strips);
}
);
// ************************************************************************* //

View File

@ -1,26 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
numberOfSubdomains 4;
method simple;
coeffs
{
n (2 2 1);
}
// ************************************************************************* //

View File

@ -1,53 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
// grad(Phi) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss linearUpwind grad(U);
div((nuEff*dev2(T(grad(U))))) Gauss linear;
div(nonlinearStress) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
// laplacian(1, Phi) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
// ************************************************************************* //

View File

@ -1,96 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-06;
relTol 0.1;
smoother GaussSeidel;
}
U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-08;
relTol 0.1;
}
/*Phi
{
solver GAMG;
smoother GaussSeidel;
tolerance 1e-08;
relTol 0.01;
}*/
Phi
{
solver GAMG;
smoother DIC;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
tolerance 1e-06;
relTol 0.01;
}
}
potentialFlow
{
nNonOrthogonalCorrectors 20;
PhiRefCell 0;
PhiRefPoint 0;
PhiRefValue 0;
Phi 0;
}
cache
{
grad(U);
}
SIMPLE
{
nNonOrthogonalCorrectors 10;
residualControl
{
p 1e-5;
U 1e-5;
}
}
relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.5;
}
}
// ************************************************************************* //

View File

@ -1,39 +1,61 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from __future__ import annotations
from numpy.typing import ArrayLike
from numpy import ndarray
import os
import shutil
from .application import application
def version() -> str:
return os.environ["WM_PROJECT_VERSION"]
def foamCleanCustom(case: str = None):
rmDirs = ["0", "constant", "system", "postProcessing", "logs"]
rmDirs.extend([ "processor{}".format(n) for n in range(os.cpu_count()) ])
path = case if case else ""
def version() -> str | None:
"""Version of the current OpenFOAM installation.
for d in rmDirs:
if os.path.exists(os.path.join(path, d)):
shutil.rmtree(os.path.join(path, d))
:return:
Version string or None if installation is not found.
"""
return os.environ.get("WM_PROJECT_VERSION")
def foamClean(case: str = None):
rmDirs = ["0", "constant", "system"]
path = case if case else ""
for d in rmDirs:
if os.path.exists(os.path.join(path, d)):
shutil.rmtree(os.path.join(path, d))
def template(header: dict, content: dict) -> str:
"""Render FoamFile with current template.
application("foamCleanTutorials", useMPI = False, case = case, stderr = True)
:param header:
Header block with the FoamFile metadata.
:param content:
Content block of the FoamFile.
:return:
Generated string of the whole FoamFile.
"""
limit = 78
desc = [
"/*--------------------------------*- C++ -*----------------------------------*\\",
"| ========= | |",
"| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |",
"| \\\\ / O peration |",
"| \\\\ / A nd | |",
"| \\\\/ M anipulation | |",
"\\*---------------------------------------------------------------------------*/"
]
desc[3] += " Version: {}".format(version() or "missed")
desc[3] += " " * (limit - len(desc[3])) + "|"
afterheader = "// " + 37 * "* " + "//"
endfile = "// " + 73 * "*" + " //"
def uniform(value) -> str:
if type(value) == list or type(value) == tuple:
return "\n".join([*desc, header, afterheader, content, endfile])
def uniform(value: ArrayLike | float | int) -> str:
"""Convert value to the OpenFOAM uniform representation.
:param value:
Vector or scalar value.
:return:
Uniform string representation.
"""
if type(value) in [list, tuple, ndarray]:
return f"uniform ({ value[0] } { value[1] } { value[2] })"
elif type(value) == int or type(value) == float:
elif type(value) in [int, float]:
return f"uniform { value }"
else:

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.

View File

@ -1,25 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import logging
logger = logging.getLogger("anisotropy")
try:
import GEOM
from salome.geom import geomBuilder
except ImportError:
logger.debug("Trying to get SALOME geometry modules outside SALOME environment. Modules won't be imported.")
if globals().get("geomBuilder"):
geompy = geomBuilder.New()
else:
geompy = None
def getGeom():
return geompy

View File

@ -1,173 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
import logging
logger = logging.getLogger("anisotropy")
try:
import SMESH
from salome.smesh import smeshBuilder
except ImportError:
logger.debug("Trying to get SALOME mesh modules outside SALOME environment. Modules won't be imported.")
if globals().get("smeshBuilder"):
smesh = smeshBuilder.New()
else:
smesh = None
import enum
class Fineness(enum.Enum):
VeryCoarse = 0
Coarse = 1
Moderate = 2
Fine = 3
VeryFine = 4
Custom = 5
def getSmesh():
return smesh
def updateParams(old, new: dict):
old.SetMaxSize(new.get("maxSize", old.GetMaxSize()))
old.SetMinSize(new.get("minSize", old.GetMinSize()))
old.SetFineness(new.get("fineness", old.GetFineness()))
old.SetGrowthRate(new.get("growthRate", old.GetGrowthRate()))
old.SetNbSegPerEdge(new.get("nbSegPerEdge", old.GetNbSegPerEdge()))
old.SetNbSegPerRadius(new.get("nbSegPerRadius", old.GetNbSegPerRadius()))
old.SetChordalErrorEnabled(new.get("chordalErrorEnabled", old.GetChordalErrorEnabled()))
old.SetChordalError(new.get("chordalError", old.GetChordalError()))
old.SetSecondOrder(new.get("secondOrder", old.GetSecondOrder()))
old.SetOptimize(new.get("optimize", old.GetOptimize()))
old.SetQuadAllowed(new.get("quadAllowed", old.GetQuadAllowed()))
old.SetUseSurfaceCurvature(new.get("useSurfaceCurvature", old.GetUseSurfaceCurvature()))
old.SetFuseEdges(new.get("fuseEdges", old.GetFuseEdges()))
old.SetCheckChartBoundary(new.get("checkChartBoundary", old.GetCheckChartBoundary()))
class Mesh(object):
def __init__(self, shape, name = ""):
self.name = name if name else shape.GetName()
self.mesh = smesh.Mesh(shape, self.name)
self.geom = shape
self.algo = None
self.params = None
self.viscousLayers = None
self.submeshes = []
def Tetrahedron(self, **kwargs):
self.algo = self.mesh.Tetrahedron(algo = smeshBuilder.NETGEN_1D2D3D)
self.params = self.algo.Parameters()
self.params = updateParams(self.params, kwargs)
def _extrusionMethod(self, key: str):
return dict(
SURF_OFFSET_SMOOTH = smeshBuilder.SURF_OFFSET_SMOOTH,
FACE_OFFSET = smeshBuilder.FACE_OFFSET,
NODE_OFFSET = smeshBuilder.NODE_OFFSET,
).get(key, smeshBuilder.SURF_OFFSET_SMOOTH)
def ViscousLayers(self,
thickness = 1,
numberOfLayers = 1,
stretchFactor = 0,
faces = [],
isFacesToIgnore = True,
extrMethod = "SURF_OFFSET_SMOOTH",
**kwargs
):
self.viscousLayers = self.algo.ViscousLayers(
thickness,
numberOfLayers,
stretchFactor,
faces,
isFacesToIgnore,
self._extrusionMethod(extrMethod)
)
def Triangle(self, subshape, **kwargs):
submesh = Submesh(self.mesh, subshape)
submesh.algo = self.mesh.Triangle(algo = smeshBuilder.NETGEN_1D2D, geom = subshape)
submesh.mesh = submesh.algo.subm
submesh.params = submesh.algo.Parameters()
submesh.params = updateParams(submesh.params, kwargs)
self.submeshes.append(submesh)
def assignGroups(self, withPrefix = True):
prefix = "smesh_" if withPrefix else ""
for group in self.mesh.geompyD.GetGroups(self.geom):
if group.GetName():
self.mesh.GroupOnGeom(group, f"{ prefix }{ group.GetName() }", SMESH.FACE)
def compute(self):
isDone = self.mesh.Compute()
returncode = int(not isDone)
err = self.mesh.GetComputeErrors()
return "", 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 exportUNV(self, path):
returncode = 0
error = ""
try:
self.mesh.ExportUNV(path)
except Exception as e:
error = e.details.text
returncode = 1
return returncode, error
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()
class Submesh(object):
def __init__(self, father, subshape, name = ""):
self.name = name if name else subshape.GetName()
self.mesh = None
self.geom = subshape
self.algo = None
self.params = None

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

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

View File

@ -1,230 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from math import pi, sqrt
from anisotropy.salomepl import geometry
class BodyCentered(object):
def __init__(self, **kwargs):
self.direction = kwargs.get("direction", [1, 0, 0])
self.theta = kwargs.get("theta", 0.01)
self.L = kwargs.get("L", 1)
self.r0 = kwargs.get("r0", self.L * sqrt(3) / 4)
self.radius = kwargs.get("radius", self.r0 / (1 - self.theta))
self.filletsEnabled = kwargs.get("filletsEnabled", False)
self.fillets = kwargs.get("fillets", 0)
self.volumeCell = None
def build(self):
geompy = geometry.getGeom()
###
# 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 = geompy.MakeVertex(0, 0, 0)
spos1 = (0, 0, 0)
spos2 = (0, 0, 0)
###
# Bounding box
##
if self.direction == [1, 0, 0]:
sk = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, diag)
elif self.direction == [0, 0, 1]:
sk = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, zh)
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVertex(0, 0, 0)
spos1 = (0, 0, 0)
spos2 = (0, 0, 0)
###
# Bounding box
##
sk = geompy.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = geompy.MakeFaceWires([sk.wire()], False)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVectorDXDYDZ(1, 0, 0)
oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), self.radius)
lattice1 = geompy.MakeMultiTranslation2D(grain, ox, self.L, xn, oy, self.L, yn)
lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, self.L, zn)
#grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius)
#lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
#lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn)
lattice2 = geompy.MakeTranslation(lattice1, 0.5 * self.L, 0.5 * self.L, 0.5 * self.L)
grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
grains = geompy.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = geompy.MakeScaleTransform(grains, oo, 1 / scale)
grains = geompy.MakeFilletAll(grains, self.fillets * scale)
###
# Groups
##
shape = geompy.MakeCutList(poreCell, [grains])
shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "bodyCentered")
sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
geompy.UnionIDs(sall,
geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
inletshape = geompy.MakeCutList(inletface, [grains])
inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
geompy.UnionList(inlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
outletshape = geompy.MakeCutList(outletface, [grains])
outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
geompy.UnionList(outlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
symetry = []
for (n, face) in enumerate(symetryface):
name = "symetry" + str(n)
symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
symetryshape = geompy.MakeCutList(face, [grains])
symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
geompy.UnionList(symetry[n], geompy.SubShapeAll(
geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
groups = []
groups.append(inlet)
groups.append(outlet)
groups.extend(symetry)
if self.filletsEnabled:
strips = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "strips")
shapeShell = geompy.ExtractShapes(shape, geompy.ShapeType["SHELL"], True)
stripsShape = geompy.MakeCutList(shapeShell[0], groups + [grainsOrigin])
geompy.UnionList(strips, geompy.SubShapeAll(
geompy.GetInPlace(shape, stripsShape, True), geompy.ShapeType["FACE"]))
groups.append(strips)
wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
groups.append(wall)
return shape, groups

View File

@ -1,230 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from math import pi, sqrt
from anisotropy.salomepl import geometry
class FaceCentered(object):
def __init__(self, **kwargs):
self.direction = kwargs.get("direction", [1, 0, 0])
self.theta = kwargs.get("theta", 0.01)
self.L = kwargs.get("L", 1)
self.r0 = kwargs.get("r0", self.L * sqrt(2) / 4)
self.radius = kwargs.get("radius", self.r0 / (1 - self.theta))
self.filletsEnabled = kwargs.get("filletsEnabled", False)
self.fillets = kwargs.get("fillets", 0)
self.volumeCell = None
def build(self):
geompy = geometry.getGeom()
###
# 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 = geompy.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 = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, length)
elif self.direction == [0, 0, 1]:
sk = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, 2 * zh)
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVertex(0, 0, 0)
spos1 = (-width * (xn - 1), 0, -width * (zn - 2))
spos2 = (-width * xn, 0, -width * (zn - 1))
###
# Bounding box
##
sk = geompy.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = geompy.MakeFaceWires([sk.wire()], False)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVectorDXDYDZ(1, 0, 0)
oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
xy = geompy.MakeVectorDXDYDZ(1, 1, 0)
xmy = geompy.MakeVectorDXDYDZ(1, -1, 0)
grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), self.radius)
lattice1 = geompy.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn)
lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, self.L, zn - 1)
grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), self.radius)
lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1)
lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, self.L, zn)
grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
grains = geompy.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = geompy.MakeScaleTransform(grains, oo, 1 / scale)
grains = geompy.MakeFilletAll(grains, self.fillets * scale)
###
# Groups
##
shape = geompy.MakeCutList(poreCell, [grains])
shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "faceCentered")
sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
geompy.UnionIDs(sall,
geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
inletshape = geompy.MakeCutList(inletface, [grains])
inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
geompy.UnionList(inlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
outletshape = geompy.MakeCutList(outletface, [grains])
outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
geompy.UnionList(outlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
symetry = []
for (n, face) in enumerate(symetryface):
name = "symetry" + str(n)
symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
symetryshape = geompy.MakeCutList(face, [grains])
symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
geompy.UnionList(symetry[n], geompy.SubShapeAll(
geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
groups = []
groups.append(inlet)
groups.append(outlet)
groups.extend(symetry)
if self.filletsEnabled:
strips = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "strips")
shapeShell = geompy.ExtractShapes(shape, geompy.ShapeType["SHELL"], True)
stripsShape = geompy.MakeCutList(shapeShell[0], groups + [grainsOrigin])
geompy.UnionList(strips, geompy.SubShapeAll(
geompy.GetInPlace(shape, stripsShape, True), geompy.ShapeType["FACE"]))
groups.append(strips)
wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
groups.append(wall)
return shape, groups

View File

@ -1,215 +0,0 @@
# -*- coding: utf-8 -*-
# This file is part of anisotropy.
# License: GNU GPL version 3, see the file "LICENSE" for details.
from math import pi, sqrt
from anisotropy.salomepl import geometry
class Simple(object):
def __init__(self, **kwargs):
self.direction = kwargs.get("direction", [1, 0, 0])
self.theta = kwargs.get("theta", 0.01)
self.r0 = kwargs.get("r0", 1)
self.L = kwargs.get("L", 2 * self.r0)
self.radius = kwargs.get("radius", self.r0 / (1 - self.theta))
self.filletsEnabled = kwargs.get("filletsEnabled", False)
self.fillets = kwargs.get("fillets", 0)
self.volumeCell = None
def build(self):
geompy = geometry.getGeom()
###
# 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 = geompy.MakeVertex(0, 0, 0)
###
# Bounding box
##
if self.direction == [1, 0, 0]:
sk = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, width)
elif self.direction == [0, 0, 1]:
sk = geompy.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 = geompy.MakeFaceWires([sk.wire()], True)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, height)
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVertex(0, 0, 0)
###
# Bounding box
##
sk = geompy.Sketcher3D()
for p in point:
sk.addPointsAbsolute(*p)
inletface = geompy.MakeFaceWires([sk.wire()], False)
vecflow = geompy.GetNormal(inletface)
poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
[_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
inletface = geompy.MakeScaleTransform(inletface, oo, scale)
poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
faces = geompy.ExtractShapes(poreCell, geompy.ShapeType["FACE"], False)
symetryface = []
for face in faces:
norm = geompy.GetNormal(face)
angle = round(geompy.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 = geompy.MakeVectorDXDYDZ(1, 0, 0)
oy = geompy.MakeVectorDXDYDZ(0, 1, 0)
oz = geompy.MakeVectorDXDYDZ(0, 0, 1)
grain = geompy.MakeSphereR(self.radius)
lattice = geompy.MakeMultiTranslation2D(grain, ox, self.L, xn, oy, self.L, yn)
lattice = geompy.MakeMultiTranslation1D(lattice, oz, self.L, zn)
grains = geompy.ExtractShapes(lattice, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
grains = geompy.MakeScaleTransform(grains, oo, scale)
grainsOrigin = None
if self.filletsEnabled:
grainsOrigin = geompy.MakeScaleTransform(grains, oo, 1 / scale)
grains = geompy.MakeFilletAll(grains, self.fillets * scale)
###
# Groups
##
shape = geompy.MakeCutList(poreCell, [grains])
shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "simple")
sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"])
geompy.UnionIDs(sall,
geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"]))
inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet")
inletshape = geompy.MakeCutList(inletface, [grains])
inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale)
geompy.UnionList(inlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"]))
outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet")
outletshape = geompy.MakeCutList(outletface, [grains])
outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale)
geompy.UnionList(outlet, geompy.SubShapeAll(
geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"]))
symetry = []
for (n, face) in enumerate(symetryface):
name = "symetry" + str(n)
symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name))
symetryshape = geompy.MakeCutList(face, [grains])
symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale)
geompy.UnionList(symetry[n], geompy.SubShapeAll(
geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"]))
groups = []
groups.append(inlet)
groups.append(outlet)
groups.extend(symetry)
if self.filletsEnabled:
strips = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "strips")
shapeShell = geompy.ExtractShapes(shape, geompy.ShapeType["SHELL"], True)
stripsShape = geompy.MakeCutList(shapeShell[0], groups + [grainsOrigin])
geompy.UnionList(strips, geompy.SubShapeAll(
geompy.GetInPlace(shape, stripsShape, True), geompy.ShapeType["FACE"]))
groups.append(strips)
wall = geompy.CutListOfGroups([sall], groups, theName = "wall")
groups.append(wall)
return shape, groups

View File

@ -0,0 +1,21 @@
# -*- coding: utf-8 -*-
"""
Shaping is a library for using OCC shapes, provides more convient
functionality with power NumPy and Python OOP and contains interesing
primitives.
"""
from . import utils
from .shape import Shape
from .periodic import Periodic
from . import primitives
__all__ = [
"utils",
"primitives",
"Shape",
"Periodic"
]

View File

@ -0,0 +1,69 @@
# -*- coding: utf-8 -*-
import numpy as np
from . import Shape
class Periodic(Shape):
def __init__(
self,
alpha: float = None,
r0: float = 1,
L: float = None,
filletsEnabled: bool = True,
gamma = None,
**kwargs
):
"""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:
If True, calculate fillets beetween spheres.
"""
Shape.__init__(self)
# Parameters
self.alpha = alpha
self.r0 = r0
self.L = L
# for lattice
# self.theta = 0.5 * pi
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
self.maximize = 1e+2
self.minimize = 1e-2
@property
def radius(self) -> float:
"""Spheres radius
:return:
Radius.
"""
return self.r0 / (1 - self.alpha)
@property
def fillets(self):
"""Fillets radius beetween spheres.
:return:
Radius.
"""
analytical = self.r0 * np.sqrt(2) / np.sqrt(1 - np.cos(self.gamma)) - self.radius
rTol = 3
minRound = np.fix(10 ** rTol * analytical) * 10 ** -rTol
return minRound * self.filletsScale

View File

@ -0,0 +1,510 @@
# -*- 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)
pc.length = 0.
# 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))
pc.lattice = np.sum(spheres)
if pc.filletsEnabled:
pc.lattice = pc.lattice.Scale(zero, pc.maximize)
pc.lattice = (
pc.lattice
.MakeFillet(pc.lattice.edges, pc.fillets * pc.maximize)
.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))
pc.length = extr
# 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
pc.length = 0.
# 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))
pc.lattice = np.sum(spheres)
if pc.filletsEnabled:
pc.lattice = pc.lattice.Scale(zero, pc.maximize)
pc.lattice = (
pc.lattice
.MakeFillet(pc.lattice.edges, pc.fillets * pc.maximize)
.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))
pc.length = extr
# 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 = pc.r0 * 4 / np.sqrt(2)
pc.direction = np.asarray(direction)
pc.length = 0.
# 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
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)
)
pc.lattice = (
np.sum(spheres)
.Move(ng_occ.Vec(-pc.r0 * 2, -pc.r0 * 2, 0))
.Rotate(ng_occ.Axis(zero, ng_occ.Z), 45)
)
if pc.filletsEnabled:
pc.lattice = pc.lattice.Scale(zero, pc.maximize)
pc.lattice = (
pc.lattice
.MakeFillet(pc.lattice.edges, pc.fillets * pc.maximize)
.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))
pc.length = extr
# 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

121
anisotropy/shaping/shape.py Normal file
View File

@ -0,0 +1,121 @@
# -*- coding: utf-8 -*-
from __future__ import annotations
from numpy import ndarray
from os import PathLike
import numpy as np
import netgen.occ as ng_occ
import pathlib
from . import utils
class Shape(object):
def __init__(self):
"""A Shape object contains OCC shape.
"""
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 to the file.
Supported formats: step.
:param filename:
Path of the file.
"""
path = pathlib.Path(filename).resolve()
ext = path.suffix[1: ]
if ext == "step":
self.shape.WriteStep(str(path))
else:
raise NotImplementedError(f"Shape format '{ ext }' is not supported")
def read(self, filename: PathLike):
"""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 = 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
) -> 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.
"""
if group:
patches_ = {}
for idn, face in enumerate(self.shape.faces):
if shiftIndex:
idn += 1
item = idn if not prefix else prefix + str(idn)
if patches_.get(face.name):
patches_[face.name].append(item)
else:
patches_[face.name] = [ item ]
else:
patches_ = []
for idn, face in enumerate(self.shape.faces):
if shiftIndex:
idn += 1
item = idn if not prefix else prefix + str(idn)
patches_.append((item, face.name))
return patches_

View File

@ -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))

View File

@ -0,0 +1,8 @@
# -*- coding: utf-8 -*-
from .onephase import FlowOnePhase
__all__ = [
"FlowOnePhase"
]

View File

@ -0,0 +1,268 @@
# -*- coding: utf-8 -*-
import numpy as np
import anisotropy.openfoam as of
from anisotropy.openfoam import presets
from anisotropy.openfoam import commands
class FlowOnePhase:
def __init__(
self,
path: str = None,
direction: list[float] = None,
patches: dict = None,
pressureInlet: float = None,
pressureOutlet: float = None,
pressureInternal: float = None,
velocityInlet: float = None,
velocityOutlet: float = None,
velocityInternal: float = None,
density: float = None,
viscosityKinematic: float = None,
scale: list[float] = None,
**kwargs
):
self.direction = direction
self.patches = patches
self.path = path
self.case = None
self.boundaryNames = ["inlet", "outlet", "symetry", "wall"]
self.pressureInlet = pressureInlet
self.pressureOutlet = pressureOutlet
self.pressureInternal = pressureInternal
self.velocityInlet = velocityInlet
self.velocityOutlet = velocityOutlet
self.velocityInternal = velocityInternal
self.density = density
self.viscosityKinematic = viscosityKinematic
self.scale = scale
def controlDict(self) -> of.FoamFile:
ff = presets.controlDict()
ff.update(
startFrom = "latestTime",
endTime = 5000,
writeInterval = 100,
runTimeModifiable = "true"
)
return ff
def fvSchemes(self) -> of.FoamFile:
ff = presets.fvSchemes()
return ff
def fvSolution(self) -> of.FoamFile:
ff = presets.fvSolution()
ff["solvers"]["U"].update(
nSweeps = 2,
tolerance = 1e-08,
smoother = "GaussSeidel"
)
ff["solvers"]["Phi"] = dict(
solver = "GAMG",
smoother = "DIC",
cacheAgglomeration = "yes",
agglomerator = "faceAreaPair",
nCellsInCoarsestLevel = 10,
mergeLevels = 1,
tolerance = 1e-06,
relTol = 0.01
)
ff["potentialFlow"] = dict(
nNonOrthogonalCorrectors = 13,
PhiRefCell = 0,
PhiRefPoint = 0,
PhiRefValue = 0,
Phi = 0
)
ff["cache"] = { "grad(U)": None }
ff["SIMPLE"].update(
nNonOrthogonalCorrectors = 6,
residualControl = dict(
p = 1e-05,
U = 1e-05
)
)
ff["relaxationFactors"]["equations"]["U"] = 0.5
return ff
def transportProperties(self) -> of.FoamFile:
ff = presets.transportProperties()
ff.update(
nu = self.viscosityKinematic
)
return ff
def turbulenceProperties(self) -> of.FoamFile:
ff = presets.turbulenceProperties()
ff.content = dict(
simulationType = "laminar"
)
return ff
def p(self) -> of.FoamFile:
ff = presets.p()
ff["boundaryField"] = {}
for boundary in self.boundaryNames:
if boundary == "inlet":
ff["boundaryField"][boundary] = dict(
type = "fixedValue",
value = of.utils.uniform(self.pressureInlet / self.density)
)
elif boundary == "outlet":
ff["boundaryField"][boundary] = dict(
type = "fixedValue",
value = of.utils.uniform(self.pressureOutlet / self.density)
)
else:
ff["boundaryField"][boundary] = dict(
type = "zeroGradient"
)
return ff
def U_approx(self) -> of.FoamFile:
ff = presets.U()
ff["boundaryField"] = {}
for boundary in self.boundaryNames:
if boundary == "inlet":
ff["boundaryField"][boundary] = dict(
type = "fixedValue",
value = of.utils.uniform(np.array(self.direction) * -6e-5)
)
elif boundary == "outlet":
ff["boundaryField"][boundary] = dict(
type = "zeroGradient",
)
else:
ff["boundaryField"][boundary] = dict(
type = "fixedValue",
value = of.utils.uniform([0., 0., 0.])
)
return ff
def U(self) -> of.FoamFile:
ff = presets.U()
ff["boundaryField"] = {}
for boundary in self.boundaryNames:
if boundary == "inlet":
ff["boundaryField"][boundary] = dict(
type = "pressureInletVelocity",
value = of.utils.uniform(self.velocityInlet)
)
elif boundary == "outlet":
ff["boundaryField"][boundary] = dict(
type = "zeroGradient",
)
else:
ff["boundaryField"][boundary] = dict(
type = "fixedValue",
value = of.utils.uniform([0., 0., 0.])
)
return ff
def createPatchDict(self) -> of.FoamFile:
# initial 43 unnamed patches ->
# 6 named patches (inlet, outlet, wall, symetry 0 to 3 or 5) ->
# 4 inGroups (inlet, outlet, wall, symetry)
ff = presets.createPatchDict()
ff["patches"] = []
for name in self.patches.keys():
if name == "inlet":
patchGroup = "inlet"
patchType = "patch"
elif name == "outlet":
patchGroup = "outlet"
patchType = "patch"
elif name == "wall":
patchGroup = "wall"
patchType = "wall"
else:
patchGroup = "symetry"
patchType = "symetryPlane"
ff["patches"].append({
"name": name,
"patchInfo": {
"type": patchType,
"inGroups": [patchGroup]
},
"constructFrom": "patches",
"patches": self.patches[name]
})
return ff
def generate(self, approximation: bool = False, meshfile: str = "mesh.mesh"):
self.case = of.FoamCase(
path = self.path,
files = [
self.controlDict(),
self.fvSchemes(),
self.fvSolution(),
self.transportProperties(),
self.turbulenceProperties(),
self.p()
]
)
if self.patches is not None:
self.case += self.createPatchDict()
self.case += self.U_approx() if approximation else self.U()
self.case.write(self.path)
self.case.chdir()
commands.netgenNeutralToFoam(meshfile)
# TODO: contain
if self.case.contains("createPatchDict"):
commands.createPatch()
commands.checkMesh()
if self.scale is not None:
commands.transformPoints({
"scale": self.scale
})
commands.renumberMesh()
if approximation:
commands.potentialFoam()
# replace velocity for the main simulation
self.case += self.U()
self.case.write(self.path)
commands.simpleFoam()
self.case.chback()

3
data/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
*.png
*.tiff
*.sql

916
data/analyze.ipynb Normal file

File diff suppressed because one or more lines are too long

BIN
data/anisotropy.db Normal file

Binary file not shown.

57
data/anisotropy.toml Normal file
View File

@ -0,0 +1,57 @@
[[structures]]
label = "simple"
alpha = [ 0.01, 0.27,]
alphaStep = 0.005
directions = [ [ 1.0, 0.0, 0.0,], [ 0.0, 0.0, 1.0,], [ 1.0, 1.0, 1.0,],]
r0 = 1
filletsEnabled = true
pressureInlet = 1
pressureOutlet = 0
pressureInternal = 0
velocityInlet = [ 0.0, 0.0, 0.0,]
velocityInternal = [ 0.0, 0.0, 0.0,]
density = 1000
viscosity = 0.001
scale = [ 1e-5, 1e-5, 1e-5,]
[[structures]]
label = "bodyCentered"
alpha = [ 0.01, 0.16,]
alphaStep = 0.005
directions = [ [ 1.0, 0.0, 0.0,], [ 0.0, 0.0, 1.0,], [ 1.0, 1.0, 1.0,],]
r0 = 1
filletsEnabled = true
pressureInlet = 1
pressureOutlet = 0
pressureInternal = 0
velocityInlet = [ 0.0, 0.0, 0.0,]
velocityInternal = [ 0.0, 0.0, 0.0,]
density = 1000
viscosity = 0.001
scale = [ 1e-5, 1e-5, 1e-5,]
[[structures]]
label = "faceCentered"
alpha = [ 0.01, 0.12,]
alphaStep = 0.005
directions = [ [ 1.0, 0.0, 0.0,], [ 0.0, 0.0, 1.0,], [ 1.0, 1.0, 1.0,],]
r0 = 1
filletsEnabled = true
pressureInlet = 1
pressureOutlet = 0
pressureInternal = 0
velocityInlet = [ 0.0, 0.0, 0.0,]
velocityInternal = [ 0.0, 0.0, 0.0,]
density = 1000
viscosity = 0.001
scale = [ 1e-5, 1e-5, 1e-5,]
[options]
nprocs = 4
stage = "shape"
overwrite = false
database = "anisotropy.db"
build = "build"
logs = "logs"
shapefile = "shape.step"
meshfile = "mesh.mesh"

View File

@ -20,6 +20,7 @@ help:
# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
%: Makefile
@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
[ -d "$(BUILDDIR)/html" ] && cp -vf "$(BUILDDIR)/html" anisotropy/
update:
@$(SPHINXAPIDOC) -f -o "$(SOURCEDIR)" "$(PROJECTDIR)"

Some files were not shown because too many files have changed in this diff Show More