From f62ce3d5b422edcf490aba3d6e96453bae44f74a Mon Sep 17 00:00:00 2001 From: L-Nafaryus Date: Wed, 3 Mar 2021 13:00:28 +0500 Subject: [PATCH] 2 directions for boundary faces --- src/genmesh.py | 74 +++++++---------- src/simple-cubic/geometry.py | 83 ------------------- src/simple-cubic/main.py | 35 -------- src/simple-cubic/mesh.py | 34 -------- src/simpleCubic.py | 156 ++++++++++++++++++++++++++--------- src/structure.gv | 14 ---- 6 files changed, 149 insertions(+), 247 deletions(-) delete mode 100644 src/simple-cubic/geometry.py delete mode 100644 src/simple-cubic/main.py delete mode 100644 src/simple-cubic/mesh.py delete mode 100644 src/structure.gv diff --git a/src/genmesh.py b/src/genmesh.py index 7af546c..2a53135 100644 --- a/src/genmesh.py +++ b/src/genmesh.py @@ -2,52 +2,40 @@ import os import subprocess import multiprocessing -src = os.getcwd() -build = os.path.join(src, "../build") +def salome(src_path, build_path, coefficient, direction): + subprocess.run(["salome", "start", "-t", src_path, "args:{},{}".format(build_path, coefficient, direction)]) -if not os.path.exists(build): - os.makedirs(build) +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) -alpha = [] + ### + processes = [] + structures = ["simple-cubic"] #, "bc-cubic", "fc-cubic"] + directions = ["001", "100"] + coefficients = [ alpha * 0.01 for alpha in range(1, 13 + 1) ] -for n in range(0.01, 0.13 + 0.01, 0.01): - alpha.append(n) + 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) -# Structures -simpleCubic = os.path.join(src, "simple-cubic/main.py") -# Body-centered cubic -#bcCubic = os.path.join(path, "bc-cubic/main.py") -# Face-centered cubic -#fcCubic = os.path.join(path, "fc-cubic/main.py") + print("starting process") + p = multiprocessing.Process(target = salome, args = (src_path, build_path, coefficient, direction)) + processes.append(p) + p.start() -### - -processes = [] -structure = ["simple-cubic"] #, "bc-cubic", "fc-cubic"] -direction = ["1"] - -def salome(src_path, build_path, arg): - subprocess.run(["salome", "start", "-t", src_path, "args:{},{}".format(build_path, arg)]) - -for s in structure: - s_path = os.path.join(build, s) - - for d in direction: - d_path = os.path.join(s_path, d) - - for c in alpha: - src_path = os.path.join(src, "%s/main.py" % s) - build_path = os.path.join(d_path, str(c)) - - if not os.path.exists(build_path): - os.makedirs(build_path) - - print("starting process") - p = multiprocessing.Process(target = salome, args = (src_path, build_path, c)) - processes.append(p) - p.start() - - for process in processes: - process.join() + for process in processes: + process.join() diff --git a/src/simple-cubic/geometry.py b/src/simple-cubic/geometry.py deleted file mode 100644 index 66241c1..0000000 --- a/src/simple-cubic/geometry.py +++ /dev/null @@ -1,83 +0,0 @@ -import GEOM -from salome.geom import geomBuilder -geompy = geomBuilder.New() - -import math - -def create(alpha): - # xyz axes - axes = [ - geompy.MakeVectorDXDYDZ(1, 0, 0), - geompy.MakeVectorDXDYDZ(0, 1, 0), - geompy.MakeVectorDXDYDZ(0, 0, 1) - ] - - # Main 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) - - vtx = [ - geompy.MakeVertex(2, 0, 0), - geompy.MakeVertex(2, 2, 0), - geompy.MakeVertex(2, 2, 2) - ] - - line = geompy.MakeLineTwoPnt(vtx[1], vtx[2]) - - # Spheres for cutting - sphere = geompy.MakeSpherePntR(vtx[0], 1 / (1 - alpha)) - sphere = geompy.MakeMultiTranslation1D(sphere, axes[1], 2, 3) - cut = geompy.MakeCutList(box, [sphere], True) - - sphere2 = geompy.MakeTranslation(sphere, 0, 0, 2) - cut2 = geompy.MakeCutList(cut, [sphere2], True) - - sphere3 = geompy.MakeRotation(sphere, line, 90 * math.pi / 180.0) - cut3 = geompy.MakeCutList(cut2, [sphere3], True) - - sphere4 = geompy.MakeTranslation(sphere3, 0, 0, 2) - - # Main geometry - Pore = geompy.MakeCutList(cut3, [sphere4], True) - geompy.addToStudy(Pore, "PORE") - - # - # D /\ B A(1, 1, 0) \vec{n}(1, 1, 0) - # / \ B(3, 3, 0) \vec{n}(1, 1, 0) - # \ / C(3, 1, 0) \vec{n}(-1, 1, 0) - # A \/ C D(1, 3, 0) \vec{n}(-1, 1, 0) - # - - # Prepare faces - vtx.append(geompy.MakeVertex(2, 2, 2)) - vtx.append(geompy.MakeVertexWithRef(vtx[3], 0, 0, 1)) - vec2 = geompy.MakeVector(vtx[3], vtx[4]) - plane = geompy.MakePlane(vtx[3], vec2, 5) - common1 = geompy.MakeCommonList([Pore, plane], True) - plane = geompy.MakeTranslation(plane, 0, 0, -2) - common2 = geompy.MakeCommonList([Pore, plane], True) - - # TODO: make 4 planes A B C D - Pore = geompy.MakeFillet(Pore, 0.1, geompy.ShapeType["EDGE"], [24, 27, 35]) - - - # Main groups (inlet, outlet, wall) - inlet = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"]) - gip = geompy.GetInPlace(Pore, common1, True) - faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"]) - geompy.UnionList(inlet, faces) - - outlet = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"]) - gip = geompy.GetInPlace(Pore, common2, True) - faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"]) - geompy.UnionList(outlet, faces) - - PoreGroup = geompy.CreateGroup(Pore, geompy.ShapeType["FACE"]) - faces = geompy.SubShapeAllIDs(Pore, geompy.ShapeType["FACE"]) - geompy.UnionIDs(PoreGroup, faces) - - wall = geompy.CutListOfGroups([PoreGroup], [inlet, outlet]) - - - return Pore, [inlet, outlet, wall] diff --git a/src/simple-cubic/main.py b/src/simple-cubic/main.py deleted file mode 100644 index 4d29336..0000000 --- a/src/simple-cubic/main.py +++ /dev/null @@ -1,35 +0,0 @@ -import os, sys - -import salome -salome.salome_init() - -import geometry, mesh - -#alpha = [ 0.1, 0.15, 0.2 ] -build_path = str(sys.argv[1]) -coef = float(sys.argv[2]) - -#for coef in alpha: -print("alpha = {}".format(coef)) -print("Building geometry ...") - -Pore, bc = geometry.create(coef) -#geometry.geompy.addToStudy(Pore, 'Pore {}'.format(coef)) -#geometry.geompy.addToStudyInFather(Pore, bc[0], "INLET") -print("Building mesh ...") - -PoreMesh = mesh.create(Pore, bc) -isDone = PoreMesh.Compute() - -status = "Succesfully" if isDone else "Mesh is not computed" -print(status) - -try: - filename = os.path.join(build_path, "mesh.unv") - PoreMesh.ExportUNV(filename) - pass -except: - print('ExportUNV() failed. Invalid file name?') - -if salome.sg.hasDesktop(): - salome.sg.updateObjBrowser() diff --git a/src/simple-cubic/mesh.py b/src/simple-cubic/mesh.py deleted file mode 100644 index 9049e0c..0000000 --- a/src/simple-cubic/mesh.py +++ /dev/null @@ -1,34 +0,0 @@ -import SMESH, SALOMEDS -from salome.smesh import smeshBuilder -smesh = smeshBuilder.New() - - -def create(geomObj, bc): - mesh = smesh.Mesh(geomObj) - 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( 4 ) - #param.SetGrowthRate( 0.1 ) - #param.SetNbSegPerEdge( 5 ) - #param.SetNbSegPerRadius( 10 ) - param.SetQuadAllowed( 0 ) - - vlayer = netgen.ViscousLayers(0.025, 2, 1.1, [], - 1, smeshBuilder.NODE_OFFSET) - - mesh.GroupOnGeom(bc[0], 'inlet', SMESH.FACE) - mesh.GroupOnGeom(bc[1], 'outlet', SMESH.FACE) - mesh.GroupOnGeom(bc[2], 'wall', SMESH.FACE) - - return mesh - diff --git a/src/simpleCubic.py b/src/simpleCubic.py index bb45e03..2954687 100644 --- a/src/simpleCubic.py +++ b/src/simpleCubic.py @@ -15,6 +15,19 @@ class simpleCubic: salome.salome_init() def geometryCreate(self, alpha): + """ + Create the simple cubic geometry. + + 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 + + Returns: + Configured geometry. + """ + geompy = geomBuilder.New() # @@ -58,15 +71,44 @@ class simpleCubic: geompy.addToStudy(self.geometry, self.name) + return self.geometry + def boundaryCreate(self, direction): - geompy = geomBuilder.New() + """ + 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. + + Returns: + boundary (dict): + + { + "inlet": , + "outlet": , + "symetryPlane": , + "wall": + } + + """ # - # D /\ B A(1, 1, 0) \vec{n}(1, 1, 0) - # / \ B(3, 3, 0) \vec{n}(1, 1, 0) - # \ / C(3, 1, 0) \vec{n}(-1, 1, 0) - # A \/ C D(1, 3, 0) \vec{n}(-1, 1, 0) + # _____ z | + # //////| | | flow + # ////// | |___y f + # | | / / + # |____|/ /x direction [0, 0, 1] + # + # _____ z f + # / /| | / flow + # /____/ | |___y / + # |||||| / / + # ||||||/ /x direction [1, 0, 0] # + geompy = geomBuilder.New() center = geompy.MakeVertex(2, 2, 1) rot = [0, 0, 45] @@ -78,7 +120,9 @@ class simpleCubic: elif direction == "100": norm = geompy.MakeVector(center, - geompy.MakeVertexWithRef(center, 1, 0, 0)) + geompy.MakeVertexWithRef(center, + math.cos((90 + rot[2]) * math.pi / 180.0), + -math.sin((90 + rot[2]) * math.pi / 180.0), 0)) vstep = math.sqrt(2) hstep = 1 @@ -102,14 +146,13 @@ class simpleCubic: 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) - planes = geompy.ExtractShapes(box, - geompy.ShapeType["FACE"], True) + planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True) vplanes = [] hplanes = [] for plane in planes: planeNorm = geompy.GetNormal(plane) - angle = geompy.GetAngle(planeNorm, norm) + angle = abs(geompy.GetAngle(planeNorm, norm)) if angle == 0 or angle == 180: vplanes.append(plane) @@ -117,26 +160,29 @@ class simpleCubic: else: hplanes.append(plane) - zvplane1 = vplanes[0] - zvplane2 = vplanes[1] - if direction == "001": - if geompy.GetPosition(zvplane1)[3] > geompy.GetPosition(zvplane1)[3]: - inletplane = zvplane1 - outletplane = zvplane2 + z1 = geompy.GetPosition(vplanes[0])[3] + z2 = geompy.GetPosition(vplanes[1])[3] + + if z1 > z2: + inletplane = vplanes[0] + outletplane = vplanes[1] else: - inletplane = zvplane2 - outletplane = zvplane1 + inletplane = vplanes[1] + outletplane = vplanes[0] elif direction == "100": - if geompy.GetPosition(zvplane1)[1] > geompy.GetPosition(zvplane1)[1]: - inletplane = zvplane1 - outletplane = zvplane2 + x1 = geompy.GetPosition(vplanes[0])[1] + x2 = geompy.GetPosition(vplanes[1])[1] + + if x1 > x2: + inletplane = vplanes[0] + outletplane = vplanes[1] else: - inletplane = zvplane2 - outletplane = zvplane1 + inletplane = vplanes[1] + outletplane = vplanes[0] # inlet and outlet @@ -146,9 +192,6 @@ class simpleCubic: common2 = geompy.MakeCommonList([self.geometry, outletplane], True) outlet = createGroup(common2, "outlet") - geompy.addToStudy(inletplane, "inletplane") - geompy.addToStudy(outletplane, "outletplane") - # symetryPlane(s) symetryPlane = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"], "symetryPlane") @@ -173,10 +216,35 @@ class simpleCubic: return self.boundary - def meshCreate(self): + 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 , containig the parameters and boundary groups. + + """ smesh = smeshBuilder.New() - mesh = smesh.Mesh(geomObj) + mesh = smesh.Mesh(self.geometry) netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D) param = netgen.Parameters() @@ -189,25 +257,28 @@ class simpleCubic: param.SetCheckChartBoundary( 0 ) param.SetMinSize( 0.01 ) param.SetMaxSize( 0.1 ) - param.SetFineness( 4 ) + param.SetFineness(fineness) #param.SetGrowthRate( 0.1 ) #param.SetNbSegPerEdge( 5 ) #param.SetNbSegPerRadius( 10 ) param.SetQuadAllowed( 0 ) - - vlayer = netgen.ViscousLayers(0.025, 2, 1.1, [], - 1, smeshBuilder.NODE_OFFSET) - mesh.GroupOnGeom(self.boundary["inlet"], "inlet", SMESH.FACE) - mesh.GroupOnGeom(self.boundary["outlet"], "outlet", SMESH.FACE) - mesh.GroupOnGeom(self.boundary["symetryPlane"], "symetryPlane", SMESH.FACE) - mesh.GroupOnGeom(self.boundary["wall"], "wall", SMESH.FACE) + if not viscousLayers is None: + vlayer = netgen.ViscousLayers(viscousLayers["thickness"], + viscousLayers["number"], + viscousLayers["stretch"], + [self.boundary["inlet"], self.boundary["outlet"]], + 1, smeshBuilder.NODE_OFFSET) + + for name, boundary in self.boundary.items(): + mesh.GroupOnGeom(boundary, name, SMESH.FACE) self.mesh = mesh return self.mesh def meshCompute(self): + """Compute the mesh.""" status = self.mesh.Compute() if status: @@ -217,6 +288,12 @@ class simpleCubic: print("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: @@ -230,12 +307,15 @@ if __name__ == "__main__": alpha = float(sys.argv[2]) direction = str(sys.argv[3]) - #salome.salome_init() sc = simpleCubic("simpleCubic-{}-{}".format(direction, alpha)) sc.geometryCreate(alpha) sc.boundaryCreate(direction) - #sc.meshCreate() - #sc.meshCompute() + sc.meshCreate(2, { + "thickness": 0.02, + "number": 2, + "stretch": 1.1 + }) + sc.meshCompute() #sc.meshExport(build) if salome.sg.hasDesktop(): diff --git a/src/structure.gv b/src/structure.gv deleted file mode 100644 index 9e56648..0000000 --- a/src/structure.gv +++ /dev/null @@ -1,14 +0,0 @@ -digraph structure { - node [shape=plaintext] - s1 [label="src/simple-cubic/geometry.py"] - s2 [label="src/simple-cubic/mesh.py"] - s3 [label="src/simple-cubic/main.py"] - s1 -> s3 - s2 -> s3 - s4 [label="src/genmesh.py"] - s3 -> s4 - b1 [label="build/simple-cubic/direction-1/0.1/"] - s4 -> b1 [label="[alpha=0.1, direction=1, mesh.unv]"] - s5 [label="src/prefoam.py"] - s5 -> b1 [label="[0, constant, system (/src/foam/*)]"] -}