From 9063fde50f472e13382de6993d7dfd56bca711c9 Mon Sep 17 00:00:00 2001 From: L-Nafaryus Date: Sat, 6 Mar 2021 19:19:43 +0500 Subject: [PATCH] Working 3rd direction --- TODO.md | 9 +- src/baseFOAM/system/createPatchDict | 76 ++++++++++ src/foam.py | 4 +- src/genmesh.py | 4 +- src/simpleCubic.py | 208 +++++++++++++++------------- 5 files changed, 199 insertions(+), 102 deletions(-) create mode 100644 src/baseFOAM/system/createPatchDict diff --git a/TODO.md b/TODO.md index a2ffe95..8a47711 100644 --- a/TODO.md +++ b/TODO.md @@ -15,6 +15,13 @@ - [x] 3rd direction - [ ] createPatch(Dict) - [ ] views (mesh, ..) -- [ ] alpha for simpleCubic [0.01 .. 0.28] +- [x] alpha for simpleCubic [0.01 .. 0.28] - [ ] translation vector (cyclicAMI) - [ ] 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 + diff --git a/src/baseFOAM/system/createPatchDict b/src/baseFOAM/system/createPatchDict new file mode 100644 index 0000000..f74d234 --- /dev/null +++ b/src/baseFOAM/system/createPatchDict @@ -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); + } +); + +// ************************************************************************* // diff --git a/src/foam.py b/src/foam.py index 1eba295..8bdc185 100644 --- a/src/foam.py +++ b/src/foam.py @@ -68,8 +68,8 @@ if __name__ == "__main__": # Main entry structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"] - directions = ["001", "100"] - coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ] + 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: diff --git a/src/genmesh.py b/src/genmesh.py index 3d31c5b..3f60d8b 100644 --- a/src/genmesh.py +++ b/src/genmesh.py @@ -60,8 +60,8 @@ if __name__ == "__main__": # Start in parallel processes = [] structures = ["simpleCubic"] #, "bc-cubic", "fc-cubic"] - directions = ["001", "100"] - coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ] + directions = ["001"] #, "100", "111"] + coefficients = [0.1] #[ alpha * 0.01 for alpha in range(1, 13 + 1) ] port = 2810 for structure in structures: diff --git a/src/simpleCubic.py b/src/simpleCubic.py index c2791ad..5d95761 100644 --- a/src/simpleCubic.py +++ b/src/simpleCubic.py @@ -20,6 +20,8 @@ class simpleCubic: self.rombus = None self.rombusbbox = None + self.spheres = None + salome.salome_init() def geometryCreate(self, alpha): @@ -29,8 +31,8 @@ class simpleCubic: Parameters: alpha (float): Sphere intersection parameter which used for cutting spheres from box. - Radius = R_0 / (1 - alpha) - Should be from 0.01 to 0.13 + Radius = R0 / (1 - alpha) + Should be from 0.01 to 0.28 Returns: Configured geometry. @@ -38,10 +40,21 @@ class simpleCubic: geompy = geomBuilder.New() - # - R_0 = 1 - R = R_0 / (1 - alpha) - R_fillet = 0 + # Parameters + R0 = 1 + R = R0 / (1 - alpha) + 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 axes = [ @@ -75,10 +88,17 @@ class simpleCubic: sphere3 = geompy.ExtractShapes(sphere3, geompy.ShapeType["SOLID"], True) sphere = geompy.MakeFuseList(sphere + sphere2 + sphere3, True, True) - + + if not R_fillet == 0: 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.geometrybbox = box @@ -86,28 +106,34 @@ class simpleCubic: # Rombus 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) - - 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) - self.rombusbbox = rombus + Vertex_2 = geompy.MakeVertex(0, 0, 4) + Vertex_1 = geompy.MakeVertex(2, 0, 2) + Vertex_3 = geompy.MakeVertex(2, 2, 0) + Vertex_4 = geompy.MakeVertex(0, 2, 2) + Edge_1 = geompy.MakeEdge(Vertex_2, Vertex_1) + 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) - Operators = ["FixShape"] - Parameters = ["FixShape.Tolerance3d"] - Values = ["1e-7"] - PS = geompy.ProcessShape(self.rombusbbox, Operators, Parameters, Values) - self.rombusbbox = PS + + #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) + + #a3D_Sketcher_1 = sk.wire() + #Face_1 = geompy.MakeFaceWires([a3D_Sketcher_1], 1) + Vector_1 = geompy.MakeVectorDXDYDZ(1, 1, 0) + rombusbbox = geompy.MakePrismVecH(Face_1, Vector_1, round(2 * math.sqrt(2), 14)) + + geompy.addToStudy(rombusbbox, "rombusbbox") + + self.rombus = geompy.MakeCutList(rombusbbox, [sphere], True) + self.rombusbbox = rombusbbox 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() rot = [0, 0, 45] buffergeometry = self.geometry @@ -174,88 +186,75 @@ class simpleCubic: elif direction == "111": center = geompy.MakeVertex(2, 2, 2) + self.geometry = self.rombus 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), math.sqrt(2) / 2)) + 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)) vstep = math.sqrt(2) hstep = 1 + logging.info("boundaryCreate: direction = {}".format(direction)) + 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": - box = self.rombusbbox + box = self.rombus else: box = self.geometrybbox planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True) - - inletplane = None - outletplane = None + inletplane = [] + outletplane = [] hplanes = [] - n = 0 + for plane in planes: planeNorm = geompy.GetNormal(plane) - n += 1 - geompy.addToStudy(planeNorm, "normalplane-{}".format(n)) - angle = abs(geompy.GetAngle(planeNorm, norm)) - logging.info("angle = {}".format(angle)) + + angle = round(abs(geompy.GetAngle(planeNorm, norm)), 0) if angle == 0: - outletplane = plane + outletplane.append(plane) 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) + if salome.sg.hasDesktop(): 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 - common1 = geompy.MakeCommonList([self.geometry, inletplane], True) - inlet = createGroup(common1, "inlet") + def createGroup(planelist, name): + gr = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], name) + + grcomp = geompy.MakeCompound(planelist) + grcut = geompy.MakeCutList(grcomp, [self.spheres], True) - common2 = geompy.MakeCommonList([self.geometry, outletplane], True) - outlet = createGroup(common2, "outlet") - - # symetryPlane(s) - 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) + gip = geompy.GetInPlace(self.geometry, grcut, True) 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 allgroup = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"]) @@ -299,6 +298,15 @@ class simpleCubic: """ 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) @@ -320,12 +328,18 @@ class simpleCubic: 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, name, SMESH.FACE) @@ -386,11 +400,11 @@ if __name__ == "__main__": sc.boundaryCreate(direction) logging.info("Creating the mesh ...") - sc.meshCreate(2, { - "thickness": 0.02, - "number": 2, - "stretch": 1.1 - }) + sc.meshCreate(2) #, { + # "thickness": 0.001, + # "number": 1, + # "stretch": 1.1 + #}) sc.meshCompute() logging.info("Exporting the mesh ...")