Merge branch 'restruction'

Succesful restruction
This commit is contained in:
L-Nafaryus 2021-06-01 12:01:23 +05:00
commit f0d9d8b4fc
70 changed files with 780 additions and 10087 deletions

1
.gitignore vendored
View File

@ -2,6 +2,7 @@ __pycache__
build/ build/
logs/ logs/
storage/ storage/
env/
*.gz *.gz
*.xz *.xz
*.fls *.fls

56
TODO.md
View File

@ -1,56 +0,0 @@
## Usefull utils
- createPatch
- polyDualMesh
- fvSolution: cellLimited
- collapseDict: collapseEdges
- renumberMesh
- processorField
## Errors
- salome:
th. 139990926538304 -
Trace /volatile/salome/jenkins/workspace/Salome_master_CO7/SALOME-9.6.0-CO7/SOURCES/SMESH/src/SMESH/SMESH_subMesh.cxx [2005] :
NETGEN_2D3D failed on sub-shape #1 with error COMPERR_BAD_INPUT_MESH
"NgException at Volume meshing: Stop meshing since surface mesh not consistent Some edges multiple times in surface mesh"
th. 140588498282048 -
Trace /volatile/salome/jenkins/workspace/Salome_master_CO7/SALOME-9.6.0-CO7/SOURCES/SMESH/src/SMESH/SMESH_subMesh.cxx [2005] :
NETGEN_2D3D failed on sub-shape #47 with error COMPERR_WARNING
"Thickness 0.001 of viscous layers not reached, average reached thickness is 0.000928207"
th. 139986338838080 -
Trace /volatile/salome/jenkins/workspace/Salome_master_CO7/SALOME-9.6.0-CO7/SOURCES/SMESH/src/SMESH/SMESH_subMesh.cxx [2005] :
NETGEN_2D3D failed on sub-shape #1 with error COMPERR_BAD_INPUT_MESH
"NgException at Volume meshing: Stop meshing since boundary mesh is overlapping Intersecting triangles"
## 1.03.21
- [x] boundary type (wall or symetryPlane)
- [x] restruct for ways
- [x] build alpha = 0.01 .. 0.13
- [x] ! symetryPlane -> cyclicAMI
## 3.03.21
- [x] configure salome server, ports, etc.
- [ ] less processes for salome, optimization.
## 4.03.21
- [x] 3rd direction
- [x] createPatch(Dict)
- [ ] views (mesh, ..)
- [x] alpha for simpleCubic [0.01 .. 0.28]
- [x] translation vector (cyclicAMI)
- [ ] BUG: angle between the direction vector and the normal to inlet is ~1.4e-14
- [x] Another solution
- [ ] BUG: ideasUnvToFoam not working with param '-case PATH'
- [x] Temporary sulution via os.chdir(PATH)
## 6.03.21
- [ ] ERROR: MakeFuseList with alpha > 0.2
## 7.03.21
- [x] Split the symetryPlane to 4 faces
## 11.03.21
- [x] Dual test for cyclicAMI

View File

@ -1,38 +1,48 @@
import os, sys import os, sys
from collections import namedtuple
import time import time
import logging
from datetime import timedelta from datetime import timedelta
import multiprocessing
import shutil import shutil
import config ROOT = "/".join(__file__.split("/")[:-2])
from config import logger sys.path.append(os.path.abspath(ROOT))
#from src import applogger
from src import utils
from src import salome_utils
from src import foam_utils
from utils import struct
import toml
import logging
CONFIG = os.path.join(ROOT, "conf/config.toml")
config = struct(toml.load(CONFIG))
#CONFIG = os.path.abspath("../conf/config.toml")
#config = struct(toml.load(CONFIG))
LOG = os.path.join(ROOT, "logs")
if not os.path.exists(LOG):
os.makedirs(LOG)
BUILD = os.path.join(ROOT, "build")
if not os.path.exists(BUILD):
os.makedirs(BUILD)
logging.basicConfig(
level = logging.INFO,
format = config.logger.format,
handlers = [
logging.StreamHandler(),
logging.FileHandler(f"{ LOG }/{ config.logger.name }.log")
]
)
logger = logging.getLogger(config.logger.name)
def main(): def main():
if not os.path.exists(config.LOG): if checkEnv():
os.makedirs(config.LOG)
if not os.path.exists(config.BUILD):
os.makedirs(config.BUILD)
#global logger
#logger = applogger.Logger()
check = checkEnv()
if check:
return return
tasks = createTasks() tasks = createTasks()
for task in tasks: for task in tasks:
logger.fancyline() logger.info("-" * 80)
logger.info(f"""main: logger.info(f"""main:
task:\t{tasks.index(task) + 1} / {len(tasks)} task:\t{tasks.index(task) + 1} / {len(tasks)}
cpu count:\t{os.cpu_count()} cpu count:\t{os.cpu_count()}
@ -55,7 +65,7 @@ def main():
if not returncode: if not returncode:
task.flow = True task.flow = True
with open(os.path.join(config.LOG, "tasks.log"), "a") as io: with open(os.path.join(LOG, "tasks.log"), "a") as io:
idx = tasks.index(task) idx = tasks.index(task)
io.write(f"""Task {idx}: io.write(f"""Task {idx}:
structure:\t{task.structure} structure:\t{task.structure}
@ -65,67 +75,34 @@ def main():
flow:\t{task.flow}\n""") flow:\t{task.flow}\n""")
logger.info(f"Warnings: {logger.warnings}\tErrors: {logger.errors}") #logger.info(f"Warnings: {logger.warnings}\tErrors: {logger.errors}")
def checkEnv():
missed = False
try:
pythonVersion = "Python {}".format(sys.version.split(" ")[0])
salomeVersion = salome_utils.salomeVersion()
foamVersion = foam_utils.foamVersion()
except Exception:
logger.critical("Missed environment")
missed = True
else:
logger.info(f"environment:\n\t{pythonVersion}\n\t{salomeVersion}\n\t{foamVersion}")
finally:
return missed
class Task:
def __init__(self, **kwargs):
for (k, v) in kwargs.items():
setattr(self, k, v)
def createTasks(): def createTasks():
#Task = namedtuple("Task", ["structure", "theta", "fillet", "direction", "export"])
tasks = [] tasks = []
structures = [ getattr(config, s)() for s in config.structures ]
for structure in structures: for structure in config.base.__dict__.keys():
for direction in structure.directions: if getattr(config.base, structure):
for theta in structure.theta: for direction in getattr(config, structure).geometry.directions:
name = type(structure).__name__ for theta in getattr(config, structure).theta:
export = os.path.join( task = struct(
config.BUILD, structure = structure,
name, theta = theta,
"direction-{}{}{}".format(*direction), fillet = getattr(config, structure).geometry.fillet,
"theta-{}".format(theta) direction = direction,
) export = os.path.join(ROOT, f"{ BUILD }/{ structure }/direction-{ direction[0] }{ direction [1] }{ direction [2] }/theta-{ theta }"),
mesh = False,
flow = False
)
task = Task( tasks.append(task)
structure = name,
theta = theta,
fillet = structure.fillet,
direction = direction,
export = export,
mesh = False,
flow = False
)
tasks.append(task)
return tasks return tasks
from salomepl.utils import runExecute, salomeVersion
def createMesh(task): def createMesh(task):
scriptpath = os.path.join(config.ROOT, "samples/__init__.py") scriptpath = os.path.join(ROOT, "salomepl/genmesh.py")
port = 2810 port = 2810
stime = time.monotonic() stime = time.monotonic()
@ -135,23 +112,25 @@ def createMesh(task):
int(task.fillet), int(task.fillet),
"".join([str(coord) for coord in task.direction]), "".join([str(coord) for coord in task.direction]),
os.path.join(task.export, "mesh.unv"), os.path.join(task.export, "mesh.unv"),
config.ROOT ROOT
) )
returncode = salome_utils.runExecute(port, scriptpath, *args) returncode = runExecute(port, scriptpath, *args)
etime = time.monotonic() etime = time.monotonic()
logger.info("createMesh: elapsed time: {}".format(timedelta(seconds = etime - stime))) logger.info("createMesh: elapsed time: {}".format(timedelta(seconds = etime - stime)))
import openfoam
def calculate(task): def calculate(task):
foamCase = [ "0", "constant", "system" ] foamCase = [ "0", "constant", "system" ]
os.chdir(task.export) os.chdir(task.export)
foam_utils.foamClean() openfoam.foamClean()
for d in foamCase: for d in foamCase:
shutil.copytree( shutil.copytree(
os.path.join(config.ROOT, "src/cubicFoam", d), os.path.join(ROOT, "openfoam/template", d),
os.path.join(task.export, d) os.path.join(task.export, d)
) )
@ -161,42 +140,62 @@ def calculate(task):
logger.critical(f"calculate: missed 'mesh.unv'") logger.critical(f"calculate: missed 'mesh.unv'")
return return
_, returncode = foam_utils.ideasUnvToFoam("mesh.unv") _, returncode = openfoam.ideasUnvToFoam("mesh.unv")
if returncode: if returncode:
os.chdir(config.ROOT) os.chdir(config.ROOT)
return returncode return returncode
foam_utils.createPatch(dictfile = "system/createPatchDict.symetry") openfoam.createPatch(dictfile = "system/createPatchDict.symetry")
foam_utils.foamDictionary("constant/polyMesh/boundary", "entry0.defaultFaces.type", "wall") openfoam.foamDictionary("constant/polyMesh/boundary", "entry0.defaultFaces.type", "wall")
foam_utils.foamDictionary("constant/polyMesh/boundary", "entry0.defaultFaces.inGroups", "1 (wall)") openfoam.foamDictionary("constant/polyMesh/boundary", "entry0.defaultFaces.inGroups", "1 (wall)")
foam_utils.checkMesh() openfoam.checkMesh()
scale = (1e-5, 1e-5, 1e-5) scale = (1e-5, 1e-5, 1e-5)
foam_utils.transformPoints(scale) openfoam.transformPoints(scale)
foam_utils.decomposePar() openfoam.decomposePar()
foam_utils.renumberMesh() openfoam.renumberMesh()
foam_utils.potentialFoam() openfoam.potentialFoam()
for n in range(os.cpu_count()): for n in range(os.cpu_count()):
foam_utils.foamDictionary(f"processor{n}/0/U", "boundaryField.inlet.type", "pressureInletVelocity") openfoam.foamDictionary(f"processor{n}/0/U", "boundaryField.inlet.type", "pressureInletVelocity")
foam_utils.foamDictionary(f"processor{n}/0/U", "boundaryField.inlet.value", "uniform (0 0 0)") openfoam.foamDictionary(f"processor{n}/0/U", "boundaryField.inlet.value", "uniform (0 0 0)")
returncode = foam_utils.simpleFoam() returncode, out = openfoam.simpleFoam()
if out:
logger.info(out)
os.chdir(config.ROOT) os.chdir(ROOT)
etime = time.monotonic() etime = time.monotonic()
logger.info("calculate: elapsed time: {}".format(timedelta(seconds = etime - stime))) logger.info("calculate: elapsed time: {}".format(timedelta(seconds = etime - stime)))
return returncode return returncode
def checkEnv():
missed = False
try:
pythonVersion = "Python {}".format(sys.version.split(" ")[0])
salomeplVersion = salomeVersion()
openfoamVersion = openfoam.foamVersion()
except Exception as e:
logger.critical("Missed environment %s", e)
missed = True
else:
logger.info(f"environment:\n\t{pythonVersion}\n\t{salomeplVersion}\n\t{openfoamVersion}")
finally:
return missed
def postprocessing(tasks): def postprocessing(tasks):
@ -223,4 +222,3 @@ def postprocessing(tasks):
if __name__ == "__main__": if __name__ == "__main__":
main() main()

153
anisotropy/utils.py Normal file
View File

@ -0,0 +1,153 @@
import logging
from multiprocessing import Queue, Process, cpu_count
import socket
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:
for (k, v) in kwargs.items():
setattr(self, k, v)
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)
class Logger:
def __init__(self, name, logpath):
logging.basicConfig(
level = logging.INFO,
format = "%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler(logpath)
]
)
self.logger = logging.getLogger(name)
self.warnings = 0
self.errors = 0
self.criticals = 0
self.exceptions = 0
def info(self, *args):
self.logger.info(*args)
def warning(self, *args):
self.warnings += 1
self.logger.warning(*args)
def error(self, *args):
self.errors += 1
self.logger.error(*args)
def critical(self, *args):
self.criticals += 1
self.logger.critical(*args)
def exception(self, *args):
self.exceptions += 1
self.logger.exception(*args)
def fancyline(self):
self.logger.info("-" * 80)
def queue(cmd, qin, qout, *args):
while True:
# Get item from the queue
pos, var = qin.get()
# 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()
logging.info("cpu count: {}".format(np))
logging.info("var: {}".format(var))
logging.info("cmd: {}".format(cmd))
# 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

50
bin/anisotropy Executable file
View File

@ -0,0 +1,50 @@
#!/usr/bin/env bash
DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )/.."
cd ${DIR}
anisotropy-help()
{
echo "Usage:"
echo " anisotropy <command> [options]"
echo ""
echo "Commands:"
echo " clean Clean project directory via git (.gitignore)."
echo " init Activate environment and install dependencies."
echo " run Start calculations and etc."
echo " help Show help for commands."
echo ""
echo "Options:"
echo " -c, --config <path-to-file> Use custom configuration file."
}
case $1 in
clean)
git clean -fdx
;;
init)
python -m venv env
source env/bin/activate
python -m pip install --upgrade pip
python -m pip install -r requirements.txt
deactivate
;;
run)
source ${OPENFOAM}
source env/bin/activate
python ${DIR}/anisotropy/anisotropy.py
deactivate
;;
help)
anisotropy-help
exit 1
;;
*)
echo "Unknown command."
anisotropy-help
exit 1
;;
esac

105
conf/config.toml Normal file
View File

@ -0,0 +1,105 @@
[logger]
name = "anisotropy"
format = "%(levelname)s: %(message)s"
[base]
simple = true
bodyCentered = true
faceCentered = true
[simple]
theta = [0.01, 0.28]
[simple.geometry]
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
fillet = true
[simple.mesh]
fineness = 3
minSize = 0.01
maxSize = 0.1
growthRate = 0.5
nbSegPerEdge = 0.5
nbSegPerRadius = 0.5
chordalErrorEnabled = false
chordalError = -1
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
thickness = 0.005
numberOfLayers = 2
stretchFactor = 1.2
isFacesToIgnore = true
[bodyCentered]
theta = [0.01, 0.18]
[bodyCentered.geometry]
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
fillet = true
[bodyCentered.mesh]
fineness = 3
minSize = 0.005
maxSize = 0.05
growthRate = 0.5
nbSegPerEdge = 0.5
nbSegPerRadius = 0.5
chordalErrorEnabled = false
chordalError = -1
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
thickness = 0.005
numberOfLayers = 2
stretchFactor = 1.2
isFacesToIgnore = true
[faceCentered]
theta = [0.01, 0.13]
[faceCentered.geometry]
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
fillet = true
[faceCentered.mesh]
fineness = 3
minSize = 0.005
maxSize = 0.05
growthRate = 0.5
nbSegPerEdge = 0.5
nbSegPerRadius = 0.5
chordalErrorEnabled = false
chordalError = -1
secondOrder = false
optimize = true
quadAllowed = false
useSurfaceCurvature = true
fuseEdges = true
checkChartBoundary = false
thickness = 0.001
numberOfLayers = 2
stretchFactor = 1.2
isFacesToIgnore = true

166
config.py
View File

@ -1,166 +0,0 @@
import os, sys
from src import applogger
###
# Paths
##
ROOT = os.getcwd()
sys.path.append(ROOT)
LOG = os.path.join(ROOT, "logs")
BUILD = os.path.join(ROOT, "build")
###
# Logger
##
global logger
logger = applogger.Logger()
###
# Utilities
##
class Parameters:
"""
[
"minSize",
"maxSize",
"growthRate",
"nbSegPerEdge",
"nbSegPerRadius",
"chordalErrorEnabled",
"chordalError",
"secondOrder",
"optimize",
"quadAllowed",
"useSurfaceCurvature",
"fuseEdges",
"checkChartBoundary"
]
"""
def __init__(self, **kwargs):
for (k, v) in kwargs.items():
setattr(self, k, v)
class ViscousLayers(Parameters):
"""
[
"thickness",
"numberOfLayers",
"stretchFactor",
"isFacesToIgnore",
"facesToIgnore",
"extrusionMethod"
]
"""
pass
###
# Project variables
##
structures = [
#"simple",
"bodyCentered",
"faceCentered"
]
class simple:
theta = [c * 0.01 for c in range(1, 28 + 1)]
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
fillet = True
fineness = 3
parameters = Parameters(
minSize = 0.01,
maxSize = 0.1,
growthRate = 0.5,
nbSegPerEdge = 0.5,
nbSegPerRadius = 0.5,
chordalErrorEnabled = False,
chordalError = -1,
secondOrder = False,
optimize = True,
quadAllowed = False,
useSurfaceCurvature = True,
fuseEdges = True,
checkChartBoundary = False
)
viscousLayers = ViscousLayers(
thickness = 0.005, # 0.01, 0.005 for 0.28, 0.01 for prism
numberOfLayers = 2,
stretchFactor = 1.2,
isFacesToIgnore = True,
facesToIgnore = None,
extrusionMethod = None
)
class bodyCentered:
theta = [c * 0.01 for c in range(1, 18 + 1)]
directions = [
[1, 0, 0],
[0, 0, 1],
[1, 1, 1]
]
fillet = True
fineness = 3
parameters = Parameters(
minSize = 0.005,
maxSize = 0.05,
growthRate = 0.5,
nbSegPerEdge = 0.5,
nbSegPerRadius = 0.5,
chordalErrorEnabled = False,
chordalError = -1,
secondOrder = False,
optimize = True,
quadAllowed = False,
useSurfaceCurvature = True,
fuseEdges = True,
checkChartBoundary = False
)
viscousLayers = ViscousLayers(
thickness = 0.005,
numberOfLayers = 2,
stretchFactor = 1.2,
isFacesToIgnore = True,
facesToIgnore = None,
extrusionMethod = None
)
class faceCentered:
theta = [0.06, 0.13] #[c * 0.01 for c in range(1, 13 + 1)]
directions = [
#[1, 0, 0],
#[0, 0, 1],
[1, 1, 1]
]
fillet = True
fineness = 3
parameters = Parameters(
minSize = 0.005,
maxSize = 0.05,
growthRate = 0.5,
nbSegPerEdge = 0.5,
nbSegPerRadius = 0.5,
chordalErrorEnabled = False,
chordalError = -1,
secondOrder = False,
optimize = True,
quadAllowed = False,
useSurfaceCurvature = True,
fuseEdges = True,
checkChartBoundary = False
)
viscousLayers = ViscousLayers(
thickness = 0.001, # Failing on 0.13-111
numberOfLayers = 2,
stretchFactor = 1.2,
isFacesToIgnore = True,
facesToIgnore = None,
extrusionMethod = None
)

View File

@ -1,30 +1,8 @@
## relaxationFactors
| p | Iterations | U | - createPatch
| --- | --- | --- | - polyDualMesh
| 0.05 | 830 | 0.2 | - fvSolution: cellLimited
| 0.1 | 678 | 0.2 | - collapseDict: collapseEdges
| 0.2 | 505 | 0.2 | - renumberMesh
| 0.3 | 305 | 0.3 | - processorField
| 0.3 | 236 | 0.4 |
| 0.3 | 181 | 0.5 |
-->
## SIMPLE.residualControl
| p | U | Iterations |
| --- | --- | --- |
| 1e-4 | 1e-4 | 338 |
| 1e-5 | 1e-5 | 499 |
##
```math
U = \frac{\Delta p L}{mu} \frac{9}{256} (\sqrt 2 - \frac{1}{1 - \alpha})^2
L = 2 R_0
R_0 = 1
scale = 1e-5
mu = 1e-3
\Delta p = 1
```

32
openfoam/__init__.py Normal file
View File

@ -0,0 +1,32 @@
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 foamVersion, foamClean
__all__ = [
# meshConversion
"ideasUnvToFoam",
# meshManipulation
"createPatch",
"transformPoints",
"checkMesh",
"renumberMesh",
# miscellaneous
"foamDictionary",
# parallelProcessing
"decomposePar",
# solvers
"potentialFoam",
"simpleFoam",
# utils
"foamVersion",
"foamClean"
]

47
openfoam/application.py Normal file
View File

@ -0,0 +1,47 @@
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, p.returncode

View File

@ -0,0 +1,5 @@
from .application import application
def ideasUnvToFoam(mesh: str, case: str = None) -> (str, int):
return application("ideasUnvToFoam", mesh, case = case, stderr = True)

View File

@ -0,0 +1,38 @@
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: tuple, case: str = None):
scale_ = "{}".format(scale).replace(",", "")
application("transformPoints", "-scale", scale_, case = case, stderr = True)
def checkMesh(case: str = None):
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
def renumberMesh(case: str = None):
application("renumberMesh", "-parallel", "-overwrite", useMPI = True, case = case, stderr = True)

10
openfoam/miscellaneous.py Normal file
View File

@ -0,0 +1,10 @@
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

@ -0,0 +1,5 @@
from .application import application
def decomposePar(case: str = None):
application("decomposePar", case = case, stderr = True)

19
openfoam/solvers.py Normal file
View File

