From 0ade0b0621f5ab1bdf7080fceecec918f54326f7 Mon Sep 17 00:00:00 2001 From: L-Nafaryus Date: Fri, 2 Apr 2021 21:41:48 +0500 Subject: [PATCH] Fix: mesh generation (not full) --- run.py | 28 ++++++++---- samples/__init__.py | 11 +++-- samples/bodyCenteredCubic.py | 2 +- src/geometry_utils.py | 25 ++++++++--- src/mesh_utils.py | 85 ++++++++++++++++++++++++++++++------ 5 files changed, 118 insertions(+), 33 deletions(-) diff --git a/run.py b/run.py index 7c4edad..144e213 100644 --- a/run.py +++ b/run.py @@ -33,7 +33,7 @@ def createTasks(): for structure in structures: if structure == "simpleCubic": - theta = [0.01, 0.29] #[c * 0.01 for c in range(1, 29 + 1)] + theta = [] #[0.01, 0.28] #[c * 0.01 for c in range(1, 29 + 1)] elif structure == "faceCenteredCubic": theta = [0.01, 0.13] #[c * 0.01 for c in range(1, 13 + 1)] @@ -67,7 +67,7 @@ def createMesh(tasks): returncode = salome_utils.runExecute(port, scriptpath, task.structure, task.coeff, "".join([str(coord) for coord in task.direction]), os.path.join(task.saveto, "mesh.unv")) logging.info("Return code:\t{}".format(returncode)) - logging.info("-" * 80) + #logging.info("-" * 80) if returncode == 1: break @@ -76,7 +76,7 @@ def createMesh(tasks): def calculate(tasks): foamCase = [ "0", "constant", "system" ] rmDirs = ["0", "constant", "system", "postProcessing", "logs"] + [ "processor{}".format(n) for n in range(4)] - fancyline = "--------------------------------------------------------------------------------" + #fancyline = "--------------------------------------------------------------------------------" for task in tasks: @@ -92,12 +92,13 @@ def calculate(tasks): os.chdir(task.saveto) casepath = "." - logging.info(fancyline) - logging.info("""Args: + logging.info("-" * 80) + logging.info("""calculate: + task:\t{} / {} structure type:\t{} coefficient:\t{} flow direction:\t{} - path:\t{}\n""".format(task.structure, task.coeff, task.direction, task.saveto)) + path:\t{}\n""".format(tasks.index(task) + 1, len(tasks) , task.structure, task.coeff, task.direction, task.saveto)) foam_utils.ideasUnvToFoam(casepath, "mesh.unv") @@ -136,7 +137,7 @@ def calculate(tasks): os.chdir(ROOT) - logging.info(fancyline) + #logging.info(fancyline) if __name__ == "__main__": @@ -178,14 +179,23 @@ if __name__ == "__main__": if args.mesh: start_time = time.monotonic() - logging.info("Started at {}".format(timedelta(seconds=start_time))) + #logging.info("Started at {}".format(timedelta(seconds=start_time))) createMesh(tasks) end_time = time.monotonic() + logging.info("-" * 80) logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time))) if args.calc: - calculate(tasks) + 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))) + diff --git a/samples/__init__.py b/samples/__init__.py index 17f782e..0678050 100644 --- a/samples/__init__.py +++ b/samples/__init__.py @@ -69,8 +69,13 @@ def genMesh(stype, theta, flowdirection, saveto): #boundary = geometry_utils.boundaryCreate(geometry, direction, grains) boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains) - fineness = 3 - mesh = mesh_utils.meshCreate(geometry, boundary, fineness) + fineness = 1 + viscousLayers = { + "thickness": 0.001, + "number": 2, + "stretch": 1 + } + mesh = mesh_utils.meshCreate(geometry, boundary, fineness, viscousLayers) mesh_utils.meshCompute(mesh) #path = os.path.join(saveto, @@ -118,7 +123,7 @@ if __name__ == "__main__": structure type:\t{} coefficient:\t{} flow direction:\t{} - export path:\t{}\n""".format(stype, theta, flowdirection, saveto)) + export path:\t{}""".format(stype, theta, flowdirection, saveto)) #print(flowdirection) diff --git a/samples/bodyCenteredCubic.py b/samples/bodyCenteredCubic.py index 4daa068..e348b0b 100644 --- a/samples/bodyCenteredCubic.py +++ b/samples/bodyCenteredCubic.py @@ -151,7 +151,7 @@ def bodyCenteredCubic(alpha): grains = geompy.MakeFuseList(grains, False, False) R = r0 / (1 - alpha) - C1 = 0.8 + C1 = 0.6 C2 = 0.4 alpha1 = 0.01 alpha2 = 0.13 diff --git a/src/geometry_utils.py b/src/geometry_utils.py index a0ed5e7..1336f9d 100644 --- a/src/geometry_utils.py +++ b/src/geometry_utils.py @@ -59,7 +59,9 @@ def createBoundary(gobj, bcount, dvec, norm, grains): side boundaries:\t{} normal direction:\t{} angles:\t{} - side directions:\t{}""".format(dvec, bcount, norm, [ ang(n, bcount) / (2 * np.pi) * 360 for n in range(limit) ], [ v for v in vecs ])) + side directions:\t{}""".format(dvec, bcount, norm, + [ ang(n, bcount) / (2 * np.pi) * 360 for n in range(limit) ], + len(vecs))) #[ v for v in vecs ])) # flowvec = geompy.MakeVector( @@ -76,7 +78,12 @@ def createBoundary(gobj, bcount, dvec, norm, grains): # planes = geompy.ExtractShapes(gobj, geompy.ShapeType["FACE"], False) - planes = geompy.MakeCompound(planes) + #planes = geompy.MakeCompound(planes) + planes = geompy.MakeShell(planes) + planes = geompy.ProcessShape(planes, + [ "FixShape", "FixFaceSize", "DropSmallEdges", "SameParameter" ], + [ "FixShape.Tolerance3d", "FixShape.MaxTolerance3d", "FixFaceSize.Tolerance", "DropSmallEdges.Tolerance3d", "SameParameter.Tolerance3d" ], + [ "1e-7", "1", "0.05", "0.05", "1e-7" ]) planes = geompy.MakeCutList(planes, [grains], False) planes = geompy.ExtractShapes(planes, geompy.ShapeType["FACE"], False) #print("planes: ", len(planes)) @@ -89,7 +96,7 @@ def createBoundary(gobj, bcount, dvec, norm, grains): nvec = geompy.GetNormal(plane) fwang = round(geompy.GetAngle(nvec, flowvec), 0) - print("fwang = ", fwang) + #print("fwang = ", fwang) if fwang == 0: inletplanes.append(plane) @@ -99,7 +106,7 @@ def createBoundary(gobj, bcount, dvec, norm, grains): for n in range(len(symvec)): sang = round(geompy.GetAngle(nvec, symvec[n]), 0) - print("sang = ", sang) + #print("sang = ", sang) if sang == 0: if symetryplanes[n][0] == None: @@ -112,13 +119,19 @@ def createBoundary(gobj, bcount, dvec, norm, grains): symetryplanes[n][1] = [] symetryplanes[n][1].append(plane) - print("\n") + #print("\n") + + symetryplanesinfo = [] + for n in range(len(symetryplanes)): + symetryplanesinfo.append([]) + for pair in range(len(symetryplanes[n])): + symetryplanesinfo[n].append(len(symetryplanes[n][pair])) logging.info("""createBoundary: planes:\t{} inlet planes:\t{} outlet planes:\t{} - side planes:\t{}""".format(len(planes), inletplanes, outletplanes, symetryplanes)) + side planes:\t{}""".format(len(planes), len(inletplanes), len(outletplanes), symetryplanesinfo)) # boundary = {} diff --git a/src/mesh_utils.py b/src/mesh_utils.py index 6d196ce..b8d6555 100644 --- a/src/mesh_utils.py +++ b/src/mesh_utils.py @@ -39,36 +39,61 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None): 1: "Coarse", 2: "Moderate", 3: "Fine", - 4: "Very fine" + 4: "Very fine", + 5: "Custom" }[fineness] - logging.info("""meshCreate: - mesh fineness:\t{}""".format(Fineness)) - mesh = smesh.Mesh(gobj) netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D) param = netgen.Parameters() + param.SetMinSize( 0.001 ) + param.SetMaxSize( 0.1 ) + param.SetSecondOrder( 0 ) param.SetOptimize( 1 ) + param.SetQuadAllowed( 0 ) 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 ) + + # TODO: add customization + if fineness == 5: + param.SetGrowthRate( 0.1 ) + param.SetNbSegPerEdge( 5 ) + param.SetNbSegPerRadius( 10 ) + + + logging.info("""meshCreate: + fineness:\t{} + min size:\t{} + max size:\t{} + growth rate:\t{} + nb segs per edge:\t{} + nb segs per radius:\t{} + limit size by surface curvature:\t{} + quad-dominated:\t{} + second order:\t{} + optimize:\t{}""".format( + Fineness, param.GetMinSize(), param.GetMaxSize(), + param.GetGrowthRate(), param.GetNbSegPerEdge(), param.GetNbSegPerRadius(), + True if param.GetUseSurfaceCurvature() else False, + True if param.GetQuadAllowed() else False, + True if param.GetSecondOrder() else False, + True if param.GetOptimize() else False)) + + if not viscousLayers is None: logging.info("""meshCreate: - viscous layers: thickness = {} - number = {} - stretch factor = {}""".format( + viscous layers: + thickness:\t{} + number:\t{} + stretch factor:\t{}""".format( viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"])) vlayer = netgen.ViscousLayers( @@ -80,7 +105,7 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None): else: logging.info("""meshCreate: - viscous layers: disabled""") + viscous layers: disabled""") for name, b in boundary.items(): mesh.GroupOnGeom(b, "{}_".format(name), SMESH.FACE) @@ -102,6 +127,38 @@ def meshCompute(mobj): logging.info("""meshCompute: status:\t{}""".format(msg)) + if status: + omniinfo = mobj.GetMeshInfo() + keys = [ str(k) for k in omniinfo.keys() ] + vals = [ v for v in omniinfo.values() ] + info = {} + + for n in range(len(keys)): + info[keys[n]] = vals[n] + + edges = info["Entity_Edge"] + + triangles = info["Entity_Triangle"] + faces = triangles + + tetra = info["Entity_Tetra"] + prism = info["Entity_Penta"] + pyramid = info["Entity_Pyramid"] + volumes = tetra + prism + pyramid + + elements = edges + faces + volumes + + logging.info("""meshCompute: + Elements:\t{} + Edges:\t{} + Faces:\t{} + Triangles:\t{} + Volumes:\t{} + Tetrahedrons:\t{} + Prisms:\t{} + Pyramid:\t{}""".format( + elements, edges, faces, triangles, volumes, tetra, prism, pyramid)) + def meshExport(mobj, path): """