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

Working 3rd direction

Browse Source
This commit is contained in:
L-Nafaryus 2021-03-06 19:19:43 +05:00
parent 2586acfa54
commit 9063fde50f
No known key found for this signature in database
GPG Key ID: C76D8DCD2727DBB7
5 changed files with 199 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
TODO.md
View File

@ -15,6 +15,13 @@
- [x] 3rd direction - [x] 3rd direction
- [ ] createPatch(Dict) - [ ] createPatch(Dict)
- [ ] views (mesh, ..) - [ ] views (mesh, ..)
- [ ] alpha for simpleCubic [0.01 .. 0.28] - [x] alpha for simpleCubic [0.01 .. 0.28]
- [ ] translation vector (cyclicAMI) - [ ] translation vector (cyclicAMI)
- [ ] BUG: angle between the direction vector and the normal to inlet is ~1.4e-14 - [ ] BUG: angle between the direction vector and the normal to inlet is ~1.4e-14
- [x] Temporary solution
- [ ] BUG: ideasUnvToFoam not working with param '-case PATH'
- [x] Temporary sulution via os.chdir(PATH)
## 6.03.21
- [ ] ERROR: MakeFuseList with alpha > 0.2

76
src/baseFOAM/system/createPatchDict Normal file
View File

@ -0,0 +1,76 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object createPatchDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
pointSync false;
// Patches to create.
patches
(
{
name inlet;
patchInfo
{
type patch;
inGroups (inlet);
}
constructFrom patches;
patches (inlet);
}
{
name outlet;
patchInfo
{
type patch;
inGroups (ouetlet);
}
constructFrom patches;
patches (outlet);
}
{
name symetryPlane;
patchInfo
{
type symetryPlane;
inGroups (symetryPlane);
}
constructFrom patches;
patches (symetryPlane);
}
{
name wall;
patchInfo
{
type wall;
inGroups (wall);
}
constructFrom patches;
patches (wall);
}
);
// ************************************************************************* //

4
src/foam.py
View File

@ -68,8 +68,8 @@ if __name__ == "__main__":
# Main entry # Main entry
structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"] structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"]
directions = ["001", "100"] directions = ["001"] #, "100", "111"]
coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ] coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ]
for structure in structures: for structure in structures:
for direction in directions: for direction in directions:

4
src/genmesh.py
View File

@ -60,8 +60,8 @@ if __name__ == "__main__":
# Start in parallel # Start in parallel
processes = [] processes = []
structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"] structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"]
directions = ["001", "100"] directions = ["001"] #, "100", "111"]
coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ] coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ]
port = 2810 port = 2810
for structure in structures: for structure in structures:

200
src/simpleCubic.py
View File