@ -0,0 +1,19 @@
from .application import application
import re
def potentialFoam(case: str = None):
application("potentialFoam", "-parallel", useMPI = True, case = case, stderr = True)
def simpleFoam(case: str = None):
_, returncode = application("simpleFoam", "-parallel", useMPI = True, 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 returncode, out

16
openfoam/utils.py Normal file
View File

@ -0,0 +1,16 @@
import os
import shutil
def foamVersion() -> str:
return "OpenFOAM-{}".format(os.environ["WM_PROJECT_VERSION"])
def foamClean(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 ""
for d in rmDirs:
if os.path.exists(os.path.join(path, d)):
shutil.rmtree(os.path.join(path, d))

3
requirements.txt Normal file
View File

@ -0,0 +1,3 @@
numpy
pyquaternion
toml

24
salomepl/README.md Normal file
View File

@ -0,0 +1,24 @@
## netgen parameters
minSize
maxSize
growthRate
nbSegPerEdge
nbSegPerRadius
chordalErrorEnabled
chordalError
secondOrder
optimize
quadAllowed
useSurfaceCurvature
fuseEdges
checkChartBoundary
## viscous layers parameters
thickness
numberOfLayers
stretchFactor
isFacesToIgnore
facesToIgnore
extrusionMethod

0
salomepl/__init__.py Normal file
View File

View File

@ -1,5 +1,5 @@
import salome #import salome
salome.salome_init() #salome.salome_init()
import GEOM import GEOM
from salome.geom import geomBuilder from salome.geom import geomBuilder
@ -284,4 +284,3 @@ def bodyCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1,
groups.append(wall) groups.append(wall)
return shape, groups return shape, groups

View File

@ -1,5 +1,5 @@
import salome #import salome
salome.salome_init() #salome.salome_init()
import GEOM import GEOM
from salome.geom import geomBuilder from salome.geom import geomBuilder
@ -281,4 +281,3 @@ def faceCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1,
groups.append(wall) groups.append(wall)
return shape, groups return shape, groups

142
salomepl/genmesh.py Normal file
View File

@ -0,0 +1,142 @@
###
# This file executes inside salome environment
#
# salome starts at user home directory
##
import os, sys
import math
import salome
# get project path from args
ROOT = sys.argv[6]
sys.path.append(ROOT)
# site-packages from virtual env
sys.path.append(os.path.join(ROOT, "env/lib/python3.9/site-packages"))
import toml
import logging
from anisotropy.utils import struct
CONFIG = os.path.join(ROOT, "conf/config.toml")
config = struct(toml.load(CONFIG))
LOG = os.path.join(ROOT, "logs")
logging.basicConfig(
level = logging.INFO,
format = config.logger.format,
handlers = [
logging.StreamHandler(),
logging.FileHandler(f"{ LOG }/{ config.logger.name }.log")
]
)
logger = logging.getLogger(config.logger.name)
from salomepl.simple import simpleCubic, simpleHexagonalPrism
from salomepl.faceCentered import faceCenteredCubic, faceCenteredHexagonalPrism
from salomepl.bodyCentered import bodyCenteredCubic, bodyCenteredHexagonalPrism
from salomepl.geometry import getGeom
from salomepl.mesh import smeshBuilder, meshCreate, meshCompute, meshStats, meshExport
def main():
stype = str(sys.argv[1])
theta = float(sys.argv[2])
fillet = int(sys.argv[3])
flowdirection = [int(coord) for coord in sys.argv[4]]
export = str(sys.argv[5])
genmesh(stype, theta, fillet, flowdirection, export)
def genmesh(stype, theta, fillet, direction, export):
logger.info("""genMesh:
structure type:\t{}
coefficient:\t{}
fillet:\t{}
flow direction:\t{}
export path:\t{}""".format(stype, theta, fillet, direction, export))
params = (theta, fillet, direction)
salome.salome_init()
###
# Structure and mesh configurations
##
if stype == "simple":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = simpleCubic
elif direction == [1, 1, 1]:
structure = simpleHexagonalPrism
#fineness = config.simple.fineness
#parameters = config.simple.parameters
#viscousLayers = config.simple.viscousLayers
meshParameters = config.simple.mesh
elif stype == "faceCentered":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = faceCenteredCubic
elif direction == [1, 1, 1]:
structure = faceCenteredHexagonalPrism
#fineness = config.faceCentered.fineness
#parameters = config.faceCentered.parameters
#viscousLayers = config.faceCentered.viscousLayers
meshParameters = config.faceCentered.mesh
elif stype == "bodyCentered":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = bodyCenteredCubic
elif direction == [1, 1, 1]:
structure = bodyCenteredHexagonalPrism
#fineness = config.bodyCentered.fineness
#parameters = config.bodyCentered.parameters
#viscousLayers = config.bodyCentered.viscousLayers
meshParameters = config.bodyCentered.mesh
###
# Shape
##
geompy = getGeom()
shape, groups = structure(*params)
[length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06)
logger.info("""shape:
edges length:\t{}
surface area:\t{}
volume:\t{}""".format(length, surfaceArea, volume))
###
# Mesh
##
facesToIgnore = []
for group in groups:
if group.GetName() in ["inlet", "outlet"]:
facesToIgnore.append(group)
meshParameters.facesToIgnore = facesToIgnore
meshParameters.extrusionMethod = smeshBuilder.SURF_OFFSET_SMOOTH
mesh = meshCreate(shape, groups, meshParameters) #fineness, parameters, viscousLayers)
meshCompute(mesh)
meshStats(mesh)
meshExport(mesh, export)
salome.salome_close()
if __name__ == "__main__":
main()

View File

@ -1,9 +1,9 @@
import GEOM import GEOM
from salome.geom import geomBuilder from salome.geom import geomBuilder
geompy = geomBuilder.New() geompy = geomBuilder.New()
import math import math
import logging
from pyquaternion import Quaternion from pyquaternion import Quaternion
import numpy as np import numpy as np
@ -236,4 +236,3 @@ def boundaryCreate(gobj, dvec, grains):
return boundary return boundary

View File

@ -2,13 +2,14 @@ import SMESH
from salome.smesh import smeshBuilder from salome.smesh import smeshBuilder
smesh = smeshBuilder.New() smesh = smeshBuilder.New()
from config import logger import logging
logger = logging.getLogger("anisotropy")
def getSmesh(): def getSmesh():
return smesh return smesh
def meshCreate(shape, groups, fineness, parameters, viscousLayers = None): def meshCreate(shape, groups, parameters): #fineness, parameters, viscousLayers = None):
""" """
Creates a mesh from a geometry. Creates a mesh from a geometry.
@ -35,7 +36,7 @@ def meshCreate(shape, groups, fineness, parameters, viscousLayers = None):
3: "Fine", 3: "Fine",
4: "Very fine", 4: "Very fine",
5: "Custom" 5: "Custom"
}[fineness] }[parameters.fineness]
# Mesh # Mesh
mesh = smesh.Mesh(shape) mesh = smesh.Mesh(shape)
@ -45,9 +46,9 @@ def meshCreate(shape, groups, fineness, parameters, viscousLayers = None):
param = netgen.Parameters() param = netgen.Parameters()
param.SetMinSize(parameters.minSize) param.SetMinSize(parameters.minSize)
param.SetMaxSize(parameters.maxSize) param.SetMaxSize(parameters.maxSize)
param.SetFineness(fineness) param.SetFineness(parameters.fineness)
if fineness == 5: if parameters.fineness == 5:
param.SetGrowthRate(parameters.growthRate) param.SetGrowthRate(parameters.growthRate)
param.SetNbSegPerEdge(parameters.nbSegPerEdge) param.SetNbSegPerEdge(parameters.nbSegPerEdge)
param.SetNbSegPerRadius(parameters.nbSegPerRadius) param.SetNbSegPerRadius(parameters.nbSegPerRadius)
@ -93,30 +94,30 @@ def meshCreate(shape, groups, fineness, parameters, viscousLayers = None):
### ###
# Viscous layers # Viscous layers
## ##
if not viscousLayers is None: #if not viscousLayers is None:
vlayer = netgen.ViscousLayers( vlayer = netgen.ViscousLayers(
viscousLayers.thickness, parameters.thickness,
viscousLayers.numberOfLayers, parameters.numberOfLayers,
viscousLayers.stretchFactor, parameters.stretchFactor,
viscousLayers.facesToIgnore, parameters.facesToIgnore,
viscousLayers.isFacesToIgnore, parameters.isFacesToIgnore,
viscousLayers.extrusionMethod parameters.extrusionMethod
) )
logger.info("""meshCreate: logger.info("""meshCreate:
viscous layers: viscous layers:
thickness:\t{} thickness:\t{}
number:\t{} number:\t{}
stretch factor:\t{}""".format( stretch factor:\t{}""".format(
viscousLayers.thickness, parameters.thickness,
viscousLayers.numberOfLayers, parameters.numberOfLayers,
viscousLayers.stretchFactor)) parameters.stretchFactor))
else: #else:
logger.info("""meshCreate: # logger.info("""meshCreate:
viscous layers: disabled""") #viscous layers: disabled""")
return mesh return mesh
@ -191,4 +192,3 @@ def meshExport(mobj, path):
except: except:
logger.error("""meshExport: Cannot export.""") logger.error("""meshExport: Cannot export.""")

View File

@ -1,5 +1,5 @@
import salome #import salome
salome.salome_init() #salome.salome_init()
import GEOM import GEOM
from salome.geom import geomBuilder from salome.geom import geomBuilder
@ -261,4 +261,3 @@ def simpleHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]):
groups.append(wall) groups.append(wall)
return shape, groups return shape, groups

View File

@ -2,8 +2,8 @@
import subprocess import subprocess
import logging import logging
import sys, os import sys, os
import config
from config import logger logger = logging.getLogger()
def hasDesktop() -> bool: def hasDesktop() -> bool:
return salome.sg.hasDesktop() return salome.sg.hasDesktop()
@ -80,4 +80,3 @@ def remote(port, cmd):
stderr = subprocess.PIPE) stderr = subprocess.PIPE)
return p return p

View File

@ -1,117 +0,0 @@
from collections import namedtuple
import os, sys
import logging
from pyquaternion import Quaternion
import math
import salome
sys.path.append(sys.argv[6])
import config
from config import logger
#from src import applogger
from simple import simpleCubic, simpleHexagonalPrism
from faceCentered import faceCenteredCubic, faceCenteredHexagonalPrism
from bodyCentered import bodyCenteredCubic, bodyCenteredHexagonalPrism
from src import geometry_utils
from src import mesh_utils
def main():
stype = str(sys.argv[1])
theta = float(sys.argv[2])
fillet = int(sys.argv[3])
flowdirection = [int(coord) for coord in sys.argv[4]]
export = str(sys.argv[5])
genMesh(stype, theta, fillet, flowdirection, export)
def genMesh(stype, theta, fillet, direction, export):
logger.info("""genMesh:
structure type:\t{}
coefficient:\t{}
fillet:\t{}
flow direction:\t{}
export path:\t{}""".format(stype, theta, fillet, direction, export))
params = (theta, fillet, direction)
salome.salome_init()
###
# Structure and mesh configurations
##
if stype == "simple":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = simpleCubic
elif direction == [1, 1, 1]:
structure = simpleHexagonalPrism
fineness = config.simple.fineness
parameters = config.simple.parameters
viscousLayers = config.simple.viscousLayers
elif stype == "faceCentered":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = faceCenteredCubic
elif direction == [1, 1, 1]:
structure = faceCenteredHexagonalPrism
fineness = config.faceCentered.fineness
parameters = config.faceCentered.parameters
viscousLayers = config.faceCentered.viscousLayers
elif stype == "bodyCentered":
if direction in [[1, 0, 0], [0, 0, 1]]:
structure = bodyCenteredCubic
elif direction == [1, 1, 1]:
structure = bodyCenteredHexagonalPrism
fineness = config.bodyCentered.fineness
parameters = config.bodyCentered.parameters
viscousLayers = config.bodyCentered.viscousLayers
###
# Shape
##
geompy = geometry_utils.getGeom()
shape, groups = structure(*params)
[length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06)
logger.info("""shape:
edges length:\t{}
surface area:\t{}
volume:\t{}""".format(length, surfaceArea, volume))
###
# Mesh
##
facesToIgnore = []
for group in groups:
if group.GetName() in ["inlet", "outlet"]:
facesToIgnore.append(group)
viscousLayers.facesToIgnore = facesToIgnore
viscousLayers.extrusionMethod = mesh_utils.smeshBuilder.SURF_OFFSET_SMOOTH
mesh = mesh_utils.meshCreate(shape, groups, fineness, parameters, viscousLayers)
mesh_utils.meshCompute(mesh)
mesh_utils.meshStats(mesh)
mesh_utils.meshExport(mesh, export)
salome.salome_close()
if __name__ == "__main__":
main()

View File

@ -1,103 +0,0 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import math
from . import geometry_utils
import GEOM
geompy = geometry_utils.getGeom()
class StructuredGrains:
def __init__(self, radius, stackAngle, theta, layers):
self.pos = [0, 0, 0]
self.angle = [0, 0, 0]
self.radius = radius
self.theta = theta
self.layers = layers
# Parameters and dependencies
R = self.radius / (1 - self.theta)
C1 = 0.8 #fillet[0]
C2 = 0.4 #fillet[1]
self.theta1 = 0.01
self.theta2 = 0.28
Cf = C1 + (C2 - C1) / (self.theta2 - self.theta1) * (self.theta - self.theta1)
R_fillet = Cf * (self.radius * math.sqrt(2) - R)
###
stackang = [
0.5 * math.pi - stackAngle[0],
0.5 * math.pi - stackAngle[1],
0.5 * math.pi - stackAngle[2]
]
xvec = geompy.MakeVector(
geompy.MakeVertex(0, 0, 0),
geompy.MakeVertex(1, 0, 0))
yvec = geometry_utils.rotate(xvec, [0.5 * math.pi, 0, 0])
zvec = geometry_utils.rotate(xvec, [0, 0.5 * math.pi, 0])
grain = geompy.MakeSpherePntR(geompy.MakeVertex(pos[0], pos[1], pos[2]), R)
xstack = geompy.MakeMultiTranslation1D(grain, xvec, 2 * self.radius, self.layers[0])
ystack = geompy.MakeMultiTranslation1D(xgrain, yvec, 2 * self.radius, self.layers[1])
zstack = geompy.MakeMultiTranslation1D(ygrain, zvec, 2 * self.radius, self.layers[2])
# Correct position to zero
stack = geompy.MakeTranslation(zstack, -2 * self.radius, 0, 0)
self.geometry = geompy.ExtractShapes(stack, geompy.ShapeType["SOLID"], True)
self.geometry = geompy.MakeFuseList(self.geometry, False, False)
if not R_fillet == 0:
self.geometry = geompy.MakeFilletAll(self.geometry, R_fillet)
class AnisotropeCubic:
def __init__(self, scale, grains, style):
self.pos = [0, 0, 0]
self.angle = [0, 0, 0]
self.scale = scale
self.grains = grains
# Bounding box
if style == 0:
# Square
profile = (
geompy.Sketcher3D()
.addPointAbsolute(0, 0, 0)
.addPointAbsolute(0, 0, self.scale[2])
.addPointAbsolute(0, self.scale[1], self.scale[2])
.addPointAbsolute(0, self.scale[1], 0)
.addPointAbsolute(0, 0, 0)
)
face = geompy.MakeFaceWires([profile.wire()], 1)
elif style == 1:
# Rombus
profile = (
geompy.Sketcher3D()
.addPointAbsolute(self.scale[0], 0.5 * self.scale[1], 0)
.addPointAbsolute(0.5 * self.scale[0], 0, 0.5 * self.scale[2])
.addPointAbsolute(0, 0.5 * self.scale[1], self.scale[2])
.addPointAbsolute(0.5 * self.scale[0], self.scale[1], 0.5 * self.scale[2])
.addPointAbsolute(self.scale[0], 0.5 * self.scale[1], 0)
)
face = geompy.MakeFaceWires([profile.wire()], 1)
face = geompy.MakeTranslation(face,
0.5 * self.scale[1], 0, 0)
self.boundingbox = geompy.MakePrismVecH(face,
geompy.MakeVectorDXDYDZ(1, 0, 0),
self.scale[0])
# Geometry
self.geometry = geompy.MakeCutList(box, [self.grains], True)

View File

@ -1,41 +0,0 @@
import logging
import config
class Logger():
def __init__(self):
logging.basicConfig(
level = logging.INFO,
format = "%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler(f"{config.LOG}/anisotrope.log")
]
)
self.logger = logging.getLogger("anisotrope")
self.warnings = 0
self.errors = 0
self.criticals = 0
self.exceptions = 0
def info(self, *args):
self.logger.info(*args)
def warning(self, *args):
self.warnings += 1
self.logger.warning(*args)
def error(self, *args):
self.errors += 1
self.logger.error(*args)
def critical(self, *args):
self.criticals += 1
self.logger.critical(*args)
def exception(self, *args):
self.exceptions += 1
self.logger.exception(*args)
def fancyline(self):
self.logger.info("-" * 80)

View File

@ -1,126 +0,0 @@
import os, sys, shutil
import subprocess
import logging
import time
import re
from datetime import timedelta
from config import logger
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, p.returncode
def foamVersion() -> str:
return "OpenFOAM-{}".format(os.environ["WM_PROJECT_VERSION"])
def foamClean(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 ""
for d in rmDirs:
if os.path.exists(os.path.join(path, d)):
shutil.rmtree(os.path.join(path, d))
def ideasUnvToFoam(mesh: str, case: str = None) -> (str, int):
return application("ideasUnvToFoam", mesh, case = case, stderr = True)
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: tuple, case: str = None):
scale_ = "{}".format(scale).replace(",", "")
application("transformPoints", "-scale", scale_, case = case, stderr = True)
def checkMesh(case: str = None):
application("checkMesh", "-allGeometry", "-allTopology", case = case, stderr = True)
with open("checkMesh.log", "r") as io:
warnings = []
for line in io:
if re.search("\*\*\*", line):
warnings.append(line.replace("***", "").strip())
if warnings:
logger.warning("checkMesh:\n\t{}".format("\n\t".join(warnings)))
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)
def decomposePar(case: str = None):
application("decomposePar", case = case, stderr = True)
def renumberMesh(case: str = None):
application("renumberMesh", "-parallel", "-overwrite", useMPI = True, case = case, stderr = True)
def potentialFoam(case: str = None):
application("potentialFoam", "-parallel", useMPI = True, case = case, stderr = True)
def simpleFoam(case: str = None):
_, returncode = application("simpleFoam", "-parallel", useMPI = True, case = case, stderr = True)
with open("simpleFoam.log", "r") as io:
for line in io:
if re.search("solution converged", line):
logger.info("simpleFoam:\n\t{}".format(line.strip()))
return returncode

View File

@ -1,77 +0,0 @@
from multiprocessing import Queue, Process, cpu_count
import socket
import logging
def queue(cmd, qin, qout, *args):
while True:
# Get item from the queue
pos, var = qin.get()
# 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()
logging.info("cpu count: {}".format(np))
logging.info("var: {}".format(var))
logging.info("cmd: {}".format(cmd))
# 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

View File

@ -1,80 +0,0 @@
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
geompy = geomBuilder.New()
import SMESH, SALOMEDS
from salome.smesh import smeshBuilder
smesh = smeshBuilder.New()
import math
axes = [
geompy.MakeVectorDXDYDZ(1, 0, 0),
geompy.MakeVectorDXDYDZ(0, 1, 0),
geompy.MakeVectorDXDYDZ(0, 0, 1)
]
vtx = [
geompy.MakeVertex(1 / math.sqrt(2), -1 / math.sqrt(2), 0.5),
geompy.MakeVertex(-1 / math.sqrt(2), 1 / math.sqrt(2), -0.5),
geompy.MakeVertex(-1, -1, -0.5),
geompy.MakeVertex(-0.5, -0.5, 0)
]
box = geompy.MakeBoxTwoPnt(vtx[0], vtx[1])
box = geompy.MakeRotation(box, axes[2], 45 * math.pi / 180.0)
#alpha = [ x * 0.01 for x in range(1, 28 + 1) ]
alpha=[0.15]
#for coef in alpha:
spheres = geompy.MakeSpherePntR(vtx[2], math.sqrt(3) / 4 / (1 - alpha[0]))
spheres = geompy.MakeMultiTranslation2D(spheres, axes[0], 1, 3, axes[1], 1, 3)
spheres = geompy.MakeMultiTranslation1D(spheres, axes[2], 1, 2)
spheres2 = geompy.MakeSpherePntR(vtx[3], math.sqrt(3) / 4 / (1 - alpha[0]))
spheres2 = geompy.MakeMultiTranslation2D(spheres2, axes[0], 1, 3, axes[1], 1, 3)
PoreBC = geompy.MakeCutList(box, [spheres, spheres2], False)
geompy.addToStudy(PoreBC, 'PoreBC')
box2 = geompy.MakeBoxTwoPnt(vtx[0], geompy.MakeVertex(0, 0, 0))
box2 = geompy.MakeRotation(box2, axes[2], 45 * math.pi / 180.0)
Segment1_8 = geompy.MakeCommonList([PoreBC, box2], True)
geompy.addToStudy(Segment1_8, 'Segment1_8')
mesh = smesh.Mesh(PoreBC)
netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
param = netgen.Parameters()
param.SetSecondOrder( 0 )
param.SetOptimize( 1 )
param.SetChordalError( -1 )
param.SetChordalErrorEnabled( 0 )
param.SetUseSurfaceCurvature( 1 )
param.SetFuseEdges( 1 )
param.SetCheckChartBoundary( 0 )
param.SetMinSize( 0.01 )
param.SetMaxSize( 0.02 )
param.SetFineness( 5 )
param.SetGrowthRate( 0.1 )
param.SetNbSegPerEdge( 5 )
param.SetNbSegPerRadius( 10 )
param.SetQuadAllowed( 1 )
#vlayer = netgen.ViscousLayers(0.05, 3, 1.5, [15, 29, 54], 1, smeshBuilder.SURF_OFFSET_SMOOTH)
isDone = mesh.Compute()
if not isDone:
print("Mesh is not computed")
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

View File

