New: remove trash

This commit is contained in:
L-Nafaryus 2021-05-26 12:12:37 +05:00
parent 5ac5231574
commit 08de009cc4
No known key found for this signature in database
GPG Key ID: C76D8DCD2727DBB7
25 changed files with 0 additions and 9242 deletions

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