@ -20,6 +20,8 @@ class simpleCubic:
self.rombus = None self.rombus = None
self.rombusbbox = None self.rombusbbox = None
self.spheres = None
salome.salome_init() salome.salome_init()
def geometryCreate(self, alpha): def geometryCreate(self, alpha):
@ -29,8 +31,8 @@ class simpleCubic:
Parameters: Parameters:
alpha (float): Sphere intersection parameter which used for cutting spheres from box. alpha (float): Sphere intersection parameter which used for cutting spheres from box.
Radius = R_0 / (1 - alpha) Radius = R0 / (1 - alpha)
Should be from 0.01 to 0.13 Should be from 0.01 to 0.28
Returns: Returns:
Configured geometry. Configured geometry.
@ -38,10 +40,21 @@ class simpleCubic:
geompy = geomBuilder.New() geompy = geomBuilder.New()
# # Parameters
R_0 = 1 R0 = 1
R = R_0 / (1 - alpha) R = R0 / (1 - alpha)
R_fillet = 0 R_fillet = 0.01
C1 = 0.8
C2 = 0.4
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))
# xyz axes # xyz axes
axes = [ axes = [
@ -76,9 +89,16 @@ class simpleCubic:
sphere = geompy.MakeFuseList(sphere + sphere2 + sphere3, True, True) sphere = geompy.MakeFuseList(sphere + sphere2 + sphere3, True, True)
if not R_fillet == 0: if not R_fillet == 0:
sphere = geompy.MakeFilletAll(sphere, R_fillet) sphere = geompy.MakeFilletAll(sphere, R_fillet)
self.spheres = sphere
#else:
# sphere = sphere + sphere2 + sphere3 #geompy.MakeCompound(sphere + sphere2 + sphere3)
# geompy.RemoveExtraEdges(obj, True)
self.geometry = geompy.MakeCutList(box, [sphere], True) self.geometry = geompy.MakeCutList(box, [sphere], True)
self.geometrybbox = box self.geometrybbox = box
@ -87,27 +107,33 @@ class simpleCubic:
# Rombus # Rombus
h = 2 h = 2
sk = geompy.Sketcher3D() Vertex_2 = geompy.MakeVertex(0, 0, 4)
sk.addPointsAbsolute(0, 0, h * 2) Vertex_1 = geompy.MakeVertex(2, 0, 2)
sk.addPointsAbsolute(h, 0, h) Vertex_3 = geompy.MakeVertex(2, 2, 0)
sk.addPointsAbsolute(h, h, 0) Vertex_4 = geompy.MakeVertex(0, 2, 2)
sk.addPointsAbsolute(0, h, h) Edge_1 = geompy.MakeEdge(Vertex_2, Vertex_1)
sk.addPointsAbsolute(0, 0, h * 2) Edge_2 = geompy.MakeEdge(Vertex_1, Vertex_3)
Edge_3 = geompy.MakeEdge(Vertex_3, Vertex_4)
Edge_4 = geompy.MakeEdge(Vertex_4, Vertex_2)
Face_1 = geompy.MakeFaceWires([Edge_1, Edge_2, Edge_3, Edge_4], 1)
a3D_Sketcher_1 = sk.wire()
Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Vector_1 = geompy.MakeVectorDXDYDZ(h, h, 0)
rombus = geompy.MakePrismVecH(Face_1, Vector_1, 2 * math.sqrt(2))
geompy.addToStudy(rombus, "romb")
self.rombus = geompy.MakeCutList(rombus, [sphere], True) #sk = geompy.Sketcher3D()
self.rombusbbox = rombus #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)
Operators = ["FixShape"] #a3D_Sketcher_1 = sk.wire()
Parameters = ["FixShape.Tolerance3d"] #Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1)
Values = ["1e-7"] Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 0)
PS = geompy.ProcessShape(self.rombusbbox, Operators, Parameters, Values) rombusbbox = geompy.MakePrismVecH(Face_1, Vector_1, round(2 * math.sqrt(2), 14))
self.rombusbbox = PS
geompy.addToStudy(rombusbbox, "rombusbbox")
self.rombus = geompy.MakeCutList(rombusbbox, [sphere], True)
self.rombusbbox = rombusbbox
geompy.addToStudy(self.rombus, "rombus") geompy.addToStudy(self.rombus, "rombus")
@ -134,20 +160,6 @@ class simpleCubic:
} }
""" """
#
# _____ z |
# //////| | | flow
# ////// | |___y f
# | | / /
# |____|/ /x direction [0, 0, 1]
#
# _____ z f
# / /| | / flow
# /____/ | |___y /
# |||||| / /
# ||||||/ /x direction [1, 0, 0]
#
geompy = geomBuilder.New() geompy = geomBuilder.New()
rot = [0, 0, 45] rot = [0, 0, 45]
buffergeometry = self.geometry buffergeometry = self.geometry
@ -174,88 +186,75 @@ class simpleCubic:
elif direction == "111": elif direction == "111":
center = geompy.MakeVertex(2, 2, 2) center = geompy.MakeVertex(2, 2, 2)
self.geometry = self.rombus
norm = geompy.MakeVector(center, norm = geompy.MakeVector(center,
geompy.MakeVertexWithRef(center, geompy.MakeVertexWithRef(center, 1, 1, 1))
-math.cos((90 + rot[2]) * math.pi / 180.0), #-math.cos((90 + rot[2]) * math.pi / 180.0),
math.sin((90 + rot[2]) * math.pi / 180.0), math.sqrt(2) / 2)) #math.sin((90 + rot[2]) * math.pi / 180.0), math.sqrt(2) / 2))
vstep = math.sqrt(2) vstep = math.sqrt(2)
hstep = 1 hstep = 1
logging.info("boundaryCreate: direction = {}".format(direction))
geompy.addToStudy(norm, "normalvector") geompy.addToStudy(norm, "normalvector")
def createGroup(shape, name):
self.geometry = self.rombus
group = geompy.CreateGroup(self.geometry,
geompy.ShapeType["FACE"], name)
gip = geompy.GetInPlace(self.geometry, shape, True)
faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
geompy.UnionList(group, faces)
return group
# xyz axes
#axes = [
# geompy.MakeVectorDXDYDZ(1, 0, 0),
# geompy.MakeVectorDXDYDZ(0, 1, 0),
# geompy.MakeVectorDXDYDZ(0, 0, 1)
#]
# Bounding box
#box = geompy.MakeBoxDXDYDZ(2 * math.sqrt(2), 2 * math.sqrt(2), 2)
#box = geompy.MakeRotation(box, axes[2], 45 * math.pi / 180.0)
#box = geompy.MakeTranslation(box, 2, 0, 0)
if direction == "111": if direction == "111":
box = self.rombusbbox box = self.rombus
else: else:
box = self.geometrybbox box = self.geometrybbox
planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True) planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True)
inletplane = []
inletplane = None outletplane = []
outletplane = None
hplanes = [] hplanes = []
n = 0
for plane in planes: for plane in planes:
planeNorm = geompy.GetNormal(plane) planeNorm = geompy.GetNormal(plane)
n += 1
geompy.addToStudy(planeNorm, "normalplane-{}".format(n)) angle = round(abs(geompy.GetAngle(planeNorm, norm)), 0)
angle = abs(geompy.GetAngle(planeNorm, norm))
logging.info("angle = {}".format(angle))
if angle == 0: if angle == 0:
outletplane = plane outletplane.append(plane)
elif angle == 180: elif angle == 180:
inletplane = plane inletplane.append(plane)
else: elif direction == "111" and (angle == 109 or angle == 71):
hplanes.append(plane)
elif direction == "100" or direction == "001":
if angle == 90:
hplanes.append(plane) hplanes.append(plane)
if salome.sg.hasDesktop(): if salome.sg.hasDesktop():
salome.sg.updateObjBrowser() salome.sg.updateObjBrowser()
logging.info("hplanes = {}".format(len(hplanes))) logging.info("boundaryCreate: inletplanes = {}, outletplanes = {}, hplanes = {}".format(
len(inletplane), len(outletplane), len(hplanes)))
# inlet and outlet def createGroup(planelist, name):
common1 = geompy.MakeCommonList([self.geometry, inletplane], True) gr = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], name)
inlet = createGroup(common1, "inlet")
common2 = geompy.MakeCommonList([self.geometry, outletplane], True) grcomp = geompy.MakeCompound(planelist)
outlet = createGroup(common2, "outlet") grcut = geompy.MakeCutList(grcomp, [self.spheres], True)
# symetryPlane(s) gip = geompy.GetInPlace(self.geometry, grcut, True)
symetryPlane = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], "symetryPlane")
for plane in hplanes:
common3 = geompy.MakeCommonList([self.geometry, plane], True)
gip = geompy.GetInPlace(self.geometry, common3, True)
faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"]) faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"])
geompy.UnionList(symetryPlane, faces) geompy.UnionList(gr, faces)
return gr
# Main groups
inlet = createGroup(inletplane, "inlet")
outlet = createGroup(outletplane, "outlet")
symetryPlane = createGroup(hplanes, "symetryPlane")
# wall # wall
allgroup = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"]) allgroup = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"])
@ -299,6 +298,15 @@ class simpleCubic:
""" """
smesh = smeshBuilder.New() 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) mesh = smesh.Mesh(self.geometry)
netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D) netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
@ -320,12 +328,18 @@ class simpleCubic:
param.SetQuadAllowed( 0 ) param.SetQuadAllowed( 0 )
if not viscousLayers is None: 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"], vlayer = netgen.ViscousLayers(viscousLayers["thickness"],
viscousLayers["number"], viscousLayers["number"],
viscousLayers["stretch"], viscousLayers["stretch"],
[self.boundary["inlet"], self.boundary["outlet"]], [self.boundary["inlet"], self.boundary["outlet"]],
1, smeshBuilder.NODE_OFFSET) 1, smeshBuilder.NODE_OFFSET)
else:
logging.info("meshCreate: viscous layers are disabled")
for name, boundary in self.boundary.items(): for name, boundary in self.boundary.items():
mesh.GroupOnGeom(boundary, name, SMESH.FACE) mesh.GroupOnGeom(boundary, name, SMESH.FACE)
@ -386,11 +400,11 @@ if __name__ == "__main__":
sc.boundaryCreate(direction) sc.boundaryCreate(direction)
logging.info("Creating the mesh ...") logging.info("Creating the mesh ...")
sc.meshCreate(2, { sc.meshCreate(2) #, {
"thickness": 0.02, # "thickness": 0.001,
"number": 2, # "number": 1,
"stretch": 1.1 # "stretch": 1.1
}) #})
sc.meshCompute() sc.meshCompute()
logging.info("Exporting the mesh ...") logging.info("Exporting the mesh ...")