@ -1,80 +0,0 @@
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
geompy = geomBuilder.New()
import SMESH, SALOMEDS
from salome.smesh import smeshBuilder
smesh = smeshBuilder.New()
import math
axes = [
geompy.MakeVectorDXDYDZ(1, 0, 0),
geompy.MakeVectorDXDYDZ(0, 1, 0),
geompy.MakeVectorDXDYDZ(0, 0, 1),
geompy.MakeVectorDXDYDZ(1, 1, 0),
geompy.MakeVectorDXDYDZ(1, -1, 0)
]
vtx = [
geompy.MakeVertex(0, 0, -0.5),
geompy.MakeVertex(1 / math.sqrt(2), 1 / math.sqrt(2), 0.5),
geompy.MakeVertex(0.5, 0, 0),
geompy.MakeVertex(0.5 / math.sqrt(2), 0.5 / math.sqrt(2), 0.5)
]
box = geompy.MakeBoxTwoPnt(vtx[0], vtx[1])
box = geompy.MakeRotation(box, axes[2], -45 * math.pi / 180.0)
#alpha = [ x * 0.01 for x in range(1, 28 + 1) ]
alpha=[0.08]
#for coef in alpha:
spheres = geompy.MakeSpherePntR(vtx[0], math.sqrt(2) / 4 / (1 - alpha[0]))
spheres = geompy.MakeMultiTranslation2D(spheres, axes[3], 1 / math.sqrt(2), 2, axes[4], 1 / math.sqrt(2), 2)
spheres = geompy.MakeMultiTranslation1D(spheres, axes[2], 1, 2)
sphere2 = geompy.MakeSpherePntR(vtx[2], math.sqrt(2) / 4 / (1 - alpha[0]))
PoreFC = geompy.MakeCutList(box, [spheres, sphere2], True)
geompy.addToStudy(PoreFC, 'PoreFC')
box2 = geompy.MakeBoxTwoPnt(geompy.MakeVertex(0, 0, 0), vtx[3])
box2 = geompy.MakeRotation(box2, axes[2], -45 * math.pi / 180.0)
Segment1_8 = geompy.MakeCommonList([PoreFC, box2], True)
geompy.addToStudy(Segment1_8, 'Segment1_8')
mesh = smesh.Mesh(PoreFC)
netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
param = netgen.Parameters()
param.SetSecondOrder( 0 )
param.SetOptimize( 1 )
param.SetChordalError( -1 )
param.SetChordalErrorEnabled( 0 )
param.SetUseSurfaceCurvature( 1 )
param.SetFuseEdges( 1 )
param.SetCheckChartBoundary( 0 )
param.SetMinSize( 0.01 )
param.SetMaxSize( 0.02 )
param.SetFineness( 5 )
param.SetGrowthRate( 0.1 )
param.SetNbSegPerEdge( 5 )
param.SetNbSegPerRadius( 10 )
param.SetQuadAllowed( 1 )
#vlayer = netgen.ViscousLayers(0.05, 3, 1.5, [15, 29, 54], 1, smeshBuilder.SURF_OFFSET_SMOOTH)
isDone = mesh.Compute()
if not isDone:
print("Mesh is not computed")
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

File diff suppressed because it is too large Load Diff

View File

@ -1,502 +0,0 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import salome, GEOM, SMESH, SALOMEDS
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math
import os, sys
import logging
import time
from datetime import timedelta
class faceCenteredCubic:
def __init__(self, name = None):
self.name = name if type(name) != None else "faceCenteredCubic"
self.geometry = None
self.geometrybbox = None
self.mesh = None
self.boundary = None
self.rombus = None
self.rombusbbox = None
self.spheres = None
salome.salome_init()
def geometryCreate(self, alpha, fillet):
"""
Create the simple cubic geometry.
Parameters:
alpha (float): Sphere intersection parameter which used for cutting spheres from box.
Radius = R0 / (1 - alpha)
Should be from 0.01 to 0.28
fillet (list): Fillet coefficient.
[fillet1, fillet2]
0 <= fillet <= 1
if fillet = [0, 0] then R_fillet = 0
Returns:
Configured geometry.
"""
geompy = geomBuilder.New()
# Parameters
R0 = math.sqrt(2) / 4
R = R0 / (1 - alpha)
C1 = fillet[0]
C2 = fillet[1]
alpha1 = 0.01
alpha2 = 0.28
Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
R_fillet = Cf * (R0 * math.sqrt(2) - R)
logging.info("geometryCreate: alpha = {}".format(alpha))
logging.info("geometryCreate: R = {}".format(R))
logging.info("geometryCreate: R_fillet = {}".format(R_fillet))
# xyz axes
axes = [
geompy.MakeVectorDXDYDZ(1, 0, 0),
geompy.MakeVectorDXDYDZ(0, 1, 0),
geompy.MakeVectorDXDYDZ(0, 0, 1)
]
# Main box
size = [1 / math.sqrt(2), 1 / math.sqrt(2), 1]
angle = [0, 0, 0]
pos = [0, 0, 0]
box = geompy.MakeBoxDXDYDZ(size[0], size[1], size[2])
for n in range(3):
box = geompy.MakeRotation(box, axes[n], angle[n] * math.pi / 180.0)
box = geompy.MakeTranslation(box, pos[0], pos[1], pos[2])
#[x, y, z, _, _, _, _, _, _] = geompy.GetPosition(box)
#pos = [x, y, z]
# Spheres for cutting
sphere = geompy.MakeSpherePntR(geompy.MakeVertex(pos[0], pos[1], pos[2]), R)
sphere = geompy.MakeMultiTranslation2D(sphere, None, 2 * R0, 3, None, 2 * R0, 3)
sphere2 = geompy.MakeTranslation(sphere, 0, 0, 2 * math.sqrt(2) * R0)
sphere3 = geompy.MakeTranslation(sphere, R0, R0, 0.5)
sphere = geompy.ExtractShapes(sphere, geompy.ShapeType["SOLID"], True)
sphere2 = geompy.ExtractShapes(sphere2, geompy.ShapeType["SOLID"], True)
sphere3 = geompy.ExtractShapes(sphere3, geompy.ShapeType["SOLID"], True)
sphere = geompy.MakeFuseList(sphere + sphere2 + sphere3, False, False)
if not R_fillet == 0:
sphere = geompy.MakeFilletAll(sphere, R_fillet)
self.spheres = sphere
# Main geometry and bounding box
# geompy.RemoveExtraEdges(obj, True)
self.geometry = geompy.MakeCutList(box, [sphere], True)
self.geometrybbox = box
# Rombus and bounding box
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, size[2])
sk.addPointsAbsolute(size[2] / 2, 0, size[2] / 2)
sk.addPointsAbsolute(size[2] / 2, size[2] / 2, 0)
sk.addPointsAbsolute(0, size[2] / 2, size[2] / 2)
sk.addPointsAbsolute(0, 0, size[2])
face = geompy.MakeFaceWires([sk.wire()], 1)
rombusbbox = geompy.MakePrismVecH(face, geompy.MakeVectorDXDYDZ(1, 1, 0), size[0])
rombusbbox = geompy.MakeRotation(rombusbbox, axes[2], 45 * math.pi / 180.0)
rombusbbox = geompy.MakeTranslation(rombusbbox, size[0], 0, 0)
self.rombus = geompy.MakeCutList(rombusbbox, [sphere], True)
self.rombusbbox = rombusbbox
# Change position
self.spheres = geompy.MakeRotation(self.spheres, axes[2], -45 * math.pi / 180.0)
self.spheres = geompy.MakeTranslation(self.spheres, 0, 0.5, 0)
self.geometry = geompy.MakeRotation(self.geometry, axes[2], -45 * math.pi / 180.0)
self.geometry = geompy.MakeTranslation(self.geometry, 0, 0.5, 0)
self.geometrybbox = geompy.MakeRotation(self.geometrybbox, axes[2], -45 * math.pi / 180.0)
self.geometrybbox = geompy.MakeTranslation(self.geometrybbox, 0, 0.5, 0)
self.rombus = geompy.MakeRotation(self.rombus, axes[2], -45 * math.pi / 180.0)
self.rombus = geompy.MakeTranslation(self.rombus, 0, 0.5, 0)
self.rombusbbox = geompy.MakeRotation(self.rombusbbox, axes[2], -45 * math.pi / 180.0)
self.rombusbbox = geompy.MakeTranslation(self.rombusbbox, 0, 0.5, 0)
geompy.addToStudy(self.spheres, "spheres")
geompy.addToStudy(self.geometry, self.name)
geompy.addToStudy(self.rombus, "rombus")
geompy.addToStudy(rombusbbox, "rombusbbox")
return self.geometry
def boundaryCreate(self, direction):
"""
Create the boundary faces from the geometry.
Parameters:
direction (str): Direction of the flow.
'001' for the flow with normal vector (0, 0, 1) to face.
'100' for the flow with normal vector (1, 0, 0) to face.
'111' for (1, 1, 1)
Returns:
boundary (dict):
{
"inlet": <GEOM._objref_GEOM_Object>,
"outlet": <GEOM._objref_GEOM_Object>,
"symetryPlane": <GEOM._objref_GEOM_Object>,
"wall": <GEOM._objref_GEOM_Object>
}
"""
geompy = geomBuilder.New()
rot = [0, 0, 45]
buffergeometry = self.geometry
if direction == "001":
[x, y, z, _, _, _, _, _, _] = geompy.GetPosition(self.geometry)
center = geompy.MakeVertex(x, y, z)
norm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, 0, 0, 1))
bnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center,
-math.cos((90 + rot[2]) * math.pi / 180.0),
math.sin((90 + rot[2]) * math.pi / 180.0), 0))
vnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center,
-math.cos((0 + rot[2]) * math.pi / 180.0),
math.sin((0 + rot[2]) * math.pi / 180.0), 0))
elif direction == "100":
[x, y, z, _, _, _, _, _, _] = geompy.GetPosition(self.geometry)
center = geompy.MakeVertex(x, y, z)
norm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center,
-math.cos((90 + rot[2]) * math.pi / 180.0),
math.sin((90 + rot[2]) * math.pi / 180.0), 0))
bnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, 0, 0, 1))
vnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center,
-math.cos((0 + rot[2]) * math.pi / 180.0),
math.sin((0 + rot[2]) * math.pi / 180.0), 0))
elif direction == "111":
self.geometry = self.rombus
[x, y, z, _, _, _, _, _, _] = geompy.GetPosition(self.geometry)
center = geompy.MakeVertex(x, y, z)
norm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, 1, 1, 1))
#-math.cos((90 + rot[2]) * math.pi / 180.0),
#math.sin((90 + rot[2]) * math.pi / 180.0), math.sqrt(2) / 2))
bnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, 1, -1, 1))
# -math.cos((90 + rot[2]) * math.pi / 180.0),
# math.sin((90 + rot[2]) * math.pi / 180.0), 0))
vnorm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, -1, 1, 1))
#-math.cos((0 + rot[2]) * math.pi / 180.0),
#math.sin((0 + rot[2]) * math.pi / 180.0), 0))
logging.info("boundaryCreate: direction = {}".format(direction))
geompy.addToStudy(norm, "normalvector")
geompy.addToStudy(bnorm, "bnorm")
geompy.addToStudy(vnorm, "vnorm")
if direction == "111":
box = self.rombus
else:
box = self.geometrybbox
planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True)
inletplane = []
outletplane = []
hplanes = []
fwplanes = []
bwplanes = []
lplanes = []
rplanes = []
for plane in planes:
planeNorm = geompy.GetNormal(plane)
angle = round(abs(geompy.GetAngle(planeNorm, norm)), 0)
if angle == 0:
outletplane.append(plane)
elif angle == 180:
inletplane.append(plane)
elif direction == "111" and (angle == 109 or angle == 71):
#hplanes.append(plane)
bangle = round(abs(geompy.GetAngle(planeNorm, bnorm)), 0)
#logging.info("bangle = {}".format(bangle))
if bangle == 0:
fwplanes.append(plane)
elif bangle == 180:
bwplanes.append(plane)
vangle = round(abs(geompy.GetAngle(planeNorm, vnorm)), 0)
#logging.info("vangle = {}".format(vangle))
if vangle == 0:
lplanes.append(plane)
elif vangle == 180:
rplanes.append(plane)
elif direction == "100" or direction == "001":
if angle == 90:
#hplanes.append(plane)
bangle = round(abs(geompy.GetAngle(planeNorm, bnorm)), 0)
#logging.info("bangle = {}".format(bangle))
if bangle == 0:
fwplanes.append(plane)
elif bangle == 180:
bwplanes.append(plane)
vangle = round(abs(geompy.GetAngle(planeNorm, vnorm)), 0)
#logging.info("vangle = {}".format(vangle))
if vangle == 0:
lplanes.append(plane)
elif vangle == 180:
rplanes.append(plane)
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
logging.info("boundaryCreate: inletplanes = {}, outletplanes = {}, hplanes = {}".format(
len(inletplane), len(outletplane), len(hplanes)))
logging.info("boundaryCreate: fwplanes = {}, bwplanes = {}, lplanes = {}, rplanes = {}".format(
len(fwplanes), len(bwplanes), len(lplanes), len(rplanes)))
def createGroup(planelist, name):
gr = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], name)
grcomp = geompy.MakeCompound(planelist)
grcut = geompy.MakeCutList(grcomp, [self.spheres], True)
gip = geompy.GetInPlace(self.geometry, grcut, True)
faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
geompy.UnionList(gr, faces)
return gr
# Main groups
inlet = createGroup(inletplane, "inlet")
outlet = createGroup(outletplane, "outlet")
# Symetry planes
symetryPlaneFW = createGroup(fwplanes, "symetryPlaneFW")
symetryPlaneBW = createGroup(bwplanes, "symetryPlaneBW")
symetryPlaneL = createGroup(lplanes, "symetryPlaneL")
symetryPlaneR = createGroup(rplanes, "symetryPlaneR")
# Wall
allgroup = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"])
faces = geompy.SubShapeAllIDs(self.geometry, geompy.ShapeType["FACE"])
geompy.UnionIDs(allgroup, faces)
wall = geompy.CutListOfGroups([allgroup],
[inlet, outlet, symetryPlaneFW, symetryPlaneBW, symetryPlaneL, symetryPlaneR], "wall")
self.boundary = {
"inlet": inlet,
"outlet": outlet,
"symetryPlaneFW": symetryPlaneFW,
"symetryPlaneBW": symetryPlaneBW,
"symetryPlaneL": symetryPlaneL,
"symetryPlaneR": symetryPlaneR,
"wall": wall
}
return self.boundary
def meshCreate(self, fineness, viscousLayers=None):
"""
Creates a mesh from a geometry.
Parameters:
fineness (int): Fineness of mesh.
0 - Very coarse,
1 - Coarse,
2 - Moderate,
3 - Fine,
4 - Very fine.
viscousLayers (dict or None): Defines viscous layers for mesh.
By default, inlets and outlets specified without layers.
{
"thickness": float,
"number": int,
"stretch": float
}
Returns:
Configured instance of class <SMESH.SMESH_Mesh>, containig the parameters and boundary groups.
"""
smesh = smeshBuilder.New()
Fineness = {
0: "Very coarse",
1: "Coarse",
2: "Moderate",
3: "Fine",
4: "Very fine"
}[fineness]
logging.info("meshCreate: mesh fineness - {}".format(Fineness))
mesh = smesh.Mesh(self.geometry)
netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
param = netgen.Parameters()
param.SetSecondOrder( 0 )
param.SetOptimize( 1 )
param.SetChordalError( -1 )
param.SetChordalErrorEnabled( 0 )
param.SetUseSurfaceCurvature( 1 )
param.SetFuseEdges( 1 )
param.SetCheckChartBoundary( 0 )
param.SetMinSize( 0.001 )
param.SetMaxSize( 0.1 )
param.SetFineness(fineness)
#param.SetGrowthRate( 0.1 )
#param.SetNbSegPerEdge( 5 )
#param.SetNbSegPerRadius( 10 )
param.SetQuadAllowed( 0 )
if not viscousLayers is None:
logging.info("meshCreate: viscous layers params - thickness = {}, number = {}, stretch factor = {}".format(
viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"]))
vlayer = netgen.ViscousLayers(
viscousLayers["thickness"],
viscousLayers["number"],
viscousLayers["stretch"],
[self.boundary["inlet"], self.boundary["outlet"]],
1, smeshBuilder.NODE_OFFSET)
else:
logging.info("meshCreate: viscous layers are disabled")
for name, boundary in self.boundary.items():
mesh.GroupOnGeom(boundary, "{}_".format(name), SMESH.FACE)
self.mesh = mesh
return self.mesh
def meshCompute(self):
"""Compute the mesh."""
status = self.mesh.Compute()
if status:
logging.info("Mesh succesfully computed.")
else:
logging.warning("Mesh is not computed.")
def meshExport(self, path):
"""
Export the mesh in a file in UNV format.
Parameters:
path (string): full path to the expected directory.
"""
exportpath = os.path.join(path, "{}.unv".format(self.name))
try:
self.mesh.ExportUNV(exportpath)
except:
logging.error("Cannot export mesh to '{}'".format(exportpath))
if __name__ == "__main__":
# Arguments
buildpath = str(sys.argv[1])
alpha = float(sys.argv[2])
direction = str(sys.argv[3])
name = "faceCenteredCubic-{}-{}".format(direction, alpha)
# Logger
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler(os.path.join(buildpath, "{}.log".format(name)))
])
start_time = time.monotonic()
# Simple cubic
fcc = faceCenteredCubic(name)
logging.info("Creating the geometry ...")
fcc.geometryCreate(alpha, [0, 0])
logging.info("Extracting boundaries ...")
fcc.boundaryCreate(direction)
logging.info("Creating the mesh ...")
fcc.meshCreate(2, {
"thickness": 0.01,
"number": 2,
"stretch": 1.1
})
fcc.meshCompute()
logging.info("Exporting the mesh ...")
fcc.meshExport(buildpath)
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
logging.info("Done.")
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

View File

@ -1,91 +0,0 @@
import os
import subprocess
import multiprocessing
import logging
import time
from datetime import timedelta
def salome(port, src_path, build_path, coefficient, direction):
"""
Run salome and terminate after.
Parameters:
port (int): Port for the process.
src_path (str): Path to the execution script.
build_path (str): Output path.
"""
logging.info("Starting SALOME on port {}.".format(port))
salomelog = open("{}/salome.log".format(build_path), "a")
subprocess.run([
"salome", "start",
"--port", str(port),
"-t", src_path,
"args:{},{},{}".format(build_path, coefficient, direction)
],
stdout=salomelog,
stderr=salomelog)
logging.info("Terminating SALOME on port {}.".format(port))
subprocess.run([
"salome", "kill",
str(port)
],
stdout=salomelog,
stderr=salomelog)
salomelog.close()
if __name__ == "__main__":
# Get main paths
project = os.getcwd()
src = os.path.join(project, "src")
build = os.path.join(project, "build")
if not os.path.exists(build):
os.makedirs(build)
# Logger
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/genmesh.log".format(build))
])
start_time = time.monotonic()
# Start in parallel
processes = []
structures = ["simpleCubic", "faceCenteredCubic"] #, "bodyCenteredCubic"]
directions = ["001"] #, "100", "111"]
coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ]
port = 2810
for structure in structures:
for direction in directions:
for coefficient in coefficients:
src_path = os.path.join(src, "{}.py".format(structure))
build_path = os.path.join(build,
structure,
"direction-{}".format(direction),
"alpha-{}".format(coefficient))
if not os.path.exists(build_path):
os.makedirs(build_path)
p = multiprocessing.Process(target = salome,
args = (port, src_path, build_path, coefficient, direction))
processes.append(p)
p.start()
logging.info("{} on port {}.".format(p, port))
port += 1
for process in processes:
process.join()
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))

View File

@ -1,119 +0,0 @@
import salome, GEOM, SMESH, SALOMEDS
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math
import os, sys
import logging
import time
from datetime import timedelta
def rotate(gobj, ang):
x = geompy.MakeVectorDXDYDZ(1, 0, 0),
y = geompy.MakeVectorDXDYDZ(0, 1, 0),
z = geompy.MakeVectorDXDYDZ(0, 0, 1)
# yaw
rotated = geompy.MakeRotation(gobj, z, ang[2])
# pitch
rotated = geompy.MakeRotation(rotated, y, ang[1])
# roll
rotated = geompy.MakeRotation(rotated, x, ang[0])
return rotated
def createGroup(gobj, planelist, grains, name):
gr = geompy.CreateGroup(gobj, geompy.ShapeType["FACE"], name)
grcomp = geompy.MakeCompound(planelist)
grcut = geompy.MakeCutList(grcomp, [grains], True)
gip = geompy.GetInPlace(gobj, grcut, True)
faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
geompy.UnionList(gr, faces)
return gr
def boundaryCreate(gobj, dvec, grains):
geompy = geomBuilder.New()
xvec = geompy.MakeVector(
geompy.MakeVertex(0, 0, 0),
geompy.MakeVertex(dvec[0], dvec[1], dvec[2]))
xvec = rotate(dvec, self.angle)
yvec = rotate(xvec, [0.5 * math.pi, 0, 0])
zvec = rotate(xvec, [0, 0.5 * math.pi, 0])
logging.info("boundaryCreate: dvec = {}".format(dvec))
planes = geompy.ExtractShapes(gobj, geompy.ShapeType["FACE"], True)
inletplanes = []
outletplanes = []
uplanes = []
fwplanes = []
bwplanes = []
lplanes = []
rplanes = []
for plane in planes:
nvec = geompy.GetNormal(plane)
xang = geompy.GetAngle(nvec, xvec)
yang = geompy.GetAngle(nvec, yvec)
zang = geompy.GetAngle(nvec, zvec)
if xang == 0:
inletplanes.append(plane)
elif xang == 180:
outletplanes.append(plane)
elif yang == 0:
fwplanes.append(plane)
elif yang == 180:
bwplanes.append(plane)
elif zang == 0:
lplanes.append(plane)
elif zang == 180:
rplanes.append(plane)
logging.info(
"boundaryCreate: inletplanes = {}, outletplanes = {}, hplanes = {}".format(
len(inletplane), len(outletplane), len(hplanes)))
logging.info(
"boundaryCreate: fwplanes = {}, bwplanes = {}, lplanes = {}, rplanes = {}".format(
len(fwplanes), len(bwplanes), len(lplanes), len(rplanes)))
# Main groups
inlet = createGroup(gobj, inletplane, grains, "inlet")
outlet = createGroup(gobj, grains, outletplane, "outlet")
symetryPlaneFW = createGroup(gobj, fwplanes, grains, "symetryPlaneFW")
symetryPlaneBW = createGroup(gobj, bwplanes, grains, "symetryPlaneBW")
symetryPlaneL = createGroup(gobj, lplanes, grains, "symetryPlaneL")
symetryPlaneR = createGroup(gobj, rplanes, grains, "symetryPlaneR")
# wall
allgroup = geompy.CreateGroup(gobj, geompy.ShapeType["FACE"])
faces = geompy.SubShapeAllIDs(gobj, geompy.ShapeType["FACE"])
geompy.UnionIDs(allgroup, faces)
wall = geompy.CutListOfGroups([allgroup],
[inlet, outlet, symetryPlaneFW, symetryPlaneBW, symetryPlaneL, symetryPlaneR], "wall")
boundary = {
"inlet": inlet,
"outlet": outlet,
"symetryPlaneFW": symetryPlaneFW,
"symetryPlaneBW": symetryPlaneBW,
"symetryPlaneL": symetryPlaneL,
"symetryPlaneR": symetryPlaneR,
"wall": wall
}
return boundary

View File

@ -1,110 +0,0 @@
import salome, GEOM, SMESH, SALOMEDS
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math
import os, sys
import logging
import time
from datetime import timedelta
def meshCreate(gobj, boundary, fineness, viscousLayers=None):
"""
Creates a mesh from a geometry.
Parameters:
fineness (int): Fineness of mesh.
0 - Very coarse,
1 - Coarse,
2 - Moderate,
3 - Fine,
4 - Very fine.
viscousLayers (dict or None): Defines viscous layers for mesh.
By default, inlets and outlets specified without layers.
{
"thickness": float,
"number": int,
"stretch": float
}
Returns:
Configured instance of class <SMESH.SMESH_Mesh>, containig the parameters and boundary groups.
"""
smesh = smeshBuilder.New()
Fineness = {
0: "Very coarse",
1: "Coarse",
2: "Moderate",
3: "Fine",
4: "Very fine"
}[fineness]
logging.info("meshCreate: mesh fineness - {}".format(Fineness))
mesh = smesh.Mesh(gobj)
netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
param = netgen.Parameters()
param.SetSecondOrder( 0 )
param.SetOptimize( 1 )
param.SetChordalError( -1 )
param.SetChordalErrorEnabled( 0 )
param.SetUseSurfaceCurvature( 1 )
param.SetFuseEdges( 1 )
param.SetCheckChartBoundary( 0 )
param.SetMinSize( 0.01 )
param.SetMaxSize( 0.1 )
param.SetFineness(fineness)
#param.SetGrowthRate( 0.1 )
#param.SetNbSegPerEdge( 5 )
#param.SetNbSegPerRadius( 10 )
param.SetQuadAllowed( 0 )
if not viscousLayers is None:
logging.info("meshCreate: viscous layers params - thickness = {}, number = {}, stretch factor = {}".format(
viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"]))
vlayer = netgen.ViscousLayers(
viscousLayers["thickness"],
viscousLayers["number"],
viscousLayers["stretch"],
[boundary["inlet"], boundary["outlet"]],
1, smeshBuilder.NODE_OFFSET)
else:
logging.info("meshCreate: viscous layers are disabled")
for name, b in boundary.items():
mesh.GroupOnGeom(b, "{}_".format(name), SMESH.FACE)
return mesh
def meshCompute(mobj):
"""Compute the mesh."""
status = mobj.Compute()
if status:
logging.info("Mesh succesfully computed.")
else:
logging.warning("Mesh is not computed.")
def meshExport(mobj, path):
"""
Export the mesh in a file in UNV format.
Parameters:
path (string): full path to the expected directory.
"""
exportpath = os.path.join(path, "{}.unv".format(mobj.name))
try:
mobj.ExportUNV(exportpath)
except:
logging.error("Cannot export mesh to '{}'".format(exportpath))

View File

@ -1,190 +0,0 @@
# state file generated using paraview version 5.8.1
# ----------------------------------------------------------------
# setup views used in the visualization
# ----------------------------------------------------------------
# trace generated using paraview version 5.8.1
#
# To ensure correct image size when batch processing, please search
# for and uncomment the line `# renderView*.ViewSize = [*,*]`
#### import the simple module from the paraview
from paraview.simple import *
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# Create a new 'Render View'
renderView1 = CreateView('RenderView')
renderView1.ViewSize = [1518, 833]
renderView1.AxesGrid = 'GridAxes3DActor'
renderView1.CenterOfRotation = [1.9999999835818016e-05, 2.0000001541120582e-05, 9.999999747378752e-06]
renderView1.StereoType = 'Crystal Eyes'
renderView1.CameraPosition = [1.9999999835818016e-05, 2.0000001541120582e-05, 0.0001558028029300748]
renderView1.CameraFocalPoint = [1.9999999835818016e-05, 2.0000001541120582e-05, 9.999999747378752e-06]
renderView1.CameraFocalDisk = 1.0
renderView1.CameraParallelScale = 3.773654229301616e-05
renderView1.CameraParallelProjection = 1
# init the 'GridAxes3DActor' selected for 'AxesGrid'
renderView1.AxesGrid.Visibility = 1
SetActiveView(None)
# ----------------------------------------------------------------
# setup view layouts
# ----------------------------------------------------------------
# create new layout object 'Layout #1'
layout1 = CreateLayout(name='Layout #1')
layout1.AssignView(0, renderView1)
# ----------------------------------------------------------------
# restore active view
SetActiveView(renderView1)
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# setup the data processing pipelines
# ----------------------------------------------------------------
# create a new 'OpenFOAMReader'
simplefoam = OpenFOAMReader(FileName='/home/nafaryus/projects/anisotrope-cube/build/simple/coeff-0.01/direction-100/simple.foam')
simplefoam.CaseType = 'Decomposed Case'
simplefoam.MeshRegions = ['internalMesh']
simplefoam.CellArrays = ['U', 'p']
# create a new 'Clip'
clip1 = Clip(Input=simplefoam)
clip1.ClipType = 'Plane'
clip1.HyperTreeGridClipper = 'Plane'
clip1.Scalars = ['POINTS', 'p']
clip1.Value = 0.000504692076901847
# init the 'Plane' selected for 'ClipType'
clip1.ClipType.Origin = [1.999999909685357e-05, 2.0000001029529813e-05, 9.999999747378752e-06]
clip1.ClipType.Normal = [0.0, 0.0, 1.0]
# init the 'Plane' selected for 'HyperTreeGridClipper'
clip1.HyperTreeGridClipper.Origin = [1.999999909685357e-05, 2.0000001029529813e-05, 9.999999747378752e-06]
# ----------------------------------------------------------------
# setup the visualization in view 'renderView1'
# ----------------------------------------------------------------
# show data from simplefoam
simplefoamDisplay = Show(simplefoam, renderView1, 'UnstructuredGridRepresentation')
# get color transfer function/color map for 'p'
pLUT = GetColorTransferFunction('p')
pLUT.RGBPoints = [-2.29675952141406e-05, 0.231373, 0.298039, 0.752941, 0.000504692076901847, 0.865003, 0.865003, 0.865003, 0.0010323517490178347, 0.705882, 0.0156863, 0.14902]
pLUT.ScalarRangeInitialized = 1.0
# get opacity transfer function/opacity map for 'p'
pPWF = GetOpacityTransferFunction('p')
pPWF.Points = [-2.29675952141406e-05, 0.0, 0.5, 0.0, 0.0010323517490178347, 1.0, 0.5, 0.0]
pPWF.ScalarRangeInitialized = 1
# trace defaults for the display properties.
simplefoamDisplay.Representation = 'Surface'
simplefoamDisplay.ColorArrayName = ['POINTS', 'p']
simplefoamDisplay.LookupTable = pLUT
simplefoamDisplay.OSPRayScaleArray = 'p'
simplefoamDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
simplefoamDisplay.SelectOrientationVectors = 'U'
simplefoamDisplay.ScaleFactor = 3.740525119155791e-06
simplefoamDisplay.SelectScaleArray = 'p'
simplefoamDisplay.GlyphType = 'Arrow'
simplefoamDisplay.GlyphTableIndexArray = 'p'
simplefoamDisplay.GaussianRadius = 1.8702625595778954e-07
simplefoamDisplay.SetScaleArray = ['POINTS', 'p']
simplefoamDisplay.ScaleTransferFunction = 'PiecewiseFunction'
simplefoamDisplay.OpacityArray = ['POINTS', 'p']
simplefoamDisplay.OpacityTransferFunction = 'PiecewiseFunction'
simplefoamDisplay.DataAxesGrid = 'GridAxesRepresentation'
simplefoamDisplay.PolarAxes = 'PolarAxesRepresentation'
simplefoamDisplay.ScalarOpacityFunction = pPWF
simplefoamDisplay.ScalarOpacityUnitDistance = 9.756969141785036e-07
simplefoamDisplay.ExtractedBlockIndex = 1
# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
simplefoamDisplay.ScaleTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 0.0010000000474974513, 1.0, 0.5, 0.0]
# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
simplefoamDisplay.OpacityTransferFunction.Points = [0.0, 0.0, 0.5, 0.0, 0.0010000000474974513, 1.0, 0.5, 0.0]
# show data from clip1
clip1Display = Show(clip1, renderView1, 'UnstructuredGridRepresentation')
# get color transfer function/color map for 'U'
uLUT = GetColorTransferFunction('U')
uLUT.RGBPoints = [0.0, 0.231373, 0.298039, 0.752941, 0.00010124565929526629, 0.865003, 0.865003, 0.865003, 0.00020249131859053257, 0.705882, 0.0156863, 0.14902]
uLUT.ScalarRangeInitialized = 1.0
# get opacity transfer function/opacity map for 'U'
uPWF = GetOpacityTransferFunction('U')
uPWF.Points = [0.0, 0.0, 0.5, 0.0, 0.00020249131859053257, 1.0, 0.5, 0.0]
uPWF.ScalarRangeInitialized = 1
# trace defaults for the display properties.
clip1Display.Representation = 'Surface'
clip1Display.ColorArrayName = ['POINTS', 'U']
clip1Display.LookupTable = uLUT
clip1Display.OSPRayScaleArray = 'p'
clip1Display.OSPRayScaleFunction = 'PiecewiseFunction'
clip1Display.SelectOrientationVectors = 'U'
clip1Display.ScaleFactor = 3.740524778095278e-06
clip1Display.SelectScaleArray = 'p'
clip1Display.GlyphType = 'Arrow'
clip1Display.GlyphTableIndexArray = 'p'
clip1Display.GaussianRadius = 1.870262389047639e-07
clip1Display.SetScaleArray = ['POINTS', 'p']
clip1Display.ScaleTransferFunction = 'PiecewiseFunction'
clip1Display.OpacityArray = ['POINTS', 'p']
clip1Display.OpacityTransferFunction = 'PiecewiseFunction'
clip1Display.DataAxesGrid = 'GridAxesRepresentation'
clip1Display.PolarAxes = 'PolarAxesRepresentation'
clip1Display.ScalarOpacityFunction = uPWF
clip1Display.ScalarOpacityUnitDistance = 1.1332065866368908e-06
clip1Display.ExtractedBlockIndex = 1
# init the 'PiecewiseFunction' selected for 'ScaleTransferFunction'
clip1Display.ScaleTransferFunction.Points = [-2.29675952141406e-05, 0.0, 0.5, 0.0, 0.0010323517490178347, 1.0, 0.5, 0.0]
# init the 'PiecewiseFunction' selected for 'OpacityTransferFunction'
clip1Display.OpacityTransferFunction.Points = [-2.29675952141406e-05, 0.0, 0.5, 0.0, 0.0010323517490178347, 1.0, 0.5, 0.0]
# setup the color legend parameters for each legend in this view
# get color legend/bar for pLUT in view renderView1
pLUTColorBar = GetScalarBar(pLUT, renderView1)
pLUTColorBar.Title = 'p'
pLUTColorBar.ComponentTitle = ''
# set color bar visibility
pLUTColorBar.Visibility = 0
# get color legend/bar for uLUT in view renderView1
uLUTColorBar = GetScalarBar(uLUT, renderView1)
uLUTColorBar.Title = 'U'
uLUTColorBar.ComponentTitle = 'Magnitude'
# set color bar visibility
uLUTColorBar.Visibility = 1
# hide data in view
Hide(simplefoam, renderView1)
# show color legend
clip1Display.SetScalarBarVisibility(renderView1, True)
# ----------------------------------------------------------------
# setup color maps and opacity mapes used in the visualization
# note: the Get..() functions create a new object, if needed
# ----------------------------------------------------------------
# ----------------------------------------------------------------
# finally, restore active source
SetActiveSource(clip1Display)
# ----------------------------------------------------------------
SaveScreenshot("test.png", renderView1)

View File

@ -1,287 +0,0 @@
import os, sys
from collections import namedtuple
import time
import logging
from datetime import timedelta
import multiprocessing
ROOT = os.getcwd()
sys.path.append(ROOT)
LOG = os.path.join(ROOT, "logs")
BUILD = os.path.join(ROOT, "build")
from src import utils
from src import salome_utils
if __name__ == "__main__":
if not os.path.exists(LOG):
os.makedirs(LOG)
if not os.path.exists(BUILD):
os.makedirs(BUILD)
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/cubic.log".format(LOG))
])
start_time = time.monotonic()
logging.info("Started at {}".format(timedelta(seconds=start_time)))
#nprocs = os.cpu_count()
#port = 2810
#salomeServers = []
#salomeServer = namedtuple("salomeServer", ["process", "port"])
#for n in range(nprocs):
# while not utils.portIsFree:
# port += 1
# p = multiprocessing.Process(target = salome_utils.startServer, args = (port,))
#s = salomeServer(salome_utils.startServer(port), port)
# s = salomeServer(p, port)
# salomeServers.append(s)
# port += 1
#var = []
# TODO: pass it to samples in namedtuples
# get sample.directions and etc ..
structures = [
"simpleCubic",
"bodyCenteredCubic",
"faceCenteredCubic"
]
directions = [
[1, 0, 0],
[0, 0, 1],
#[1, 1, 1]
]
#cmd = """python -c
#\"import sys;
#sys.append('{}');
#import samples;
#samples.genMesh('{}', {}, {}, '{}');
#\"
#""".replace("\n", "")
Task = namedtuple("Task", ["structure", "coeff", "direction", "saveto"])
tasks = []
for structure in structures:
if structure == "simpleCubic":
theta = [c * 0.01 for c in range(1, 14)]
elif structure == "faceCenteredCubic":
theta = []
elif structure == "bodyCenteredCubic":
theta = []
for coeff in theta:
for direction in directions:
saveto = os.path.join(BUILD, structure, "coeff-{}".format(coeff),
"direction-{}{}{}".format(direction[0], direction[1], direction[2]))
if not os.path.exists(saveto):
os.makedirs(saveto)
direction_ = "".join([str(coord) for coord in direction])
t = Task(structure, coeff, direction_, os.path.join(saveto, "mesh.unv"))
tasks.append(t)
#tasks.append(cmd.format(ROOT, structure, coeff, direction, BUILD))
logging.info("tasks: {}".format(len(tasks)))
#n = 0
#for cmd in tasks:
# var.append((salomeServer[n].port, cmd))
# n = n + 1 if n < 3 else 0
#for s in salomeServers:
# s.process.start()
scriptpath = os.path.join(ROOT, "samples/__init__.py")
port = 2810
for task in tasks:
returncode = salome_utils.runExecute(port, scriptpath, task.structure, task.coeff, task.direction, task.saveto)
logging.info("Return code: {}".format(returncode))
if returncode == 1:
break
#res = utils.parallel(nprocs, var, salome_utils.remote)
#print(res)
#for s in salomeServers:
# s.process.join()
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
"""
if __name__ == "__main__":
###
# SALOME
#
# Get main paths
project = os.getcwd()
src = os.path.join(project, "src")
build = os.path.join(project, "build")
if not os.path.exists(build):
os.makedirs(build)
# Logger
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/genmesh.log".format(build))
])
start_time = time.monotonic()
# Start in parallel
processes = []
structures = ["simpleCubic", "faceCenteredCubic", "bodyCenteredCubic"]
directions = ["001"] #, "100", "111"]
coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ]
port = 2810
for structure in structures:
for direction in directions:
for coefficient in coefficients:
src_path = os.path.join(src, "{}.py".format(structure))
build_path = os.path.join(build,
structure,
"direction-{}".format(direction),
"alpha-{}".format(coefficient))
if not os.path.exists(build_path):
os.makedirs(build_path)
p = multiprocessing.Process(target = salome,
args = (port, src_path, build_path, coefficient, direction))
processes.append(p)
p.start()
logging.info("{} on port {}.".format(p, port))
port += 1
for process in processes:
process.join()
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
###
# FOAM
#
# Get main paths
project = os.getcwd()
src = os.path.join(project, "src")
build = os.path.join(project, "build")
if not os.path.exists(build):
os.makedirs(build)
# Logger
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/foam.log".format(build))
])
start_time = time.monotonic()
# Main entry
structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"]
directions = ["001"] #, "100", "111"]
coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ]
for structure in structures:
for direction in directions:
for coefficient in coefficients:
foamCase = [ "0", "constant", "system" ]
src_path = os.path.join(src, "{}Foam".format(structure))
build_path = os.path.join(build,
structure,
"direction-{}".format(direction),
"alpha-{}".format(coefficient))
logging.info("Entry with parameters: {}, direction = {}, alpha = {}".format(structure, direction, coefficient))
logging.info("Copying baseFOAM case ...")
for d in foamCase:
if not os.path.exists(os.path.join(build_path, d)):
shutil.copytree(os.path.join(src_path, d),
os.path.join(build_path, d))
os.chdir(build_path)
case_path = "."
logging.info("Importing mesh to foam ...")
ideasUnvToFoam(case_path, "{}-{}-{}.unv".format(structure, direction, coefficient))
logging.info("Creating patches ...")
createPatch(case_path)
logging.info("Scaling mesh ...")
transformPoints(case_path, "(1e-5 1e-5 1e-5)")
logging.info("Checking mesh ...")
checkMesh(case_path)
#logging.info("Changing mesh boundaries types ...")
#foamDictionarySet(case_path, "constant/polyMesh/boundary", "entry0.wall.type", "wall")
#foamDictionarySet(case_path, "constant/polyMesh/boundary", "entry0.symetryPlane.type", "symetryPlane")
logging.info("Decomposing case ...")
decomposePar(case_path)
logging.info("Evaluating initial approximation via potentialFoam ...")
potentialFoam(case_path)
logging.info("Preparing boundaryFields for simpleFoam ...")
for n in range(4):
foamDictionarySet(case_path, "processor{}/0/U".format(n),
"boundaryField.inlet.type", "pressureInletVelocity")
foamDictionarySet(case_path, "processor{}/0/U",
"boundaryField.inlet.value", "uniform (0 0 0)")
logging.info("Calculating ...")
simpleFoam(case_path)
os.chdir(project)
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
"""

View File

@ -1,27 +0,0 @@
#!/usr/bin/env bash
ideasUnvToFoam mesh.unv
createPatch -overwrite
#patchSummary
transformPoints -scale '1e-5'
foamDictionary constant/polyMesh/boundary -entry entry0.cyclicPlaneL.separationVector -set '(0 1e-5 0)'
foamDictionary constant/polyMesh/boundary -entry entry0.cyclicPlaneR.separationVector -set '(0 -1e-5 0)'
checkMesh -allGeometry -allTopology | tee -a checkMesh.log
#
decomposePar
mpirun -np 4 --oversubscribe potentialFoam -parallel | tee -a potentialFoam.log
foamDictionary processor0/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor1/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor2/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor3/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor0/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor1/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor2/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor3/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
mpirun -np 4 --oversubscribe simpleFoam -parallel | tee -a simpleFoam.log

View File

@ -1,27 +0,0 @@
#!/usr/bin/env bash
ideasUnvToFoam mesh.unv
createPatch -overwrite
#patchSummary
transformPoints -scale '1e-5'
foamDictionary constant/polyMesh/boundary -entry entry0.cyclicPlaneL.separationVector -set '(0 1e-5 0)'
foamDictionary constant/polyMesh/boundary -entry entry0.cyclicPlaneR.separationVector -set '(0 -1e-5 0)'
checkMesh -allGeometry -allTopology | tee -a checkMesh.log
#
decomposePar
mpirun -np 4 --oversubscribe potentialFoam -parallel | tee -a potentialFoam.log
foamDictionary processor0/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor1/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor2/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor3/0/U -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary processor0/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor1/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor2/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
foamDictionary processor3/0/U -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
mpirun -np 4 --oversubscribe simpleFoam -parallel | tee -a simpleFoam.log

View File

@ -1,100 +0,0 @@
from collections import namedtuple
import os, sys
import logging
from pyquaternion import Quaternion
import math
ROOT = "/home/nafaryus/projects/anisotrope-cube"
sys.path.append(ROOT)
LOG = os.path.join(ROOT, "logs")
import salome
from simpleCubic import simpleCubic
from faceCenteredCubic import faceCenteredCubic
from bodyCenteredCubic import bodyCenteredCubic
from src import geometry_utils
from src import mesh_utils
def genMesh(stype, theta, fillet, flowdirection, saveto):
_G = globals()
structure = _G.get(stype)
if structure:
salome.salome_init()
grains, geometry1, geometry2 = structure(theta, fillet)
geometry = geometry1
if flowdirection == [1, 1, 1]:
geometry = geometry2
norm = [-1, 1, 0]
bcount = 6
# initial angle
angle = math.pi / 6
v1 = Quaternion(axis = norm, angle = math.pi / 2).rotate(flowdirection)
normvec = Quaternion(axis = flowdirection, angle = angle).rotate(v1)
direction = [1, 1, 1]
if flowdirection == [1, 0, 0]:
normvec = [0, 0, 1]
bcount = 4
direction = [1, 1, 0]
if flowdirection == [0, 0, 1]:
normvec = [1, 1, 0]
bcount = 4
direction = [0, 0, 1]
#
#geometry = geometry_utils.geompy.RemoveExtraEdges(geometry, False)
#
#boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains)
fineness = 1
viscousLayers = {
"thickness": 0.0001,
"number": 3,
"stretch": 1.2
}
mesh = mesh_utils.meshCreate(geometry, boundary, fineness, viscousLayers)
mesh_utils.meshCompute(mesh)
mesh_utils.meshExport(mesh, saveto)
salome.salome_close()
else:
raise Exception("Unknown sample function")
if __name__ == "__main__":
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/cubic.log".format(LOG))
])
stype = str(sys.argv[1])
theta = float(sys.argv[2])
fillet = True if int(sys.argv[3]) == 1 else False
flowdirection = [int(coord) for coord in sys.argv[4]]
saveto = str(sys.argv[5])
logging.info("""genMesh:
structure type:\t{}
coefficient:\t{}
fillet:\t{}
flow direction:\t{}
export path:\t{}""".format(stype, theta, fillet, flowdirection, saveto))
genMesh(stype, theta, fillet, flowdirection, saveto)

View File

@ -1,186 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def bodyCenteredCubic(alpha, fillet = False):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
V_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
V_2 = geompy.MakeVectorDXDYDZ(1, -1, 0)
L = 1.0
r0 = L*math.sqrt(3)/4
d1 = L*math.sqrt(2)
Pi_4 = 45*math.pi/180.0
n = 3
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, 2*L) # Vertex_1 of Rombus
sk.addPointsAbsolute(L, 0, L) # Vertex_2 of Rombus
sk.addPointsAbsolute(L, L, 0) # Vertex_3 of Rombus
sk.addPointsAbsolute(0, L, L) # Vertex_4 of Rombus
sk.addPointsAbsolute(0, 0, 2*L) # Vertex_1 of Rombus
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Extrusion_1 = geompy.MakePrismDXDYDZ(Face_1, L/2, L/2, L/2)
sk2 = geompy.Sketcher3D()
sk2.addPointsAbsolute(L/3, L/3, 4*L/3) # Vertex_1 of Triangle
sk2.addPointsAbsolute(L, 0, L) # Vertex_2 of Triangle
sk2.addPointsAbsolute(2*L/3, 2*L/3, 2*L/3) # Vertex_3 of Triangle
sk2.addPointsAbsolute(L/3, L/3, 4*L/3) # Vertex_1 of Triangle
a3D_Sketcher_2 = sk2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Extrusion_2 = geompy.MakePrismDXDYDZ(Face_2, L/2, L/2, L/2)
sk3 = geompy.Sketcher3D()
sk3.addPointsAbsolute(L/3, L/3, 4*L/3) # Vertex_1 of Hexagon
sk3.addPointsAbsolute(L, 0, L) # Vertex_2 of Hexagon
sk3.addPointsAbsolute(4*L/3, L/3, L/3) # Vertex_3 of Hexagon
#sk3.addPointsAbsolute(2*L/3, 2*L/3, 2*L/3) # Vertex_3 of Hexagon
sk3.addPointsAbsolute(L, L, 0) # Vertex_4 of Hexagon
sk3.addPointsAbsolute(L/3, 4*L/3, L/3) # Vertex_5 of Hexagon
sk3.addPointsAbsolute(0, L, L) # Vertex_6 of Hexagon
sk3.addPointsAbsolute(L/3, L/3, 4*L/3) # Vertex_1 of Hexagon
a3D_Sketcher_3 = sk3.wire()
Face_3 = geompy.MakeFaceWires([a3D_Sketcher_3], 1)
Extrusion_3 = geompy.MakePrismDXDYDZ(Face_3, L/2, L/2, L/2)
a3D_Sketcher_4 = geompy.MakeTranslation(a3D_Sketcher_3, -L/4, -L/4, -L/4)
Face_4 = geompy.MakeFaceWires([a3D_Sketcher_4], 1)
#Extrusion_4 = geompy.MakePrismDXDYDZ(Face_4, (n-1.5)*L, (n-1.5)*L, (n-1.5)*L) # Extrusion_3Direction
Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 1)
#Extrusion_4 = geompy.MakePrismVecH(Face_4, Vector_1, 4*(n-1.5)*r0) # Extrusion_3Direction
Extrusion_4 = geompy.MakePrismVecH(Face_4, Vector_1, 4*(n-2.0)*r0) # Extrusion_3Direction
Box_1 = geompy.MakeBoxDXDYDZ(d1, d1, L)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, Pi_4)
Translation_1 = geompy.MakeTranslation(Rotation_1, L, 0, 0)
Translation_2 = geompy.MakeTranslation(Translation_1, 0, 0, L/4)
Box_2 = geompy.MakeBoxDXDYDZ(d1/2, d1/2, L)
Rotation_2 = geompy.MakeRotation(Box_2, OZ, Pi_4)
Translation_3 = geompy.MakeTranslation(Rotation_2, L, 0, L/4)
Box_3 = geompy.MakeBoxDXDYDZ(d1/2, d1/2, L/2)
Rotation_3 = geompy.MakeRotation(Box_3, OZ, Pi_4)
Translation_4 = geompy.MakeTranslation(Rotation_3, L, 0, L/4)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( V_1, 'V_1' )
geompy.addToStudy( V_2, 'V_2' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Translation_1, 'Translation_1' )
geompy.addToStudy( Translation_2, 'Translation_2' )
geompy.addToStudy( Box_2, 'Box_2' )
geompy.addToStudy( Rotation_2, 'Rotation_2' )
geompy.addToStudy( Translation_3, 'Translation_3' )
geompy.addToStudy( Rotation_3, 'Rotation_3' )
geompy.addToStudy( Translation_4, 'Translation_4' )
#geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
#geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
#geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
geompy.addToStudy( a3D_Sketcher_3, 'a3D_Sketcher_3' )
#geompy.addToStudy( Face_3, 'Face_3' )
#geompy.addToStudy( Extrusion_3, 'Extrusion_3' )
geompy.addToStudy( a3D_Sketcher_4, 'a3D_Sketcher_4' )
#geompy.addToStudy( Face_4, 'Face_4' )
geompy.addToStudy( Extrusion_4, 'Extrusion_4' )
Radius = r0/(1-alpha)
Sphere_1 = geompy.MakeSphereR(Radius)
Down = geompy.MakeMultiTranslation2D(Sphere_1, OX, L, n, OY, L, n)
Up = geompy.MakeTranslation(Down, 0, 0, (n-1)*L)
Up_Down = geompy.MakeMultiTranslation1D(Down, OZ, L, n)
Cut_1a = geompy.MakeCutList(Translation_1, [Up_Down])
Cut_1b = geompy.MakeCutList(Translation_2, [Up_Down])
Cut_2a = geompy.MakeCutList(Translation_3, [Up_Down])
Cut_2b = geompy.MakeCutList(Translation_4, [Up_Down])
Cut_3a = geompy.MakeCutList(Extrusion_1, [Up_Down])
Cut_3b = geompy.MakeCutList(Extrusion_2, [Up_Down])
Cut_3c = geompy.MakeCutList(Extrusion_4, [Up_Down])
Middle = geompy.MakeTranslation(Up_Down, L/2, L/2, L/2)
Pore_1a = geompy.MakeCutList(Cut_1a, [Middle])
Pore_1b = geompy.MakeCutList(Cut_1b, [Middle])
Pore_2a = geompy.MakeCutList(Cut_2a, [Middle])
Pore_2b = geompy.MakeCutList(Cut_2b, [Middle])
Pore_3a = geompy.MakeCutList(Cut_3a, [Middle])
Pore_3b = geompy.MakeCutList(Cut_3b, [Middle])
Pore_3c = geompy.MakeCutList(Cut_3c, [Middle])
geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Down, 'Down_' )
geompy.addToStudy( Up, 'Up_' )
geompy.addToStudy( Up_Down, 'Up_Down_' )
geompy.addToStudy( Cut_1a, 'Cut_1a_' )
geompy.addToStudy( Cut_1b, 'Cut_1b_' )
geompy.addToStudy( Middle, 'Middle_' )
geompy.addToStudy( Pore_1a, 'Pore_1a_' )
geompy.addToStudy( Pore_1b, 'Pore_1b_' )
geompy.addToStudy( Pore_2a, 'Pore_2a_' )
geompy.addToStudy( Pore_2b, 'Pore_2b_' )
geompy.addToStudy( Pore_3a, 'Pore_3a_' )
geompy.addToStudy( Pore_3b, 'Pore_3b_' )
geompy.addToStudy( Pore_3c, 'Pore_3c_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(Middle, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
geometry1 = Pore_1a
geometry2 = Pore_3c
if fillet:
R = r0 / (1 - alpha)
C1 = 0.8
C2 = 0.05
alpha1 = 0.01
alpha2 = 0.18
Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
R_fillet = Cf * (R - r0)
# Scaling up
scale = 100
grains = geompy.MakeScaleTransform(grains, O, scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, scale)
#
grains = geompy.MakeFilletAll(grains, R_fillet * scale)
geometry1 = geompy.MakeCutList(geometry1, [grains], True)
geometry2 = geompy.MakeCutList(geometry2, [grains], True)
# Scaling down
grains = geompy.MakeScaleTransform(grains, O, 1 / scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, 1 / scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, 1 / scale)
#
geompy.addToStudy(grains, "grains")
geompy.addToStudy(geometry1, "geometry1")
geompy.addToStudy(geometry2, "geometry2")
return grains, geometry1, geometry2

View File

@ -1,181 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def faceCenteredCubic(alpha, fillet = False):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
V = geompy.MakeVectorDXDYDZ(1, 1, 1)
V_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
V_2 = geompy.MakeVectorDXDYDZ(1, -1, 0)
L = 1.0
r0 = L*math.sqrt(2)/4
a=2*r0
b=L/2
#c=b/2
d = L*math.sqrt(3)
h = d/3
Pi_4 = 45*math.pi/180.0
n = 3 # number of cubes in line
Vertex_1 = geompy.MakeVertex(-a, -a, -b)
Vertex_2 = geompy.MakeVertex(a, a, b)
Vertex_3 = geompy.MakeVertex(-b*(n-1), 0, -b*(n-2)) # Center of Sphere_1
Vertex_4 = geompy.MakeVertex(-b*n, 0, -b*(n-3)) # Center of Sphere_2
Vertex_5 = geompy.MakeVertex(0, 0, -b)
sk = geompy.Sketcher3D() # Rombus
sk.addPointsAbsolute(-b, -b, b) # Vertex_1 of Rombus
sk.addPointsAbsolute(0, -b, 0) # Vertex_2 of Rombus
sk.addPointsAbsolute(0, 0, -b) # Vertex_3 of Rombus
sk.addPointsAbsolute(-b, 0, 0) # Vertex_4 of Rombus
sk.addPointsAbsolute(-b, -b, b) # Vertex_1 of Rombus
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Face_1_Up = geompy.MakeTranslation(Face_1, b, b, 0)
Rhombohedron = geompy.MakeHexa2Faces(Face_1, Face_1_Up)
sk_2 = geompy.Sketcher3D() # Triangle
sk_2.addPointsAbsolute(-2*b/3, -2*b/3, b/3) # Vertex_1 of Triangle
sk_2.addPointsAbsolute(0, -b, 0) # Vertex_2 of Triangle
sk_2.addPointsAbsolute(-b/3, -b/3, -b/3) # Vertex_3 of Triangle
sk_2.addPointsAbsolute(-2*b/3, -2*b/3, b/3) # Vertex_1 of Triangle
a3D_Sketcher_2 = sk_2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Extrusion_2 = geompy.MakePrismVecH(Face_2, V, h)
sk_3 = geompy.Sketcher3D() # Hexagon
sk_3.addPointsAbsolute(-2*b/3, -2*b/3, b/3) # Vertex_1 of Hexagon
sk_3.addPointsAbsolute(0, -b, 0) # Vertex_2 of Hexagon
sk_3.addPointsAbsolute(b/3, -2*b/3, -2*b/3) # Vertex_3 of Hexagon
sk_3.addPointsAbsolute(0, 0, -b) # Vertex_4 of Hexagon
sk_3.addPointsAbsolute(-2*b/3, b/3, -2*b/3) # Vertex_5 of Hexagon
sk_3.addPointsAbsolute(-b, 0, 0) # Vertex_6 of Hexagon
sk_3.addPointsAbsolute(-2*b/3, -2*b/3, b/3) # Vertex_1 of Hexagon
a3D_Sketcher_3 = sk_3.wire()
Face_3 = geompy.MakeFaceWires([a3D_Sketcher_3], 1)
Extrusion_3 = geompy.MakePrismVecH(Face_3, V, h)
Translation_1 = geompy.MakeTranslation(a3D_Sketcher_3, -(n-2)*L/3, -(n-2)*L/3, -(n-2)*L/3)
Face_4 = geompy.MakeFaceWires([Translation_1], 1)
Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 1)
Extrusion_4 = geompy.MakePrismVecH(Face_4, Vector_1, (2*n-3)/3.0*d) # Extrusion_3Direction
Box_1 = geompy.MakeBoxTwoPnt(Vertex_1, Vertex_2)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, Pi_4)
Box_2 = geompy.MakeBoxTwoPnt(Vertex_5, Vertex_2)
Rotation_2 = geompy.MakeRotation(Box_2, OZ, Pi_4)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( V_1, 'V_1' )
geompy.addToStudy( V_2, 'V_2' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Box_2, 'Box_2' )
geompy.addToStudy( Rotation_2, 'Rotation_2_' )
geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Face_1_Up, 'Face_1_Up' )
geompy.addToStudy( Rhombohedron, 'Rhombohedron' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
geompy.addToStudy( a3D_Sketcher_3, 'a3D_Sketcher_3' )
geompy.addToStudy( Face_3, 'Face_3' )
geompy.addToStudy( Extrusion_3, 'Extrusion_3' )
geompy.addToStudy( Face_4, 'Face_4' )
geompy.addToStudy( Extrusion_4, 'Extrusion_4' )
Radius = r0/(1-alpha)
Sphere_1 = geompy.MakeSpherePntR(Vertex_3, Radius)
Down = geompy.MakeMultiTranslation2D(Sphere_1, V_1, a, n, V_2, a, n)
Up_Down = geompy.MakeMultiTranslation1D(Down, OZ, 1, n-1)
Cut_1 = geompy.MakeCutList(Rotation_1, [Up_Down], True)
Sphere_2 = geompy.MakeSpherePntR(Vertex_4, Radius)
Middle = geompy.MakeMultiTranslation2D(Sphere_2, V_1, a, n+1, V_2, a, n+1)
Middle_Up = geompy.MakeMultiTranslation1D(Middle, OZ, 1, n-2)
Cut_2 = geompy.MakeCutList(Cut_1, [Middle_Up], True)
Common = geompy.MakeCommonList([Cut_2, Rotation_2], True)
Pore_3 = geompy.MakeCommonList([Rhombohedron, Cut_2], True)
Cut_3 = geompy.MakeCutList(Extrusion_2, [Up_Down], True)
Cut_4 = geompy.MakeCutList(Cut_3, [Middle_Up], True)
Cut_5 = geompy.MakeCutList(Extrusion_3, [Up_Down], True)
Cut_6 = geompy.MakeCutList(Cut_5, [Middle_Up], True)
Cut_7 = geompy.MakeCutList(Extrusion_4, [Up_Down], True)
Cut_8 = geompy.MakeCutList(Cut_7, [Middle_Up], True)
#geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Down, 'Down_' )
geompy.addToStudy( Up_Down, 'Up_Down_' )
geompy.addToStudy( Cut_1, 'Cut_1_' )
geompy.addToStudy( Middle_Up, 'Middle_Up_' )
geompy.addToStudy( Cut_2, 'Pore_' )
geompy.addToStudy( Common, 'Pore_2_' )
geompy.addToStudy( Pore_3, 'Pore_3_' )
geompy.addToStudy( Cut_4, 'Pore_3a_' )
geompy.addToStudy( Cut_6, 'Pore_3b_' )
geompy.addToStudy( Cut_8, 'Pore_3c_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(Middle_Up, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
geometry1 = Common
geometry2 = Cut_8
if fillet:
R = r0 / (1 - alpha)
C1 = 0.8
C2 = 0.05
alpha1 = 0.01
alpha2 = 0.13
Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
R_fillet = Cf * (R - r0)
# Scaling up
scale = 100
grains = geompy.MakeScaleTransform(grains, O, scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, scale)
#
grains = geompy.MakeFilletAll(grains, R_fillet * scale)
geometry1 = geompy.MakeCutList(geometry1, [grains], True)
geometry2 = geompy.MakeCutList(geometry2, [grains], True)
# Scaling down
grains = geompy.MakeScaleTransform(grains, O, 1 / scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, 1 / scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, 1 / scale)
#
geompy.addToStudy(grains, "grains")
geompy.addToStudy(geometry1, "geometry1")
geompy.addToStudy(geometry2, "geometry2")
return grains, geometry1, geometry2

View File

@ -1,234 +0,0 @@
import os, sys
from collections import namedtuple
import time
import logging
from datetime import timedelta
import multiprocessing
import shutil
ROOT = os.getcwd()
sys.path.append(ROOT)
LOG = os.path.join(ROOT, "logs")
BUILD = os.path.join(ROOT, "build")
from src import utils
from src import salome_utils
from src import foam_utils
def createTasks():
###
# Control variables
##
structures = [
"simple",
#"bodyCentered",
#"faceCentered"
]
directions = [
#[1, 0, 0],
[0, 0, 1],
#[1, 1, 1]
]
fillet = 1
###
# Tasks
##
Task = namedtuple("Task", ["structure", "coeff", "fillet", "direction", "saveto"])
tasks = []
for structure in structures:
if structure == "simple":
theta = [0.01] #[c * 0.01 for c in range(1, 28 + 1)]
#theta = [ 0.01, 0.28 ]
elif structure == "faceCentered":
#theta = [c * 0.01 for c in range(1, 13 + 1)]
theta = [ 0.01, 0.13 ]
elif structure == "bodyCentered":
#theta = [c * 0.01 for c in range(1, 18 + 1)]
theta = [ 0.01, 0.13, 0.14, 0.18 ]
for coeff in theta:
for direction in directions:
saveto = os.path.join(BUILD, structure, "coeff-{}".format(coeff),
"direction-{}{}{}".format(direction[0], direction[1], direction[2]))
if not os.path.exists(saveto):
os.makedirs(saveto)
t = Task(structure, coeff, fillet, direction, saveto)
tasks.append(t)
return tasks
def createMesh(tasks):
scriptpath = os.path.join(ROOT, "samples/__init__.py")
port = 2810
errors = 0
for task in tasks:
logging.info("-" * 80)
logging.info("""createMesh:
task:\t{} / {}""".format(tasks.index(task) + 1, len(tasks)))
start_time = time.monotonic()
returncode = salome_utils.runExecute(port, scriptpath, task.structure, task.coeff, task.fillet, "".join([str(coord) for coord in task.direction]), os.path.join(task.saveto, "mesh.unv"))
end_time = time.monotonic()
logging.info("createMesh: elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
logging.info("salome: return code:\t{}".format(returncode))
if returncode == 1:
#break
errors += 1
pass
return errors
def calculate(tasks):
foamCase = [ "0", "constant", "system" ]
rmDirs = ["0", "constant", "system", "postProcessing", "logs"] + [ "processor{}".format(n) for n in range(4)]
for task in tasks:
for d in rmDirs:
if os.path.exists(os.path.join(task.saveto, d)):
shutil.rmtree(os.path.join(task.saveto, d))
for d in foamCase:
if not os.path.exists(os.path.join(task.saveto, d)):
shutil.copytree(os.path.join(ROOT, "src/cubicFoam", d),
os.path.join(task.saveto, d))
os.chdir(task.saveto)
casepath = "."
logging.info("-" * 80)
logging.info("""calculate:
task:\t{} / {}
structure type:\t{}
coefficient:\t{}
flow direction:\t{}
path:\t{}""".format(
tasks.index(task) + 1,
len(tasks) ,
task.structure,
task.coeff,
task.direction,
task.saveto
))
foam_utils.ideasUnvToFoam("mesh.unv")
foam_utils.createPatch(dictfile = "system/createPatchDict.symetry")
foam_utils.checkMesh()
scale = (1e-5, 1e-5, 1e-5)
foam_utils.transformPoints(scale)
foam_utils.decomposePar()
foam_utils.renumberMesh()
foam_utils.potentialFoam()
for n in range(4):
foam_utils.foamDictionary("processor{}/0/U".format(n),
"boundaryField.inlet.type", "pressureInletVelocity")
foam_utils.foamDictionary("processor{}/0/U".format(n),
"boundaryField.inlet.value", "uniform (0 0 0)")
foam_utils.simpleFoam()
os.chdir(ROOT)
def postprocessing(tasks):
surfaceFieldValue = {}
porosity = {}
for task in tasks:
direction = "direction-{}{}{}".format(task.direction[0], task.direction[1], task.direction[2])
path = os.path.join(BUILD, task.structure, "postProcessing", direction)
surfaceFieldValuePath = os.path.join(task.saveto, "postProcessing", "")
if not os.path.exists(path):
os.makedirs(path)
surfaceFieldValues = [ line.strip().split() for line in open(surfaceFieldValuePath, "r").readlines() ]
with open(os.path.join(path, "porosity.dat")) as io:
io.write("{}\t{}".format(task.coeff, surfaceFieldValues[-1][1]))
if __name__ == "__main__":
if not os.path.exists(LOG):
os.makedirs(LOG)
if not os.path.exists(BUILD):
os.makedirs(BUILD)
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler("{}/cubic.log".format(LOG))
])
# TODO: add force arg
Args = namedtuple("Args", ["mesh", "calc", "pp"])
if len(sys.argv) > 1:
action = sys.argv[1]
if action == "mesh":
args = Args(True, False, False)
elif action == "calc":
args = Args(False, True, False)
elif action == "pp":
args = Args(False, False, True)
elif action == "all":
args = Args(True, True, True)
else:
args = Args(True, True, True)
tasks = createTasks()
logging.info("Tasks: {}".format(len(tasks)))
if args.mesh:
start_time = time.monotonic()
#logging.info("Started at {}".format(timedelta(seconds=start_time)))
errors = createMesh(tasks)
end_time = time.monotonic()
logging.info("-" * 80)
logging.info("Elapsed time:\t{}".format(timedelta(seconds=end_time - start_time)))
logging.info("Errors:\t{}".format(errors))
if args.calc:
start_time = time.monotonic()
#logging.info("Started at {}".format(timedelta(seconds=start_time)))
calculate(tasks)
end_time = time.monotonic()
logging.info("-" * 80)
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
if args.pp:
postprocessing(tasks)

View File

@ -1,150 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def simpleCubic(alpha, fillet = False):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
r0 = 1.0
L = 2*r0
h = L
h2 = 2*h
d1 = L*math.sqrt(2)
d2 = L*math.sqrt(3)
d = r0/math.sqrt(3)
pi_4 = 45*math.pi/180.0
pi_2 = 90*math.pi/180.0
n = 3 # number of cubes in line
Box_1 = geompy.MakeBoxDXDYDZ(d1, d1, h)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, pi_4)
Translation_2 = geompy.MakeTranslation(Rotation_1, h, 0, 0)
Vertex_1 = geompy.MakeVertex(h, 0, 0)
Vertex_2 = geompy.MakeVertex(h, h, 0)
Vertex_3 = geompy.MakeVertex(h, h, h)
Line_1 = geompy.MakeLineTwoPnt(Vertex_2, Vertex_3)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, h2) # Rombus
sk.addPointsAbsolute(h, 0, h) # Vertex_2 of Rombus
sk.addPointsAbsolute(h, h, 0) # Vertex_3 of Rombus
sk.addPointsAbsolute(0, h, h) # Vertex_4 of Rombus
sk.addPointsAbsolute(0, 0, h2) # Vertex_1 of Rombus
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
Extrusion_1 = geompy.MakePrismVecH(Face_1, Vector_1, d1)
sk_2 = geompy.Sketcher3D() # Hexagon
sk_2.addPointsAbsolute(L, L, L) # Vertex_1 of Hexagon
sk_2.addPointsAbsolute(5*L/3, 2*L/3, 2*L/3) # Vertex_2 of Hexagon
sk_2.addPointsAbsolute(2*L, L, 0) # Vertex_3 of Hexagon
sk_2.addPointsAbsolute(5*L/3, 5*L/3, -L/3) # Vertex_4 of Hexagon
sk_2.addPointsAbsolute(L, 2*L, 0) # Vertex_5 of Hexagon
sk_2.addPointsAbsolute(2*L/3, 5*L/3, 2*L/3) # Vertex_6 of Hexagon
sk_2.addPointsAbsolute(L, L, L) # Vertex_1 of Hexagon
a3D_Sketcher_2 = sk_2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Vector_2 = geompy.MakeVectorDXDYDZ(1, 1, 1)
Extrusion_2 = geompy.MakePrismVecH(Face_2, Vector_2, d2/3)
Translation_3 = geompy.MakeTranslation(a3D_Sketcher_2, -L-L/6, -L-L/6, 0-L/6)
Face_3 = geompy.MakeFaceWires([Translation_3], 1)
Vector_2 = geompy.MakeVectorDXDYDZ(1, 1, 1)
#Extrusion_3 = geompy.MakePrismVecH(Face_3, Vector_2, (n-4.0/3)*d2) # Extrusion_3Direction
Extrusion_3 = geompy.MakePrismVecH(Face_3, Vector_2, (n-2.0)*d2) # Extrusion_3Direction
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( Vertex_1, 'Vertex_1' )
geompy.addToStudy( Vertex_2, 'Vertex_2' )
geompy.addToStudy( Vertex_3, 'Vertex_3' )
geompy.addToStudy( Line_1, 'Line_1' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Translation_2, 'Translation_2' )
geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Vector_1, 'Vector_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Vector_2, 'Vector_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
geompy.addToStudy( Translation_3, 'a3D_Sketcher_3' )
geompy.addToStudy( Extrusion_3, 'Extrusion_3' )
Sphere_1 = geompy.MakeSphereR(r0/(1-alpha))
Multi_Translation_2 = geompy.MakeMultiTranslation2D(Sphere_1, OX, L, n, OY, L, n)
Multi_Translation_3 = geompy.MakeMultiTranslation1D(Multi_Translation_2, OZ, L, n)
Cut_1 = geompy.MakeCutList(Translation_2, [Multi_Translation_3])
Cut_2 = geompy.MakeCutList(Extrusion_1, [Multi_Translation_3])
Cut_3 = geompy.MakeCutList(Extrusion_2, [Multi_Translation_3])
Cut_V = geompy.MakeCutList(Cut_1, [Cut_3])
#Cut_2.SetColor(SALOMEDS.Color(0,0,1))
Cut_4 = geompy.MakeCutList(Extrusion_3, [Multi_Translation_3])
geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Multi_Translation_2, 'Multi-Translation_2_' )
geompy.addToStudy( Multi_Translation_3, 'Multi-Translation_3_' )
geompy.addToStudy( Cut_1, 'Pore1_' )
geompy.addToStudy( Cut_2, 'Pore2_' )
geompy.addToStudy( Cut_3, 'Pore3_' )
geompy.addToStudy( Cut_V, 'Cut_V_' )
geompy.addToStudy( Cut_4, 'Pore4_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
grains = geompy.ExtractShapes(Multi_Translation_3, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
geometry1 = Cut_1
geometry2 = Cut_4
if fillet:
R = r0 / (1 - alpha)
C1 = 0.8
C2 = 0.05
alpha1 = 0.01
alpha2 = 0.28
Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
R_fillet = Cf * (R - r0)
# Scaling up
scale = 100
grains = geompy.MakeScaleTransform(grains, O, scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, scale)
#
grains = geompy.MakeFilletAll(grains, R_fillet * scale)
geometry1 = geompy.MakeCutList(geometry1, [grains], True)
geometry2 = geompy.MakeCutList(geometry2, [grains], True)
# Scaling down
grains = geompy.MakeScaleTransform(grains, O, 1 / scale)
geometry1 = geompy.MakeScaleTransform(geometry1, O, 1 / scale)
geometry2 = geompy.MakeScaleTransform(geometry2, O, 1 / scale)
#
geompy.addToStudy(grains, "grains")
geompy.addToStudy(geometry1, "geometry1")
geompy.addToStudy(geometry2, "geometry2")
return grains, geometry1, geometry2

View File

@ -1,127 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def bodyCenteredCubic(alpha):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
V_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
V_2 = geompy.MakeVectorDXDYDZ(1, -1, 0)
L = 1.0
L2 = L/2
r0 = L*math.sqrt(3)/4
a = L*math.sqrt(2)
a2 = a/2
h4 = L/4
Pi_4 = 45*math.pi/180.0
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, 2*L) #Vertex_1 of Rombus
sk.addPointsAbsolute(L, 0, L) #Vertex_2 of Rombus
sk.addPointsAbsolute(L, L, 0) #Vertex_3 of Rombus
sk.addPointsAbsolute(0, L, L) #Vertex_4 of Rombus
sk.addPointsAbsolute(0, 0, 2*L) #Vertex_1 of Rombus
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Extrusion_1 = geompy.MakePrismDXDYDZ(Face_1, L2, L2, L2)
sk2 = geompy.Sketcher3D()
sk2.addPointsAbsolute(L/3, L/3, 4*L/3) #Vertex_1 of Rombus
sk2.addPointsAbsolute(L, 0, L) #Vertex_2 of Rombus
sk2.addPointsAbsolute(2*L/3, 2*L/3, 2*L/3) #Vertex_3 of Rombus
#sk2.addPointsAbsolute(0, L, L) #Vertex_4 of Rombus
sk2.addPointsAbsolute(L/3, L/3, 4*L/3) #Vertex_1 of Rombus
a3D_Sketcher_2 = sk2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Extrusion_2 = geompy.MakePrismDXDYDZ(Face_2, L2, L2, L2)
Box_1 = geompy.MakeBoxDXDYDZ(a, a, L)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, Pi_4)
Translation_1 = geompy.MakeTranslation(Rotation_1, L, 0, 0)
Translation_2 = geompy.MakeTranslation(Translation_1, 0, 0, h4)
Box_2 = geompy.MakeBoxDXDYDZ(a2, a2, L)
Rotation_2 = geompy.MakeRotation(Box_2, OZ, Pi_4)
Translation_3 = geompy.MakeTranslation(Rotation_2, L, 0, h4)
Box_3 = geompy.MakeBoxDXDYDZ(a2, a2, L2)
Rotation_3 = geompy.MakeRotation(Box_3, OZ, Pi_4)
Translation_4 = geompy.MakeTranslation(Rotation_3, L, 0, h4)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( V_1, 'V_1' )
geompy.addToStudy( V_2, 'V_2' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Translation_1, 'Translation_1' )
geompy.addToStudy( Translation_2, 'Translation_2' )
geompy.addToStudy( Box_2, 'Box_2' )
geompy.addToStudy( Rotation_2, 'Rotation_2' )
geompy.addToStudy( Translation_3, 'Translation_3' )
geompy.addToStudy( Rotation_3, 'Rotation_3' )
geompy.addToStudy( Translation_4, 'Translation_4' )
#geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
#geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
#geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_1' )
#geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
Radius = r0/(1-alpha)
Sphere_1 = geompy.MakeSphereR(Radius)
Down = geompy.MakeMultiTranslation2D(Sphere_1, OX, L, 3, OY, L, 3)
Up = geompy.MakeTranslation(Down, 0, 0, 2*L)
Up_Down = geompy.MakeMultiTranslation1D(Down, OZ, L, 3)
Cut_1a = geompy.MakeCutList(Translation_1, [Up_Down])
Cut_1b = geompy.MakeCutList(Translation_2, [Up_Down])
Cut_2a = geompy.MakeCutList(Translation_3, [Up_Down])
Cut_2b = geompy.MakeCutList(Translation_4, [Up_Down])
Cut_3a = geompy.MakeCutList(Extrusion_1, [Up_Down])
Cut_3b = geompy.MakeCutList(Extrusion_2, [Up_Down])
Middle = geompy.MakeTranslation(Up_Down, L2, L2, L2)
Pore_1a = geompy.MakeCutList(Cut_1a, [Middle])
Pore_1b = geompy.MakeCutList(Cut_1b, [Middle])
Pore_2a = geompy.MakeCutList(Cut_2a, [Middle])
Pore_2b = geompy.MakeCutList(Cut_2b, [Middle])
Pore_3a = geompy.MakeCutList(Cut_3a, [Middle])
Pore_3b = geompy.MakeCutList(Cut_3b, [Middle])
geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Down, 'Down_' )
geompy.addToStudy( Up, 'Up_' )
geompy.addToStudy( Up_Down, 'Up_Down_' )
geompy.addToStudy( Cut_1a, 'Cut_1a_' )
geompy.addToStudy( Cut_1b, 'Cut_1b_' )
geompy.addToStudy( Middle, 'Middle_' )
geompy.addToStudy( Pore_1a, 'Pore_1a_' )
geompy.addToStudy( Pore_1b, 'Pore_1b_' )
geompy.addToStudy( Pore_2a, 'Pore_2a_' )
geompy.addToStudy( Pore_2b, 'Pore_2b_' )
geompy.addToStudy( Pore_3a, 'Pore_3a_' )
geompy.addToStudy( Pore_3b, 'Pore_3b_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(Middle, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
return grains, Pore_1a, Pore_3a

View File

@ -1,135 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def faceCenteredCubic(alpha):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
V = geompy.MakeVectorDXDYDZ(1, 1, 1)
V_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
V_2 = geompy.MakeVectorDXDYDZ(1, -1, 0)
L = 1.0
r0 = L*math.sqrt(2)/4
a=2*r0
b=L/2
#c=b/2
h=L/math.sqrt(3)
Pi_4 = 45*math.pi/180.0
Vertex_1 = geompy.MakeVertex(-a, -a, -b)
Vertex_2 = geompy.MakeVertex(a, a, b)
Vertex_3 = geompy.MakeVertex(-1, 0, -b) # Center of Sphere_1
Vertex_4 = geompy.MakeVertex(-b, 0, -0)
Vertex_5 = geompy.MakeVertex(0, 0, -b)
sk = geompy.Sketcher3D() # Rombus
sk.addPointsAbsolute(-b, -b, b) #Vertex_1 of Rombus
sk.addPointsAbsolute(0, -b, 0) #Vertex_2 of Rombus
sk.addPointsAbsolute(0, 0, -b) #Vertex_3 of Rombus
sk.addPointsAbsolute(-b, 0, 0) #Vertex_4 of Rombus
sk.addPointsAbsolute(-b, -b, b) #Vertex_1 of Rombus
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Face_1_Up = geompy.MakeTranslation(Face_1, b, b, 0)
Rhombohedron = geompy.MakeHexa2Faces(Face_1, Face_1_Up)
sk_2 = geompy.Sketcher3D() # Triangle
sk_2.addPointsAbsolute(-2*b/3, -2*b/3, b/3) #Vertex_1 of Triangle
sk_2.addPointsAbsolute(0, -b, 0) #Vertex_2 of Triangle
sk_2.addPointsAbsolute(-b/3, -b/3, -b/3) #Vertex_3 of Triangle
sk_2.addPointsAbsolute(-2*b/3, -2*b/3, b/3) #Vertex_1 of Triangle
a3D_Sketcher_2 = sk_2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Extrusion_2 = geompy.MakePrismVecH(Face_2, V, h)
sk_3 = geompy.Sketcher3D() # Hexagon
sk_3.addPointsAbsolute(-2*b/3, -2*b/3, b/3) #Vertex_1 of Hexagon
sk_3.addPointsAbsolute(0, -b, 0) #Vertex_2 of Hexagon
sk_3.addPointsAbsolute(b/3, -2*b/3, -2*b/3) #Vertex_3 of Hexagon
sk_3.addPointsAbsolute(0, 0, -b) #Vertex_4 of Hexagon
sk_3.addPointsAbsolute(-2*b/3, b/3, -2*b/3) #Vertex_5 of Hexagon
sk_3.addPointsAbsolute(-b, 0, 0) #Vertex_6 of Hexagon
sk_3.addPointsAbsolute(-2*b/3, -2*b/3, b/3) #Vertex_1 of Hexagon
a3D_Sketcher_3 = sk_3.wire()
Face_3 = geompy.MakeFaceWires([a3D_Sketcher_3], 1)
Extrusion_3 = geompy.MakePrismVecH(Face_3, V, h)
Box_1 = geompy.MakeBoxTwoPnt(Vertex_1, Vertex_2)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, Pi_4)
Box_2 = geompy.MakeBoxTwoPnt(Vertex_5, Vertex_2)
Rotation_2 = geompy.MakeRotation(Box_2, OZ, Pi_4)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( V_1, 'V_1' )
geompy.addToStudy( V_2, 'V_2' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Box_2, 'Box_2' )
geompy.addToStudy( Rotation_2, 'Rotation_2_' )
geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Face_1_Up, 'Face_1_Up' )
geompy.addToStudy( Rhombohedron, 'Rhombohedron' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
geompy.addToStudy( a3D_Sketcher_3, 'a3D_Sketcher_3' )
geompy.addToStudy( Face_3, 'Face_3' )
geompy.addToStudy( Extrusion_3, 'Extrusion_3' )
Radius = r0/(1-alpha)
Sphere_1 = geompy.MakeSpherePntR(Vertex_3, Radius)
Down = geompy.MakeMultiTranslation2D(Sphere_1, V_1, a, 3, V_2, a, 3)
Up_Down = geompy.MakeMultiTranslation1D(Down, OZ, 1, 2)
Cut_1 = geompy.MakeCutList(Rotation_1, [Up_Down], True)
Sphere_2 = geompy.MakeSpherePntR(Vertex_4, Radius)
Middle = geompy.MakeMultiTranslation2D(Sphere_2, V_1, a, 2, V_2, a, 2)
Cut_2 = geompy.MakeCutList(Cut_1, [Middle], True)
Common = geompy.MakeCommonList([Cut_2, Rotation_2], True)
Pore_3 = geompy.MakeCommonList([Rhombohedron, Cut_2], True)
Cut_3 = geompy.MakeCutList(Extrusion_2, [Up_Down], True)
Cut_4 = geompy.MakeCutList(Cut_3, [Middle], True)
Cut_5 = geompy.MakeCutList(Extrusion_3, [Up_Down], True)
Cut_6 = geompy.MakeCutList(Cut_5, [Middle], True)
#geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Down, 'Down_' )
geompy.addToStudy( Up_Down, 'Up_Down_' )
geompy.addToStudy( Cut_1, 'Cut_1_' )
geompy.addToStudy( Middle, 'Middle_' )
geompy.addToStudy( Cut_2, 'Pore_' )
geompy.addToStudy( Common, 'Pore_2_' )
geompy.addToStudy( Pore_3, 'Pore_3_' )
geompy.addToStudy( Cut_4, 'Pore_3a_' )
geompy.addToStudy( Cut_6, 'Pore_3b_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
grains += geompy.ExtractShapes(Middle, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
return grains, Common, Pore_3

View File

@ -1,110 +0,0 @@
import sys
import salome
salome.salome_init()
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
def simpleCubic(alpha):
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
r0 = 1.0
L = 2*r0
h = L
h2 = 2*h
d1 = L*math.sqrt(2)
d2 = L*math.sqrt(3)
d = 1/math.sqrt(3)
pi_4 = 45*math.pi/180.0
pi_2 = 90*math.pi/180.0
Box_1 = geompy.MakeBoxDXDYDZ(d1, d1, h)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, pi_4)
Translation_2 = geompy.MakeTranslation(Rotation_1, h, 0, 0)
Vertex_1 = geompy.MakeVertex(h, 0, 0)
Vertex_2 = geompy.MakeVertex(h, h, 0)
Vertex_3 = geompy.MakeVertex(h, h, h)
Line_1 = geompy.MakeLineTwoPnt(Vertex_2, Vertex_3)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, h2) #Vertex_1 of Rombus
sk.addPointsAbsolute(h, 0, h) #Vertex_2 of Rombus
sk.addPointsAbsolute(h, h, 0) #Vertex_3 of Rombus
sk.addPointsAbsolute(0, h, h) #Vertex_4 of Rombus
sk.addPointsAbsolute(0, 0, h2) #Vertex_1 of Rombus
sk_2 = geompy.Sketcher3D()
sk_2.addPointsAbsolute(L, L, L) #Vertex_1 of Hexagon
sk_2.addPointsAbsolute(5*L/3, 2*L/3, 2*L/3) #Vertex_2 of Hexagon
sk_2.addPointsAbsolute(2*L, L, 0) #Vertex_3 of Hexagon
sk_2.addPointsAbsolute(5*L/3, 5*L/3, -L/3) #Vertex_4 of Hexagon
sk_2.addPointsAbsolute(L, 2*L, 0) #Vertex_5 of Hexagon
sk_2.addPointsAbsolute(2*L/3, 5*L/3, 2*L/3) #Vertex_6 of Hexagon
sk_2.addPointsAbsolute(L, L, L) #Vertex_1 of Hexagon
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
Extrusion_1 = geompy.MakePrismVecH(Face_1, Vector_1, d1)
a3D_Sketcher_2 = sk_2.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
Vector_2 = geompy.MakeVectorDXDYDZ(1, 1, 1)
Extrusion_2 = geompy.MakePrismVecH(Face_2, Vector_2, d2/3)
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( Vertex_1, 'Vertex_1' )
geompy.addToStudy( Vertex_2, 'Vertex_2' )
geompy.addToStudy( Vertex_3, 'Vertex_3' )
geompy.addToStudy( Line_1, 'Line_1' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Translation_2, 'Translation_2' )
geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Vector_1, 'Vector_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Vector_2, 'Vector_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
Sphere_1 = geompy.MakeSphereR(r0/(1-alpha))
Multi_Translation_2 = geompy.MakeMultiTranslation2D(Sphere_1, OX, 2, 3, OY, 2, 3)
Multi_Translation_3 = geompy.MakeMultiTranslation1D(Multi_Translation_2, OZ, 2, 3)
Cut_1 = geompy.MakeCutList(Translation_2, [Multi_Translation_3])
Cut_2 = geompy.MakeCutList(Extrusion_1, [Multi_Translation_3])
Cut_3 = geompy.MakeCutList(Extrusion_2, [Multi_Translation_3])
Cut_V = geompy.MakeCutList(Cut_1, [Cut_3])
#Cut_2.SetColor(SALOMEDS.Color(0,0,1))
geompy.addToStudy( Sphere_1, 'Sphere_' )
geompy.addToStudy( Multi_Translation_2, 'Multi-Translation_2_' )
geompy.addToStudy( Multi_Translation_3, 'Multi-Translation_3_' )
geompy.addToStudy( Cut_1, 'Pore1_' )
geompy.addToStudy( Cut_2, 'Pore2_' )
geompy.addToStudy( Cut_3, 'Pore3_' )
geompy.addToStudy( Cut_V, 'Cut_V_' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()
# Preparation
#Cut_1 = geompy.MakeRotation(Cut_1, OZ, -pi_4)
grains = geompy.ExtractShapes(Multi_Translation_3, geompy.ShapeType["SOLID"], True)
grains = geompy.MakeFuseList(grains, False, False)
return Multi_Translation_3, Cut_1, Cut_2

View File

@ -1,84 +0,0 @@
from src import geometry_utils
import GEOM
from src import mesh_utils
import SMESH
from src import anisotropeCubic
import salome
import math
import samples
def genMesh(ctype, theta, flowdirection):
_G = globals()
#for g in _G:
func = _G.get(ctype)
if func:
salome.salome_init()
func(theta, flowdirection)
salome.salome_close()
else:
raise Exception("Unknown type of the sample function")
def simpleCubic(theta, flowdirection):
#radius = 1
#stackAngle = [0.5 * math.pi, 0.5 * math.pi, 0.5 * math.pi]
theta = theta if theta else 0.1
#layers = [3, 3, 3]
#grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
#scale = [2 * math.sqrt(2), 2 * math.sqrt(2), 2]
flowdirection = flowdirection if flowdirection else [1, 0, 0]
#style = 0
#cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
grains, cubic, _ = samples.simpleCubic(theta)
boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
fineness = 3
mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
mesh_utils.meshCompute(mesh)
def bodyCenteredCubic(theta, flowdirection):
#radius = 1
#stackAngle = [0.5 * math.pi, 0.25 * math.pi, 0.25 * math.pi]
theta = theta if theta else 0.1
#layers = [3, 3, 3]
#grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
#scale = [2 / math.sqrt(2), 2 / math.sqrt(2), 1]
flowdirection = flowdirection if flowdirection else [1, 0, 0]
#style = 0
#cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
grains, cubic, _ = samples.bodyCenteredCubic(theta)
boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
fineness = 3
mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
mesh_utils.meshCompute(mesh)
def faceCenteredCubic(theta, flowdirection):
#radius = 1
#stackAngle = [0.5 * math.pi, 0.5 * math.pi, 0.5 * math.pi]
theta = theta if theta else 0.1
#layers = [3, 3, 3]
#grains = anisotropeCubic.StructuredGrains(radius, stackAngle, theta, layers)
#scale = [1 / math.sqrt(2), 1 / math.sqrt(2), 1]
flowdirection = flowdirection if flowdirection else [1, 0, 0]
#style = 0
#cubic = anisotropeCubic.AnisotropeCubic(scale, grains, style)
grains, cubic, _ = samples.faceCenteredCubic(theta)
boundary = geometry_utils.boundaryCreate(cubic, flowdirection, grains)
fineness = 3
mesh = mesh_utils.meshCreate(cubic, boundary, fineness)
mesh_utils.meshCompute(mesh)

View File

@ -1,181 +0,0 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import salome, GEOM, SMESH, SALOMEDS
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math
import os, sys
import logging
import time
from datetime import timedelta
from .gutils import *
from .mutils import *
class simpleCubic:
def __init__(self, alpha, fillet, name = None):
self.name = name if type(name) != None else "simpleCubic"
self.scale = None
self.angle = None
self.pos = None
self.geometry = None
self.geometrybbox = None
self.mesh = None
self.boundary = None
self.geometry2 = None
self.geometry2bbox = None
self.grains = None
salome.salome_init()
"""
Create the simple cubic geometry.
Parameters:
alpha (float): Sphere intersection parameter which used for cutting spheres from box.
Radius = R0 / (1 - alpha)
Should be from 0.01 to 0.28
fillet (list): Fillet coefficient.
[fillet1, fillet2]
0 <= fillet <= 1
if fillet = [0, 0] then R_fillet = 0
Returns:
Configured geometry.
"""
geompy = geomBuilder.New()
# Parameters
R0 = 1
R = R0 / (1 - alpha)
C1 = fillet[0]
C2 = fillet[1]
alpha1 = 0.01
alpha2 = 0.28
Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
R_fillet = Cf * (R0 * math.sqrt(2) - R)
logging.info("geometryCreate: alpha = {}".format(alpha))
logging.info("geometryCreate: R_fillet = {}".format(R_fillet))
self.scale = [2 * math.sqrt(2), 2 * math.sqrt(2), 2]
self.angle = [0, 0, 0.25 * math.pi]
self.pos = [2, 0, 0]
# Box
box = geompy.MakeBoxDXDYDZ(scale[0], scale[1], scale[2])
box = rotate(box, angle)
box = geompy.MakeTranslation(box, pos[0], pos[1], pos[2])
self.geometrybbox = box
# Grains
stackang = [
0.5 * math.pi - stackAng[0],
0.5 * math.pi - stackAng[1],
0.5 * math.pi - stackAng[2]
]
xvec = geompy.MakeVector(
geompy.MakeVertex(0, 0, 0),
geompy.MakeVertex(1, 0, 0))
yvec = rotate(xvec, [0.5 * math.pi, 0, 0])
zvec = rotate(xvec, [0, 0.5 * math.pi, 0])
hvec = rotate(xvec, [stackang[0], 0, 0])
vvec = rotate(zvec, [0, stackang[1], stackang[2]])
grain = geompy.MakeSpherePntR(geompy.MakeVertex(pos[0], pos[1], pos[2]), R)
hstack = geompy.MakeMultiTranslation1D(grain, hvec, 2 * R0, 3)
vstack = geompy.MakeMultiTranslation1D(hstack, vvec, 2 * Ro, 3)
stack = geompy.MakeTranslation(vstack, -2 * R0, 0, 0)
self.grains = geompy.ExtractShapes(stack, geompy.ShapeType["SOLID"], True)
self.grains = geompy.MakeFuseList(self.grains, False, False)
if not R_fillet == 0:
self.grains = geompy.MakeFilletAll(self.grains, R_fillet)
# Geometry 1
self.geometry = geompy.MakeCutList(box, [self.grains], True)
# Rhombohedron
h = 2
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0, 0, h * 2)
sk.addPointsAbsolute(h, 0, h)
sk.addPointsAbsolute(h, h, 0)
sk.addPointsAbsolute(0, h, h)
sk.addPointsAbsolute(0, 0, h * 2)
rhombus = sk.wire()
rhombus = geompy.MakeFaceWires([rombus], 1)
vec = geompy.MakeVectorDXDYDZ(1, 1, 0)
rhombohedron = geompy.MakePrismVecH(rombus, vec, self.scale[0])
self.geometry2bbox = rhombohedron
# Geometry 2
self.geometry2 = geompy.MakeCutList(rhombohedron, [self.grains], True)
# Debug study
geompy.addToStudy(self.grains, "grains")
geompy.addToStudy(self.geometry, self.name)
geompy.addToStudy(self.geometrybbox, "bbox 1")
geompy.addToStudy(self.geometry2, "geometry 2")
geompy.addToStudy(self.geometry2bbox, "bbox 2")
if __name__ == "__main__":
# Arguments
buildpath = str(sys.argv[1])
alpha = float(sys.argv[2])
direction = str(sys.argv[3])
name = "simpleCubic-{}-{}".format(direction, alpha)
# Logger
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s: %(message)s",
handlers = [
logging.StreamHandler(),
logging.FileHandler(os.path.join(buildpath, "{}.log".format(name)))
])
start_time = time.monotonic()
# Simple cubic
logging.info("Creating the geometry ...")
sc = simpleCubic(alpha, [0, 0], name)
logging.info("Extracting boundaries ...")
boundary = sc.boundaryCreate(sc.geometry, direction, sc.grains)
logging.info("Creating the mesh ...")
mesh = sc.meshCreate(sc.geometry, boundary, 2) #, {
# "thickness": 0.001,
# "number": 1,
# "stretch": 1.1
#})
meshCompute(mesh)
logging.info("Exporting the mesh ...")
meshExport(mesh, buildpath)
end_time = time.monotonic()
logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time)))
logging.info("Done.")
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

View File

@ -1,314 +0,0 @@
#!/usr/bin/env python
###
### This file is generated automatically by SALOME v9.6.0 with dump python functionality
###
import sys
import salome
salome.salome_init()
import salome_notebook
notebook = salome_notebook.NoteBook()
sys.path.insert(0, r'/home/nafaryus/projects/anisotrope-cube/temp')
###
### GEOM component
###
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
geomObj_1 = geompy.MakeVertex(0, 0, 0)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(-0.6666667, -0.6666667, -0.1666667)
sk.addPointsAbsolute(-0.3333333, -0.8333333, -0.3333333)
sk.addPointsAbsolute(-0.1666667, -0.6666667, -0.6666667)
sk.addPointsAbsolute(-0.3333333, -0.3333333, -0.8333333)
sk.addPointsAbsolute(-0.6666667, -0.1666667, -0.6666667)
sk.addPointsAbsolute(-0.8333333, -0.3333333, -0.3333333)
sk.addPointsAbsolute(-0.6666667, -0.6666667, -0.1666667)
geomObj_2 = sk.wire()
geomObj_3 = geompy.MakeFaceWires([geomObj_2], 0)
geomObj_4 = geompy.GetNormal(geomObj_3)
geomObj_5 = geompy.MakePrismVecH(geomObj_3, geomObj_4, 1.732050807568877)
geomObj_6 = geompy.MakeScaleTransform(geomObj_3, geomObj_1, 100)
geomObj_7 = geompy.MakeScaleTransform(geomObj_5, geomObj_1, 100)
[geomObj_8,geomObj_9,geomObj_10,geomObj_11,geomObj_12,geomObj_13,geomObj_14,geomObj_15] = geompy.ExtractShapes(geomObj_7, geompy.ShapeType["FACE"], False)
geomObj_16 = geompy.GetNormal(geomObj_8)
geomObj_17 = geompy.GetNormal(geomObj_9)
geomObj_18 = geompy.GetNormal(geomObj_10)
geomObj_19 = geompy.GetNormal(geomObj_11)
geomObj_20 = geompy.GetNormal(geomObj_12)
geomObj_21 = geompy.GetNormal(geomObj_13)
geomObj_22 = geompy.GetNormal(geomObj_14)
geomObj_23 = geompy.GetNormal(geomObj_15)
geomObj_24 = geompy.MakeVectorDXDYDZ(1, 0, 0)
geomObj_25 = geompy.MakeVectorDXDYDZ(0, 1, 0)
geomObj_26 = geompy.MakeVectorDXDYDZ(0, 0, 1)
geomObj_27 = geompy.MakeVectorDXDYDZ(1, 1, 0)
geomObj_28 = geompy.MakeVectorDXDYDZ(1, -1, 0)
geomObj_29 = geompy.MakeVertex(-1, 0, -0.5)
geomObj_30 = geompy.MakeSpherePntR(geomObj_29, 0.3761206282907168)
geomObj_31 = geompy.MakeMultiTranslation2D(geomObj_30, geomObj_27, 0.7071067811865476, 3, geomObj_28, 0.7071067811865476, 3)
geomObj_32 = geompy.MakeMultiTranslation1D(geomObj_31, geomObj_26, 1, 2)
[geomObj_33,geomObj_34,geomObj_35,geomObj_36,geomObj_37,geomObj_38,geomObj_39,geomObj_40,geomObj_41,geomObj_42,geomObj_43,geomObj_44,geomObj_45,geomObj_46,geomObj_47,geomObj_48,geomObj_49,geomObj_50] = geompy.ExtractShapes(geomObj_32, geompy.ShapeType["SOLID"], True)
geomObj_51 = geompy.MakeVertex(-1.5, 0, -1)
geomObj_52 = geompy.MakeSpherePntR(geomObj_51, 0.3761206282907168)
geomObj_53 = geompy.MakeMultiTranslation2D(geomObj_52, geomObj_27, 0.7071067811865476, 4, geomObj_28, 0.7071067811865476, 4)
geomObj_54 = geompy.MakeMultiTranslation1D(geomObj_53, geomObj_26, 1, 3)
[geomObj_55,geomObj_56,geomObj_57,geomObj_58,geomObj_59,geomObj_60,geomObj_61,geomObj_62,geomObj_63,geomObj_64,geomObj_65,geomObj_66,geomObj_67,geomObj_68,geomObj_69,geomObj_70,geomObj_71,geomObj_72,geomObj_73,geomObj_74,geomObj_75,geomObj_76,geomObj_77,geomObj_78,geomObj_79,geomObj_80,geomObj_81,geomObj_82,geomObj_83,geomObj_84,geomObj_85,geomObj_86,geomObj_87,geomObj_88,geomObj_89,geomObj_90,geomObj_91,geomObj_92,geomObj_93,geomObj_94,geomObj_95,geomObj_96,geomObj_97,geomObj_98,geomObj_99,geomObj_100,geomObj_101,geomObj_102] = geompy.ExtractShapes(geomObj_54, geompy.ShapeType["SOLID"], True)
geomObj_103 = geompy.MakeFuseList([geomObj_33, geomObj_34, geomObj_35, geomObj_36, geomObj_37, geomObj_38, geomObj_39, geomObj_40, geomObj_41, geomObj_42, geomObj_43, geomObj_44, geomObj_45, geomObj_46, geomObj_47, geomObj_48, geomObj_49, geomObj_50, geomObj_55, geomObj_56, geomObj_57, geomObj_58, geomObj_59, geomObj_60, geomObj_61, geomObj_62, geomObj_63, geomObj_64, geomObj_65, geomObj_66, geomObj_67, geomObj_68, geomObj_69, geomObj_70, geomObj_71, geomObj_72, geomObj_73, geomObj_74, geomObj_75, geomObj_76, geomObj_77, geomObj_78, geomObj_79, geomObj_80, geomObj_81, geomObj_82, geomObj_83, geomObj_84, geomObj_85, geomObj_86, geomObj_87, geomObj_88, geomObj_89, geomObj_90, geomObj_91, geomObj_92, geomObj_93, geomObj_94, geomObj_95, geomObj_96, geomObj_97, geomObj_98, geomObj_99, geomObj_100, geomObj_101, geomObj_102], False, False)
geomObj_104 = geompy.MakeScaleTransform(geomObj_103, geomObj_1, 100)
geomObj_105 = geompy.MakeFilletAll(geomObj_104, 0.6170130261493888)
geomObj_106 = geompy.MakeCutList(geomObj_7, [geomObj_105])
faceCenteredCubic = geompy.MakeScaleTransform(geomObj_106, geomObj_1, 0.01)
geomObj_107 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geompy.UnionIDs(geomObj_107, [3, 21, 30, 35, 42, 59, 84, 89, 118, 131, 136, 177, 180, 183, 188, 203, 230, 233, 236, 247, 259, 261, 266, 270, 272, 277, 281, 296, 301, 304, 306, 311, 316, 322, 324, 338, 353, 357, 359, 362, 365, 370, 375, 380, 383, 388, 391, 398, 421, 424, 427, 438, 443, 448, 455, 460, 465, 470, 475, 480, 487, 492, 497, 504, 517, 534, 557, 561, 563, 566, 569, 577, 579, 583, 586, 588, 591, 595, 598, 602, 604])
inlet = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_108 = geompy.MakeCutList(geomObj_6, [geomObj_105])
geomObj_109 = geompy.MakeScaleTransform(geomObj_108, geomObj_1, 0.01)
geomObj_110 = geompy.GetInPlace(faceCenteredCubic, geomObj_109, True)
[geomObj_111,geomObj_112,geomObj_113,geomObj_114,geomObj_115,geomObj_116] = geompy.SubShapeAll(geomObj_110, geompy.ShapeType["FACE"])
geompy.UnionList(inlet, [geomObj_111, geomObj_112, geomObj_113, geomObj_114, geomObj_115, geomObj_116])
outlet = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_117 = geompy.MakeCutList(geomObj_15, [geomObj_105])
geomObj_118 = geompy.MakeScaleTransform(geomObj_117, geomObj_1, 0.01)
geomObj_119 = geompy.GetInPlace(faceCenteredCubic, geomObj_118, True)
[geomObj_120,geomObj_121,geomObj_122,geomObj_123,geomObj_124,geomObj_125] = geompy.SubShapeAll(geomObj_119, geompy.ShapeType["FACE"])
geompy.UnionList(outlet, [geomObj_120, geomObj_121, geomObj_122, geomObj_123, geomObj_124, geomObj_125])
symetry0 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_126 = geompy.MakeCutList(geomObj_8, [geomObj_105])
geomObj_127 = geompy.MakeScaleTransform(geomObj_126, geomObj_1, 0.01)
geomObj_128 = geompy.GetInPlace(faceCenteredCubic, geomObj_127, True)
[geomObj_129,geomObj_130] = geompy.SubShapeAll(geomObj_128, geompy.ShapeType["FACE"])
geompy.UnionList(symetry0, [geomObj_129, geomObj_130])
symetry1 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_131 = geompy.MakeCutList(geomObj_9, [geomObj_105])
geomObj_132 = geompy.MakeScaleTransform(geomObj_131, geomObj_1, 0.01)
geomObj_133 = geompy.GetInPlace(faceCenteredCubic, geomObj_132, True)
[geomObj_134,geomObj_135] = geompy.SubShapeAll(geomObj_133, geompy.ShapeType["FACE"])
geompy.UnionList(symetry1, [geomObj_134, geomObj_135])
symetry2 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_136 = geompy.MakeCutList(geomObj_10, [geomObj_105])
geomObj_137 = geompy.MakeScaleTransform(geomObj_136, geomObj_1, 0.01)
geomObj_138 = geompy.GetInPlace(faceCenteredCubic, geomObj_137, True)
[geomObj_139,geomObj_140] = geompy.SubShapeAll(geomObj_138, geompy.ShapeType["FACE"])
geompy.UnionList(symetry2, [geomObj_139, geomObj_140])
symetry3 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_141 = geompy.MakeCutList(geomObj_11, [geomObj_105])
geomObj_142 = geompy.MakeScaleTransform(geomObj_141, geomObj_1, 0.01)
geomObj_143 = geompy.GetInPlace(faceCenteredCubic, geomObj_142, True)
[geomObj_144,geomObj_145] = geompy.SubShapeAll(geomObj_143, geompy.ShapeType["FACE"])
geompy.UnionList(symetry3, [geomObj_144, geomObj_145])
symetry4 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_146 = geompy.MakeCutList(geomObj_12, [geomObj_105])
geomObj_147 = geompy.MakeScaleTransform(geomObj_146, geomObj_1, 0.01)
geomObj_148 = geompy.GetInPlace(faceCenteredCubic, geomObj_147, True)
[geomObj_149,geomObj_150] = geompy.SubShapeAll(geomObj_148, geompy.ShapeType["FACE"])
geompy.UnionList(symetry4, [geomObj_149, geomObj_150])
symetry5 = geompy.CreateGroup(faceCenteredCubic, geompy.ShapeType["FACE"])
geomObj_151 = geompy.MakeCutList(geomObj_13, [geomObj_105])
geomObj_152 = geompy.MakeScaleTransform(geomObj_151, geomObj_1, 0.01)
geomObj_153 = geompy.GetInPlace(faceCenteredCubic, geomObj_152, True)
[geomObj_154,geomObj_155] = geompy.SubShapeAll(geomObj_153, geompy.ShapeType["FACE"])
geompy.UnionList(symetry5, [geomObj_154, geomObj_155])
wall = geompy.CutListOfGroups([geomObj_107], [inlet, outlet, symetry0, symetry1, symetry2, symetry3, symetry4, symetry5])
[geomObj_107, inlet, geomObj_110, outlet, geomObj_119, symetry0, geomObj_128, symetry1, geomObj_133, symetry2, geomObj_138, symetry3, geomObj_143, symetry4, geomObj_148, symetry5, geomObj_153, wall] = geompy.GetExistingSubObjects(faceCenteredCubic, False)
[geomObj_107, inlet, geomObj_110, outlet, geomObj_119, symetry0, geomObj_128, symetry1, geomObj_133, symetry2, geomObj_138, symetry3, geomObj_143, symetry4, geomObj_148, symetry5, geomObj_153, wall] = geompy.GetExistingSubObjects(faceCenteredCubic, False)
[geomObj_107, inlet, geomObj_110, outlet, geomObj_119, symetry0, geomObj_128, symetry1, geomObj_133, symetry2, geomObj_138, symetry3, geomObj_143, symetry4, geomObj_148, symetry5, geomObj_153, wall] = geompy.GetExistingSubObjects(faceCenteredCubic, False)
geompy.addToStudy( faceCenteredCubic, 'faceCenteredCubic' )
geompy.addToStudyInFather( faceCenteredCubic, inlet, 'inlet' )
geompy.addToStudyInFather( faceCenteredCubic, outlet, 'outlet' )
geompy.addToStudyInFather( faceCenteredCubic, symetry0, 'symetry0' )
geompy.addToStudyInFather( faceCenteredCubic, symetry1, 'symetry1' )
geompy.addToStudyInFather( faceCenteredCubic, symetry2, 'symetry2' )
geompy.addToStudyInFather( faceCenteredCubic, symetry3, 'symetry3' )
geompy.addToStudyInFather( faceCenteredCubic, symetry4, 'symetry4' )
geompy.addToStudyInFather( faceCenteredCubic, symetry5, 'symetry5' )
geompy.addToStudyInFather( faceCenteredCubic, wall, 'wall' )
###
### SMESH component
###
import SMESH, SALOMEDS
from salome.smesh import smeshBuilder
smesh = smeshBuilder.New()
#smesh.SetEnablePublish( False ) # Set to False to avoid publish in study if not needed or in some particular situations:
# multiples meshes built in parallel, complex and numerous mesh edition (performance)
Mesh_1 = smesh.Mesh(faceCenteredCubic)
NETGEN_1D_2D_3D = Mesh_1.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
NETGEN_3D_Parameters_1 = NETGEN_1D_2D_3D.Parameters()
NETGEN_3D_Parameters_1.SetMaxSize( 0.05 )
NETGEN_3D_Parameters_1.SetMinSize( 0.005 )
NETGEN_3D_Parameters_1.SetSecondOrder( 0 )
NETGEN_3D_Parameters_1.SetOptimize( 1 )
NETGEN_3D_Parameters_1.SetFineness( 1 )
NETGEN_3D_Parameters_1.SetChordalError( -1 )
NETGEN_3D_Parameters_1.SetChordalErrorEnabled( 0 )
NETGEN_3D_Parameters_1.SetUseSurfaceCurvature( 1 )
NETGEN_3D_Parameters_1.SetFuseEdges( 1 )
NETGEN_3D_Parameters_1.SetQuadAllowed( 0 )
NETGEN_3D_Parameters_1.SetCheckChartBoundary( 88 )
Viscous_Layers_1 = NETGEN_1D_2D_3D.ViscousLayers(0.001,2,1.2,[ 21, 35, 316, 183, 365, 380, 475, 460, 497, 448, 595, 579 ],1,smeshBuilder.SURF_OFFSET_SMOOTH)
#Group_1 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
inlet_1 = Mesh_1.GroupOnGeom(inlet,'inlet',SMESH.FACE)
#Group_3 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
outlet_1 = Mesh_1.GroupOnGeom(outlet,'outlet',SMESH.FACE)
#Group_5 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry0_1 = Mesh_1.GroupOnGeom(symetry0,'symetry0',SMESH.FACE)
#Group_7 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry1_1 = Mesh_1.GroupOnGeom(symetry1,'symetry1',SMESH.FACE)
#Group_9 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry2_1 = Mesh_1.GroupOnGeom(symetry2,'symetry2',SMESH.FACE)
#Group_11 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry3_1 = Mesh_1.GroupOnGeom(symetry3,'symetry3',SMESH.FACE)
#Group_13 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry4_1 = Mesh_1.GroupOnGeom(symetry4,'symetry4',SMESH.FACE)
#Group_15 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry5_1 = Mesh_1.GroupOnGeom(symetry5,'symetry5',SMESH.FACE)
#Group_17 = Mesh_1.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
wall_1 = Mesh_1.GroupOnGeom(wall,'wall',SMESH.FACE)
isDone = Mesh_1.Compute()
[ Group_1, inlet_1, Group_3, outlet_1, Group_5, symetry0_1, Group_7, symetry1_1, Group_9, symetry2_1, Group_11, symetry3_1, Group_13, symetry4_1, Group_15, symetry5_1, Group_17, wall_1 ] = Mesh_1.GetGroups()
Mesh_4 = smesh.Mesh()
Mesh_2 = smesh.Mesh(faceCenteredCubic)
status = Mesh_2.AddHypothesis(NETGEN_3D_Parameters_1)
status = Mesh_2.AddHypothesis(Viscous_Layers_1)
NETGEN_1D_2D_3D_1 = Mesh_2.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
#Group_1_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
inlet_2 = Mesh_2.GroupOnGeom(inlet,'inlet',SMESH.FACE)
#Group_3_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
outlet_2 = Mesh_2.GroupOnGeom(outlet,'outlet',SMESH.FACE)
#Group_5_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry0_2 = Mesh_2.GroupOnGeom(symetry0,'symetry0',SMESH.FACE)
#Group_7_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry1_2 = Mesh_2.GroupOnGeom(symetry1,'symetry1',SMESH.FACE)
#Group_9_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry2_2 = Mesh_2.GroupOnGeom(symetry2,'symetry2',SMESH.FACE)
#Group_11_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry3_2 = Mesh_2.GroupOnGeom(symetry3,'symetry3',SMESH.FACE)
#Group_13_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry4_2 = Mesh_2.GroupOnGeom(symetry4,'symetry4',SMESH.FACE)
#Group_15_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
symetry5_2 = Mesh_2.GroupOnGeom(symetry5,'symetry5',SMESH.FACE)
#Group_17_1 = Mesh_2.GroupOnGeom(__NOT__Published__Object__,'',SMESH.FACE)
wall_2 = Mesh_2.GroupOnGeom(wall,'wall',SMESH.FACE)
isDone = Mesh_2.Compute()
[ Group_1_1, inlet_2, Group_3_1, outlet_2, Group_5_1, symetry0_2, Group_7_1, symetry1_2, Group_9_1, symetry2_2, Group_11_1, symetry3_2, Group_13_1, symetry4_2, Group_15_1, symetry5_2, Group_17_1, wall_2 ] = Mesh_2.GetGroups()
[ Group_1, inlet_1, Group_3, outlet_1, Group_5, symetry0_1, Group_7, symetry1_1, Group_9, symetry2_1, Group_11, symetry3_1, Group_13, symetry4_1, Group_15, symetry5_1, Group_17, wall_1 ] = Mesh_1.GetGroups()
[ Group_1_1, inlet_2, Group_3_1, outlet_2, Group_5_1, symetry0_2, Group_7_1, symetry1_2, Group_9_1, symetry2_2, Group_11_1, symetry3_2, Group_13_1, symetry4_2, Group_15_1, symetry5_2, Group_17_1, wall_2 ] = Mesh_2.GetGroups()
aCriteria = []
aCriterion = smesh.GetCriterion(SMESH.VOLUME,SMESH.FT_ElemGeomType,SMESH.FT_Undefined,SMESH.Geom_PYRAMID)
aCriteria.append(aCriterion)
[ Group_1_1, inlet_2, Group_3_1, outlet_2, Group_5_1, symetry0_2, Group_7_1, symetry1_2, Group_9_1, symetry2_2, Group_11_1, symetry3_2, Group_13_1, symetry4_2, Group_15_1, symetry5_2, Group_17_1, wall_2 ] = Mesh_2.GetGroups()
Mesh_1.Clear()
Mesh_2.Clear()
Mesh_3 = smesh.Mesh(faceCenteredCubic)
status = Mesh_3.AddHypothesis(NETGEN_3D_Parameters_1)
status = Mesh_3.AddHypothesis(Viscous_Layers_1)
NETGEN_1D_2D_3D_2 = Mesh_3.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
inlet_3 = Mesh_3.GroupOnGeom(inlet,'inlet',SMESH.FACE)
outlet_3 = Mesh_3.GroupOnGeom(outlet,'outlet',SMESH.FACE)
symetry0_3 = Mesh_3.GroupOnGeom(symetry0,'symetry0',SMESH.FACE)
symetry1_3 = Mesh_3.GroupOnGeom(symetry1,'symetry1',SMESH.FACE)
symetry2_3 = Mesh_3.GroupOnGeom(symetry2,'symetry2',SMESH.FACE)
symetry3_3 = Mesh_3.GroupOnGeom(symetry3,'symetry3',SMESH.FACE)
symetry4_3 = Mesh_3.GroupOnGeom(symetry4,'symetry4',SMESH.FACE)
symetry5_3 = Mesh_3.GroupOnGeom(symetry5,'symetry5',SMESH.FACE)
wall_3 = Mesh_3.GroupOnGeom(wall,'wall',SMESH.FACE)
isDone = Mesh_3.Compute()
[ smeshObj_1, inlet_3, smeshObj_2, outlet_3, smeshObj_3, symetry0_3, smeshObj_4, symetry1_3, smeshObj_5, symetry2_3, smeshObj_6, symetry3_3, smeshObj_7, symetry4_3, smeshObj_8, symetry5_3, smeshObj_9, wall_3 ] = Mesh_3.GetGroups()
aCriteria = []
aCriterion = smesh.GetCriterion(SMESH.VOLUME,SMESH.FT_ElemGeomType,SMESH.FT_Undefined,SMESH.Geom_PYRAMID)
aCriteria.append(aCriterion)
Mesh_3.SplitVolumesIntoTetra( Mesh_3.GetIDSource([ 189580, 222102, 189581, 222103, 146658, 146659, 106364, 106365, 146364, 146365, 105266, 105267, 162840, 162841, 189808, 189809, 105172, 105173, 180526, 180527, 154066, 154067, 104882, 104883, 123066, 123067, 123068, 123069, 180640, 222050, 180641, 222051, 162664, 180238, 162665, 180239, 106312, 106313, 154188, 154189, 122780, 122781, 162784, 146726, 162785, 222274, 146727, 162786, 222275, 162787, 189756, 189757, 189458, 189459, 222394, 222395 ], SMESH.VOLUME), 1 )
[ smeshObj_1, inlet_3, smeshObj_2, outlet_3, smeshObj_3, symetry0_3, smeshObj_4, symetry1_3, smeshObj_5, symetry2_3, smeshObj_6, symetry3_3, smeshObj_7, symetry4_3, smeshObj_8, symetry5_3, smeshObj_9, wall_3 ] = Mesh_3.GetGroups()
## some objects were removed
aStudyBuilder = salome.myStudy.NewBuilder()
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_8))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_9))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_6))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_7))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_5))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_3))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_4))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_1))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
SO = salome.myStudy.FindObjectIOR(salome.myStudy.ConvertObjectToIOR(smeshObj_2))
if SO: aStudyBuilder.RemoveObjectWithChildren(SO)
## Set names of Mesh objects
smesh.SetName(NETGEN_1D_2D_3D.GetAlgorithm(), 'NETGEN 1D-2D-3D')
smesh.SetName(symetry1_1, 'symetry1')
smesh.SetName(Group_9, 'Group_9')
smesh.SetName(Viscous_Layers_1, 'Viscous Layers_1')
smesh.SetName(NETGEN_3D_Parameters_1, 'NETGEN 3D Parameters_1')
smesh.SetName(Group_1, 'Group_1')
smesh.SetName(inlet_1, 'inlet')
smesh.SetName(Group_3, 'Group_3')
smesh.SetName(outlet_1, 'outlet')
smesh.SetName(Group_5, 'Group_5')
smesh.SetName(symetry0_1, 'symetry0')
smesh.SetName(Group_7, 'Group_7')
smesh.SetName(Mesh_1.GetMesh(), 'Mesh_1')
smesh.SetName(Mesh_2.GetMesh(), 'Mesh_2')
smesh.SetName(Mesh_4.GetMesh(), 'Mesh_4')
smesh.SetName(Mesh_3.GetMesh(), 'Mesh_3')
smesh.SetName(wall_3, 'wall')
smesh.SetName(symetry4_3, 'symetry4')
smesh.SetName(symetry5_3, 'symetry5')
smesh.SetName(symetry2_3, 'symetry2')
smesh.SetName(wall_1, 'wall')
smesh.SetName(symetry3_3, 'symetry3')
smesh.SetName(Group_11, 'Group_11')
smesh.SetName(symetry2_1, 'symetry2')
smesh.SetName(Group_13, 'Group_13')
smesh.SetName(symetry3_1, 'symetry3')
smesh.SetName(Group_15, 'Group_15')
smesh.SetName(symetry1_2, 'symetry1')
smesh.SetName(symetry4_1, 'symetry4')
smesh.SetName(symetry1_3, 'symetry1')
smesh.SetName(Group_9_1, 'Group_9')
smesh.SetName(Group_17, 'Group_17')
smesh.SetName(symetry0_2, 'symetry0')
smesh.SetName(symetry5_1, 'symetry5')
smesh.SetName(Group_7_1, 'Group_7')
smesh.SetName(outlet_2, 'outlet')
smesh.SetName(outlet_3, 'outlet')
smesh.SetName(Group_5_1, 'Group_5')
smesh.SetName(inlet_2, 'inlet')
smesh.SetName(symetry0_3, 'symetry0')
smesh.SetName(Group_3_1, 'Group_3')
smesh.SetName(Group_1_1, 'Group_1')
smesh.SetName(inlet_3, 'inlet')
smesh.SetName(wall_2, 'wall')
smesh.SetName(Group_17_1, 'Group_17')
smesh.SetName(symetry5_2, 'symetry5')
smesh.SetName(Group_15_1, 'Group_15')
smesh.SetName(symetry4_2, 'symetry4')
smesh.SetName(Group_13_1, 'Group_13')
smesh.SetName(symetry3_2, 'symetry3')
smesh.SetName(Group_11_1, 'Group_11')
smesh.SetName(symetry2_2, 'symetry2')
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

View File

@ -1,175 +0,0 @@
#!/usr/bin/env python
###
### This file is generated automatically by SALOME v9.6.0 with dump python functionality
###
import sys
import salome
salome.salome_init()
import salome_notebook
notebook = salome_notebook.NoteBook()
sys.path.insert(0, r'/home/nafaryus/projects/anisotrope-cube/build')
###
### GEOM component
###
import GEOM
from salome.geom import geomBuilder
import math
import SALOMEDS
geompy = geomBuilder.New()
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
Box_1 = geompy.MakeBoxDXDYDZ(2.82842712474619, 2.82842712474619, 2)
Rotation_1 = geompy.MakeRotation(Box_1, OZ, 45*math.pi/180.0)
Translation_2 = geompy.MakeTranslation(Rotation_1, 2, 0, 0)
Vertex_1 = geompy.MakeVertex(2, 0, 0)
Vertex_2 = geompy.MakeVertex(2, 2, 0)
Vertex_3 = geompy.MakeVertex(2, 2, 2)
Line_1 = geompy.MakeLineTwoPnt(Vertex_2, Vertex_3)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(0.0000000, 0.0000000, 4.0000000)
sk.addPointsAbsolute(2.0000000, 0.0000000, 2.0000000)
sk.addPointsAbsolute(2.0000000, 2.0000000, 0.0000000)
sk.addPointsAbsolute(0.0000000, 2.0000000, 2.0000000)
sk.addPointsAbsolute(0.0000000, 0.0000000, 4.0000000)
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
Extrusion_1 = geompy.MakePrismVecH(Face_1, Vector_1, 2.82842712474619)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(2.0000000, 2.0000000, 2.0000000)
sk.addPointsAbsolute(3.3333333, 1.3333333, 1.3333333)
sk.addPointsAbsolute(4.0000000, 2.0000000, 0.0000000)
sk.addPointsAbsolute(3.3333333, 3.3333333, -0.6666667)
sk.addPointsAbsolute(2.0000000, 4.0000000, 0.0000000)
sk.addPointsAbsolute(1.3333333, 3.3333333, 1.3333333)
sk.addPointsAbsolute(2.0000000, 2.0000000, 2.0000000)
a3D_Sketcher_2 = sk.wire()
Face_2 = geompy.MakeFaceWires([a3D_Sketcher_2], 1)
geomObj_1 = geompy.MakeVectorDXDYDZ(1, 1, 1)
Extrusion_2 = geompy.MakePrismVecH(Face_2, geomObj_1, 1.154700538379251)
a3D_Sketcher_3 = geompy.MakeTranslation(a3D_Sketcher_2, -2.333333333333333, -2.333333333333333, -0.3333333333333333)
geomObj_2 = geompy.MakeFaceWires([a3D_Sketcher_3], 1)
Vector_2 = geompy.MakeVectorDXDYDZ(1, 1, 1)
Extrusion_3 = geompy.MakePrismVecH(geomObj_2, Vector_2, 3.464101615137754)
Sphere_ = geompy.MakeSphereR(1.111111111111111)
Multi_Translation_2_ = geompy.MakeMultiTranslation2D(Sphere_, OX, 2, 3, OY, 2, 3)
Multi_Translation_3_ = geompy.MakeMultiTranslation1D(Multi_Translation_2_, OZ, 2, 3)
[geomObj_3,geomObj_4,geomObj_5,geomObj_6,geomObj_7,geomObj_8,geomObj_9,geomObj_10,geomObj_11,geomObj_12,geomObj_13,geomObj_14,geomObj_15,geomObj_16,geomObj_17,geomObj_18,geomObj_19,geomObj_20,geomObj_21,geomObj_22,geomObj_23,geomObj_24,geomObj_25,geomObj_26,geomObj_27,geomObj_28,geomObj_29] = geompy.ExtractShapes(Multi_Translation_3_, geompy.ShapeType["SOLID"], True)
Pore1_ = geompy.MakeCutList(Translation_2, [Multi_Translation_3_])
Pore2_ = geompy.MakeCutList(Extrusion_1, [Multi_Translation_3_])
Pore3_ = geompy.MakeCutList(Extrusion_2, [Multi_Translation_3_])
Cut_V_ = geompy.MakeCutList(Pore1_, [Pore3_])
Pore4_ = geompy.MakeCutList(Extrusion_3, [Multi_Translation_3_])
grains = geompy.MakeFuseList([geomObj_3, geomObj_4, geomObj_5, geomObj_6, geomObj_7, geomObj_8, geomObj_9, geomObj_10, geomObj_11, geomObj_12, geomObj_13, geomObj_14, geomObj_15, geomObj_16, geomObj_17, geomObj_18, geomObj_19, geomObj_20, geomObj_21, geomObj_22, geomObj_23, geomObj_24, geomObj_25, geomObj_26, geomObj_27, geomObj_28, geomObj_29], False, False)
geomObj_30 = geompy.MakeMarker(0, 0, 0, 1, 0, 0, 0, 1, 0)
sk = geompy.Sketcher3D()
sk.addPointsAbsolute(2, 0, 0)
sk.addPointsAbsolute(0, 2, 0)
sk.addPointsAbsolute(0, 2, 2)
sk.addPointsAbsolute(2, 0, 2)
sk.addPointsAbsolute(2, 0, 0)
a3D_Sketcher_1_1 = sk.wire()
Face_3 = geompy.MakeFaceWires([a3D_Sketcher_1_1], 0)
Vector_Normal_1 = geompy.GetNormal(Face_3)
cuban = geompy.MakePrismVecH(Face_3, Vector_Normal_1, 2)
inlet = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(inlet, [31])
outlet = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(outlet, [33])
sp1 = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(sp1, [20])
sp2 = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(sp2, [3])
sp3 = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(sp3, [27])
sp4 = geompy.CreateGroup(cuban, geompy.ShapeType["FACE"])
geompy.UnionIDs(sp4, [13])
Cut_1 = geompy.MakeCutList(cuban, [grains])
Cut_2 = geompy.MakeCutList(inlet, [grains])
[geomObj_31] = geompy.SubShapeAll(Cut_2, geompy.ShapeType["FACE"])
Cut_3 = geompy.MakeCutList(outlet, [grains])
Cut_4 = geompy.MakeCutList(sp1, [grains])
Cut_5 = geompy.MakeCutList(sp2, [grains])
Cut_6 = geompy.MakeCutList(sp3, [grains])
Cut_7 = geompy.MakeCutList(sp4, [grains])
geomObj_32 = geompy.GetInPlace(Cut_1, Cut_2, True)
[geomObj_33] = geompy.SubShapeAll(geomObj_32, geompy.ShapeType["FACE"])
[geomObj_34] = geompy.SubShapeAll(geomObj_32, geompy.ShapeType["FACE"])
inlet_1 = geompy.CreateGroup(Cut_1, geompy.ShapeType["FACE"])
geompy.UnionIDs(inlet_1, [13])
geomObj_35 = geompy.GetInPlace(Cut_1, Cut_3, True)
[geomObj_36] = geompy.SubShapeAll(geomObj_35, geompy.ShapeType["FACE"])
[geomObj_37] = geompy.SubShapeAll(geomObj_35, geompy.ShapeType["FACE"])
outlet_1 = geompy.CreateGroup(Cut_1, geompy.ShapeType["FACE"])
geompy.UnionIDs(outlet_1, [101])
geometry1 = Pore1_
geometry2 = Pore4_
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( Box_1, 'Box_1' )
geompy.addToStudy( Rotation_1, 'Rotation_1' )
geompy.addToStudy( Translation_2, 'Translation_2' )
geompy.addToStudy( Vertex_1, 'Vertex_1' )
geompy.addToStudy( Vertex_2, 'Vertex_2' )
geompy.addToStudy( Vertex_3, 'Vertex_3' )
geompy.addToStudy( Line_1, 'Line_1' )
geompy.addToStudy( a3D_Sketcher_1, 'a3D_Sketcher_1' )
geompy.addToStudy( Face_1, 'Face_1' )
geompy.addToStudy( Vector_1, 'Vector_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
geompy.addToStudy( a3D_Sketcher_2, 'a3D_Sketcher_2' )
geompy.addToStudy( Face_2, 'Face_2' )
geompy.addToStudy( Extrusion_2, 'Extrusion_2' )
geompy.addToStudy( a3D_Sketcher_3, 'a3D_Sketcher_3' )
geompy.addToStudy( Vector_2, 'Vector_2' )
geompy.addToStudy( Extrusion_3, 'Extrusion_3' )
geompy.addToStudy( Sphere_, 'Sphere_' )
geompy.addToStudy( Multi_Translation_2_, 'Multi-Translation_2_' )
geompy.addToStudy( Multi_Translation_3_, 'Multi-Translation_3_' )
geompy.addToStudy( Pore1_, 'Pore1_' )
geompy.addToStudy( Pore2_, 'Pore2_' )
geompy.addToStudy( Pore3_, 'Pore3_' )
geompy.addToStudy( Cut_V_, 'Cut_V_' )
geompy.addToStudy( Pore4_, 'Pore4_' )
geompy.addToStudy( grains, 'grains' )
geompy.addToStudy( a3D_Sketcher_1_1, '3D Sketcher_1' )
geompy.addToStudy( Face_3, 'Face_3' )
geompy.addToStudy( Vector_Normal_1, 'Vector_Normal_1' )
geompy.addToStudy( cuban, 'cuban' )
geompy.addToStudyInFather( cuban, inlet, 'inlet' )
geompy.addToStudyInFather( cuban, outlet, 'outlet' )
geompy.addToStudyInFather( cuban, sp1, 'sp1' )
geompy.addToStudyInFather( cuban, sp2, 'sp2' )
geompy.addToStudyInFather( cuban, sp3, 'sp3' )
geompy.addToStudyInFather( cuban, sp4, 'sp4' )
geompy.addToStudy( Cut_1, 'Cut_1' )
geompy.addToStudy( Cut_2, 'Cut_2' )
geompy.addToStudy( Cut_3, 'Cut_3' )
geompy.addToStudy( Cut_4, 'Cut_4' )
geompy.addToStudy( Cut_5, 'Cut_5' )
geompy.addToStudy( Cut_6, 'Cut_6' )
geompy.addToStudy( Cut_7, 'Cut_7' )
geompy.addToStudyInFather( Cut_1, inlet_1, 'inlet' )
geompy.addToStudyInFather( Cut_1, outlet_1, 'outlet' )
if salome.sg.hasDesktop():
salome.sg.updateObjBrowser()

View File

@ -1,7 +0,0 @@
#!/usr/bin/env bash
for n in {0..3}; do
foamDictionary "processor${n}/0/U" -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary "processor${n}/0/U" -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
done

View File

@ -1,37 +0,0 @@
#!/usr/bin/env bash
path=.
#build/simple-cubic/0.1
# python src/genmesh.py
# python src/prefoam.py
# Clean case directory
rm -rf postProcessing processor* *.log logs *.obj constant/polyMesh
# Export and transform mesh
ideasUnvToFoam -case $path mesh.unv
transformPoints -scale '(1e-5 1e-5 1e-5)'
#polyDualMesh 70 -overwrite
checkMesh -case $path -allGeometry -allTopology > checkMesh.log
# Change boundary type for mesh
foamDictionary -case $path constant/polyMesh/boundary -entry entry0.wall.type -set wall
# Case decomposition
decomposePar -case $path
# Initial approximation
mpirun -np 4 --oversubscribe potentialFoam -case $path -parallel | tee -a $path/potentialFoam.log
# Change boundary type for simpleFoam
for n in {0..3}; do
foamDictionary "processor${n}/0/U" -entry boundaryField.inlet.type -set pressureInletVelocity
foamDictionary "processor${n}/0/U" -entry boundaryField.inlet.value -set 'uniform (0 0 0)'
done
sleep 2
# Main calculation
mpirun -np 4 --oversubscribe simpleFoam -parallel -case $path | tee -a $path/simpleFoam.log