From fc81b84c8cf522841a62a86a44cbac2b3424d3fb Mon Sep 17 00:00:00 2001 From: L-Nafaryus Date: Mon, 12 Apr 2021 16:00:26 +0500 Subject: [PATCH] Destruction [4]: Mesh --- run.py | 22 +- samples/__init__.py | 132 ++++---- samples/bodyCentered.py | 284 ++++++++++++++++++ samples/faceCentered.py | 282 +++++++++++++++++ samples/{test.py => simple.py} | 48 ++- src/mesh_utils.py | 111 ++++--- temp/samples.old.2/__init__.py | 100 ++++++ .../samples.old.2}/bodyCenteredCubic.py | 0 .../samples.old.2}/faceCenteredCubic.py | 0 .../samples.old.2}/simpleCubic.py | 0 10 files changed, 866 insertions(+), 113 deletions(-) create mode 100644 samples/bodyCentered.py create mode 100644 samples/faceCentered.py rename samples/{test.py => simple.py} (85%) create mode 100644 temp/samples.old.2/__init__.py rename {samples => temp/samples.old.2}/bodyCenteredCubic.py (100%) rename {samples => temp/samples.old.2}/faceCenteredCubic.py (100%) rename {samples => temp/samples.old.2}/simpleCubic.py (100%) diff --git a/run.py b/run.py index 85f751d..52d71e5 100644 --- a/run.py +++ b/run.py @@ -17,10 +17,13 @@ from src import salome_utils from src import foam_utils def createTasks(): + ### + # Control variables + ## structures = [ - "simpleCubic", - "bodyCenteredCubic", - "faceCenteredCubic" + "simple", + #"bodyCentered", + #"faceCentered" ] directions = [ [1, 0, 0], @@ -29,19 +32,22 @@ def createTasks(): ] fillet = 1 + ### + # Tasks + ## Task = namedtuple("Task", ["structure", "coeff", "fillet", "direction", "saveto"]) tasks = [] for structure in structures: - if structure == "simpleCubic": - #theta = [c * 0.01 for c in range(1, 28 + 1)] - theta = [ 0.01, 0.28 ] + if structure == "simple": + theta = [c * 0.01 for c in range(1, 28 + 1)] + #theta = [ 0.01, 0.28 ] - elif structure == "faceCenteredCubic": + elif structure == "faceCentered": #theta = [c * 0.01 for c in range(1, 13 + 1)] theta = [ 0.01, 0.13 ] - elif structure == "bodyCenteredCubic": + elif structure == "bodyCentered": #theta = [c * 0.01 for c in range(1, 18 + 1)] theta = [ 0.01, 0.13, 0.14, 0.18 ] diff --git a/samples/__init__.py b/samples/__init__.py index 33d7cf6..eb5264c 100644 --- a/samples/__init__.py +++ b/samples/__init__.py @@ -11,66 +11,101 @@ LOG = os.path.join(ROOT, "logs") import salome -from simpleCubic import simpleCubic -from faceCenteredCubic import faceCenteredCubic -from bodyCenteredCubic import bodyCenteredCubic +from simple import simpleCubic, simpleHexagonalPrism +from faceCentered import faceCenteredCubic, faceCenteredHexagonalPrism +from bodyCentered import bodyCenteredCubic, bodyCenteredHexagonalPrism from src import geometry_utils from src import mesh_utils -def genMesh(stype, theta, fillet, flowdirection, saveto): - _G = globals() +def genMesh(stype, theta, fillet, direction, saveto): - structure = _G.get(stype) + logging.info("""genMesh: + structure type:\t{} + coefficient:\t{} + fillet:\t{} + flow direction:\t{} + export path:\t{}""".format(stype, theta, fillet, direction, saveto)) - if structure: - salome.salome_init() + params = (theta, fillet, direction) - grains, geometry1, geometry2 = structure(theta, fillet) - geometry = geometry1 - - if flowdirection == [1, 1, 1]: - geometry = geometry2 - norm = [-1, 1, 0] - bcount = 6 + salome.salome_init() + + ### + # Structure and mesh configurations + ## + if stype == "simple": + if direction in [[1, 0, 0], [0, 0, 1]]: + structure = simpleCubic - # 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] + elif direction == [1, 1, 1]: + structure = simpleHexagonalPrism - if flowdirection == [1, 0, 0]: - normvec = [0, 0, 1] - bcount = 4 - direction = [1, 1, 0] + elif stype == "faceCentered": + if direction in [[1, 0, 0], [0, 0, 1]]: + structure = faceCenteredCubic - if flowdirection == [0, 0, 1]: - normvec = [1, 1, 0] - bcount = 4 - direction = [0, 0, 1] - - # - geometry = geometry_utils.geompy.RemoveExtraEdges(geometry, False) + elif direction == [1, 1, 1]: + structure = faceCenteredHexagonalPrism - # - boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains) + elif stype == "bodyCentered": + if direction in [[1, 0, 0], [0, 0, 1]]: + structure = bodyCenteredCubic - fineness = 1 - viscousLayers = { - "thickness": 0.0001, - "number": 3, - "stretch": 1.2 - } - mesh = mesh_utils.meshCreate(geometry, boundary, fineness, viscousLayers) - mesh_utils.meshCompute(mesh) + elif direction == [1, 1, 1]: + structure = bodyCenteredHexagonalPrism + + ### + # Shape + ## + geompy = geometry_utils.getGeom() + shape, groups = structure(*params) + [length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06) - mesh_utils.meshExport(mesh, saveto) + logging.info("""shape: + surface area:\t{} + volume:\t{}""".format(surfaceArea, volume)) + + ### + # Mesh + ## + fineness = 1 + parameters = mesh_utils.Parameters( + minSize = 0.001, + maxSize = 0.1, + growthRate = 0.1, + nbSegPerEdge = 5, + nbSegPerRadius = 10, + chordalErrorEnabled = False, + chordalError = -1, + secondOrder = False, + optimize = True, + quadAllowed = False, + useSurfaceCurvature = True, + fuseEdges = True, + checkChartBoundary = False + ) + + facesToIgnore = [] + for group in groups: + if group.GetName() in ["inlet", "outlet"]: + facesToIgnore.append(group) - salome.salome_close() + viscousLayers = mesh_utils.ViscousLayers( + thickness = 0.001, + numberOfLayers = 3, + stretchFactor = 1.2, + isFacesToIgnore = True, + facesToIgnore = facesToIgnore, + extrusionMethod = mesh_utils.smeshBuilder.NODE_OFFSET + ) + + mesh = mesh_utils.meshCreate(shape, groups, fineness, parameters, viscousLayers) + mesh_utils.meshCompute(mesh) - else: - raise Exception("Unknown sample function") + mesh_utils.meshExport(mesh, saveto) + + salome.salome_close() if __name__ == "__main__": @@ -89,12 +124,5 @@ if __name__ == "__main__": 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) diff --git a/samples/bodyCentered.py b/samples/bodyCentered.py new file mode 100644 index 0000000..bba7331 --- /dev/null +++ b/samples/bodyCentered.py @@ -0,0 +1,284 @@ +import salome +salome.salome_init() + +import GEOM +from salome.geom import geomBuilder +geompy = geomBuilder.New() + +from math import pi, sqrt + +def bodyCenteredCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]): + + ### + # Parameters + ## + L = 1.0 + r0 = L * sqrt(3) / 4 + + radius = r0 / (1 - theta) + xn, yn, zn = 3, 3, 3 + + length = 2 * r0 + width = L / 2 + diag = L * sqrt(2) + height = L + + point = [] + xl = sqrt(diag ** 2 + diag ** 2) * 0.5 + yw = xl + zh = height + + C1, C2 = 0.8, 0.05 + theta1, theta2 = 0.01, 0.18 + Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1) + filletradius = Cf * (radius - r0) + + scale = 100 + oo = geompy.MakeVertex(0, 0, 0) + spos1 = (0, 0, 0) + spos2 = (0, 0, 0) + + ### + # Bounding box + ## + if direction == [1, 0, 0]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(xl, 0, 0) + sk.addPointsAbsolute(0, yw, 0) + sk.addPointsAbsolute(0, yw, zh) + sk.addPointsAbsolute(xl, 0, zh) + sk.addPointsAbsolute(xl, 0, 0) + + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, diag) + + elif direction == [0, 0, 1]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(0, yw, 0) + sk.addPointsAbsolute(xl, 0, 0) + sk.addPointsAbsolute(2 * xl, yw, 0) + sk.addPointsAbsolute(xl, 2 * yw, 0) + sk.addPointsAbsolute(0, yw, 0) + + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, zh) + + else: + raise Exception("The direction is not implemented") + + inletface = geompy.MakeScaleTransform(inletface, oo, scale) + cubic = geompy.MakeScaleTransform(cubic, oo, scale) + + faces = geompy.ExtractShapes(cubic, geompy.ShapeType["FACE"], False) + symetryface = [] + + for face in faces: + norm = geompy.GetNormal(face) + angle = round(geompy.GetAngle(norm, vecflow), 0) + + if (angle == 0 or angle == 180) and not face == inletface: + outletface = face + + else: + symetryface.append(face) + + ### + # Grains + ## + ox = geompy.MakeVectorDXDYDZ(1, 0, 0) + oy = geompy.MakeVectorDXDYDZ(0, 1, 0) + oz = geompy.MakeVectorDXDYDZ(0, 0, 1) + xy = geompy.MakeVectorDXDYDZ(1, 1, 0) + xmy = geompy.MakeVectorDXDYDZ(1, -1, 0) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius) + lattice1 = geompy.MakeMultiTranslation2D(grain, ox, L, xn, oy, L, yn) + lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn) + + #grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius) + #lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1) + #lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn) + lattice2 = geompy.MakeTranslation(lattice1, 0.5 * L, 0.5 * L, 0.5 * L) + + grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True) + grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True) + grains = geompy.MakeFuseList(grains, False, False) + + grains = geompy.MakeScaleTransform(grains, oo, scale) + + if fillet: + grains = geompy.MakeFilletAll(grains, filletradius * scale) + + ### + # Groups + ## + shape = geompy.MakeCutList(cubic, [grains]) + shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "bodyCenteredCubic") + + sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"]) + geompy.UnionIDs(sall, + geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"])) + + inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet") + inletshape = geompy.MakeCutList(inletface, [grains]) + inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale) + geompy.UnionList(inlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"])) + + outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet") + outletshape = geompy.MakeCutList(outletface, [grains]) + outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale) + geompy.UnionList(outlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"])) + + symetry = [] + for (n, face) in enumerate(symetryface): + name = "symetry" + str(n) + symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name)) + symetryshape = geompy.MakeCutList(face, [grains]) + symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale) + geompy.UnionList(symetry[n], geompy.SubShapeAll( + geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"])) + + groups = [] + groups.append(inlet) + groups.append(outlet) + groups.extend(symetry) + wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) + + return shape, groups + +def bodyCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]): + + ### + # Parameters + ## + L = 1.0 + r0 = L * sqrt(3) / 4 + + radius = r0 / (1 - theta) + xn, yn, zn = 3, 3, 3 + + length = 2 * r0 + width = L / 2 + diag = L * sqrt(2) + height = diag / 3 + + point = [] + xl, yw, zh = -L / 4, -L / 4, -L / 4 + point.append((L / 3 + xl, L / 3 + yw, 4 * L / 3 + zh)) + point.append((L + xl, 0 + yw, L + zh)) + point.append((4 * L / 3 + xl, L / 3 + yw, L / 3 + zh)) + point.append((L + xl, L + yw, 0 + zh)) + point.append((L / 3 + xl, 4 * L / 3 + yw, L / 3 + zh)) + point.append((0 + xl, L + yw, L + zh)) + point.append((L / 3 + xl, L / 3 + yw, 4 * L / 3 + zh)) + + C1, C2 = 0.8, 0.05 + theta1, theta2 = 0.01, 0.18 + Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1) + filletradius = Cf * (radius - r0) + + scale = 100 + oo = geompy.MakeVertex(0, 0, 0) + spos1 = (0, 0, 0) + spos2 = (0, 0, 0) + + ### + # Bounding box + ## + sk = geompy.Sketcher3D() + + for p in point: + sk.addPointsAbsolute(*p) + + inletface = geompy.MakeFaceWires([sk.wire()], False) + vecflow = geompy.GetNormal(inletface) + hexagonPrism = geompy.MakePrismVecH(inletface, vecflow, L * sqrt(3)) + + inletface = geompy.MakeScaleTransform(inletface, oo, scale) + hexagonPrism = geompy.MakeScaleTransform(hexagonPrism, oo, scale) + + faces = geompy.ExtractShapes(hexagonPrism, geompy.ShapeType["FACE"], False) + symetryface = [] + + for face in faces: + norm = geompy.GetNormal(face) + angle = round(geompy.GetAngle(norm, vecflow), 0) + + if (angle == 0 or angle == 180) and not face == inletface: + outletface = face + + else: + symetryface.append(face) + + ### + # Grains + ## + ox = geompy.MakeVectorDXDYDZ(1, 0, 0) + oy = geompy.MakeVectorDXDYDZ(0, 1, 0) + oz = geompy.MakeVectorDXDYDZ(0, 0, 1) + xy = geompy.MakeVectorDXDYDZ(1, 1, 0) + xmy = geompy.MakeVectorDXDYDZ(1, -1, 0) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius) + lattice1 = geompy.MakeMultiTranslation2D(grain, ox, L, xn, oy, L, yn) + lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn) + + #grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius) + #lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1) + #lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn) + lattice2 = geompy.MakeTranslation(lattice1, 0.5 * L, 0.5 * L, 0.5 * L) + + grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True) + grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True) + grains = geompy.MakeFuseList(grains, False, False) + + grains = geompy.MakeScaleTransform(grains, oo, scale) + + if fillet: + grains = geompy.MakeFilletAll(grains, filletradius * scale) + + ### + # Groups + ## + shape = geompy.MakeCutList(hexagonPrism, [grains]) + shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "bodyCenteredCubic") + + sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"]) + geompy.UnionIDs(sall, + geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"])) + + inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet") + inletshape = geompy.MakeCutList(inletface, [grains]) + inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale) + geompy.UnionList(inlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"])) + + outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet") + outletshape = geompy.MakeCutList(outletface, [grains]) + outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale) + geompy.UnionList(outlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"])) + + symetry = [] + for (n, face) in enumerate(symetryface): + name = "symetry" + str(n) + symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name)) + symetryshape = geompy.MakeCutList(face, [grains]) + symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale) + geompy.UnionList(symetry[n], geompy.SubShapeAll( + geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"])) + + groups = [] + groups.append(inlet) + groups.append(outlet) + groups.extend(symetry) + wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) + + return shape, groups + diff --git a/samples/faceCentered.py b/samples/faceCentered.py new file mode 100644 index 0000000..bb32cb4 --- /dev/null +++ b/samples/faceCentered.py @@ -0,0 +1,282 @@ +import salome +salome.salome_init() + +import GEOM +from salome.geom import geomBuilder +geompy = geomBuilder.New() + +from math import pi, sqrt + +def faceCenteredCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]): + + ### + # Parameters + ## + L = 1.0 + r0 = L * sqrt(2) / 4 + + radius = r0 / (1 - theta) + xn, yn, zn = 3, 3, 3 + + length = 2 * r0 + width = L / 2 + diag = L * sqrt(3) + height = diag / 3 + + point = [] + xl = sqrt(length ** 2 + length ** 2) * 0.5 + yw = xl + zh = width + + C1, C2 = 0.8, 0.05 + theta1, theta2 = 0.01, 0.13 + Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1) + filletradius = Cf * (radius - r0) + + scale = 100 + oo = geompy.MakeVertex(0, 0, 0) + spos1 = (-width * (xn - 1), 0, -width * (zn - 2)) + spos2 = (-width * xn, 0, -width * (zn - 1)) + + ### + # Bounding box + ## + if direction == [1, 0, 0]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(0, 0, -zh) + sk.addPointsAbsolute(-xl, yw, -zh) + sk.addPointsAbsolute(-xl, yw, zh) + sk.addPointsAbsolute(0, 0, zh) + sk.addPointsAbsolute(0, 0, -zh) + + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, length) + + elif direction == [0, 0, 1]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(0, 0, -zh) + sk.addPointsAbsolute(xl, yw, -zh) + sk.addPointsAbsolute(0, 2 * yw, -zh) + sk.addPointsAbsolute(-xl, yw, -zh) + sk.addPointsAbsolute(0, 0, -zh) + + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, 2 * zh) + + else: + raise Exception("The direction is not implemented") + + inletface = geompy.MakeScaleTransform(inletface, oo, scale) + cubic = geompy.MakeScaleTransform(cubic, oo, scale) + + faces = geompy.ExtractShapes(cubic, geompy.ShapeType["FACE"], False) + symetryface = [] + + for face in faces: + norm = geompy.GetNormal(face) + angle = round(geompy.GetAngle(norm, vecflow), 0) + + if (angle == 0 or angle == 180) and not face == inletface: + outletface = face + + else: + symetryface.append(face) + + ### + # Grains + ## + ox = geompy.MakeVectorDXDYDZ(1, 0, 0) + oy = geompy.MakeVectorDXDYDZ(0, 1, 0) + oz = geompy.MakeVectorDXDYDZ(0, 0, 1) + xy = geompy.MakeVectorDXDYDZ(1, 1, 0) + xmy = geompy.MakeVectorDXDYDZ(1, -1, 0) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius) + lattice1 = geompy.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn) + lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn - 1) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius) + lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1) + lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn) + + grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True) + grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True) + grains = geompy.MakeFuseList(grains, False, False) + + grains = geompy.MakeScaleTransform(grains, oo, scale) + + if fillet: + grains = geompy.MakeFilletAll(grains, filletradius * scale) + + ### + # Groups + ## + shape = geompy.MakeCutList(cubic, [grains]) + shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "faceCenteredCubic") + + sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"]) + geompy.UnionIDs(sall, + geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"])) + + inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet") + inletshape = geompy.MakeCutList(inletface, [grains]) + inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale) + geompy.UnionList(inlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"])) + + outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet") + outletshape = geompy.MakeCutList(outletface, [grains]) + outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale) + geompy.UnionList(outlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"])) + + symetry = [] + for (n, face) in enumerate(symetryface): + name = "symetry" + str(n) + symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name)) + symetryshape = geompy.MakeCutList(face, [grains]) + symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale) + geompy.UnionList(symetry[n], geompy.SubShapeAll( + geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"])) + + groups = [] + groups.append(inlet) + groups.append(outlet) + groups.extend(symetry) + wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) + + return shape, groups + +def faceCenteredHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]): + + ### + # Parameters + ## + L = 1.0 + r0 = L * sqrt(2) / 4 + + radius = r0 / (1 - theta) + xn, yn, zn = 3, 3, 3 + + length = 2 * r0 + width = L / 2 + diag = L * sqrt(3) + height = diag / 3 + + point = [] + xl, yw, zh = -(xn - 2) * L / 3, -(yn - 2) * L / 3, -(zn - 2) * L / 3 + point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh)) + point.append((0 + xl, -width + yw, 0 + zh)) + point.append((width / 3 + xl, -2 * width / 3 + yw, -2 * width / 3 + zh)) + point.append((0 + xl, 0 + yw, -width + zh)) + point.append((-2 * width / 3 + xl, width / 3 + yw, -2 * width / 3 + zh)) + point.append((-width + xl, 0 + yw, 0 + zh)) + point.append((-2 * width / 3 + xl, -2 * width / 3 + yw, width / 3 + zh)) + + C1, C2 = 0.8, 0.05 + theta1, theta2 = 0.01, 0.13 + Cf = C1 + (C2 - C1) / (theta2 - theta1) * (theta - theta1) + filletradius = Cf * (radius - r0) + + scale = 100 + oo = geompy.MakeVertex(0, 0, 0) + spos1 = (-width * (xn - 1), 0, -width * (zn - 2)) + spos2 = (-width * xn, 0, -width * (zn - 1)) + + ### + # Bounding box + ## + sk = geompy.Sketcher3D() + + for p in point: + sk.addPointsAbsolute(*p) + + inletface = geompy.MakeFaceWires([sk.wire()], False) + vecflow = geompy.GetNormal(inletface) + hexagonPrism = geompy.MakePrismVecH(inletface, vecflow, L * sqrt(3)) + + inletface = geompy.MakeScaleTransform(inletface, oo, scale) + hexagonPrism = geompy.MakeScaleTransform(hexagonPrism, oo, scale) + + faces = geompy.ExtractShapes(hexagonPrism, geompy.ShapeType["FACE"], False) + symetryface = [] + + for face in faces: + norm = geompy.GetNormal(face) + angle = round(geompy.GetAngle(norm, vecflow), 0) + + if (angle == 0 or angle == 180) and not face == inletface: + outletface = face + + else: + symetryface.append(face) + + ### + # Grains + ## + ox = geompy.MakeVectorDXDYDZ(1, 0, 0) + oy = geompy.MakeVectorDXDYDZ(0, 1, 0) + oz = geompy.MakeVectorDXDYDZ(0, 0, 1) + xy = geompy.MakeVectorDXDYDZ(1, 1, 0) + xmy = geompy.MakeVectorDXDYDZ(1, -1, 0) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos1), radius) + lattice1 = geompy.MakeMultiTranslation2D(grain, xy, length, xn, xmy, length, yn) + lattice1 = geompy.MakeMultiTranslation1D(lattice1, oz, L, zn - 1) + + grain = geompy.MakeSpherePntR(geompy.MakeVertex(*spos2), radius) + lattice2 = geompy.MakeMultiTranslation2D(grain, xy, length, xn + 1, xmy, length, yn + 1) + lattice2 = geompy.MakeMultiTranslation1D(lattice2, oz, L, zn) + + grains = geompy.ExtractShapes(lattice1, geompy.ShapeType["SOLID"], True) + grains += geompy.ExtractShapes(lattice2, geompy.ShapeType["SOLID"], True) + grains = geompy.MakeFuseList(grains, False, False) + + grains = geompy.MakeScaleTransform(grains, oo, scale) + + if fillet: + grains = geompy.MakeFilletAll(grains, filletradius * scale) + + ### + # Groups + ## + shape = geompy.MakeCutList(hexagonPrism, [grains]) + shape = geompy.MakeScaleTransform(shape, oo, 1 / scale, theName = "faceCenteredCubic") + + sall = geompy.CreateGroup(shape, geompy.ShapeType["FACE"]) + geompy.UnionIDs(sall, + geompy.SubShapeAllIDs(shape, geompy.ShapeType["FACE"])) + + inlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "inlet") + inletshape = geompy.MakeCutList(inletface, [grains]) + inletshape = geompy.MakeScaleTransform(inletshape, oo, 1 / scale) + geompy.UnionList(inlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, inletshape, True), geompy.ShapeType["FACE"])) + + outlet = geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = "outlet") + outletshape = geompy.MakeCutList(outletface, [grains]) + outletshape = geompy.MakeScaleTransform(outletshape, oo, 1 / scale) + geompy.UnionList(outlet, geompy.SubShapeAll( + geompy.GetInPlace(shape, outletshape, True), geompy.ShapeType["FACE"])) + + symetry = [] + for (n, face) in enumerate(symetryface): + name = "symetry" + str(n) + symetry.append(geompy.CreateGroup(shape, geompy.ShapeType["FACE"], theName = name)) + symetryshape = geompy.MakeCutList(face, [grains]) + symetryshape = geompy.MakeScaleTransform(symetryshape, oo, 1 / scale) + geompy.UnionList(symetry[n], geompy.SubShapeAll( + geompy.GetInPlace(shape, symetryshape, True), geompy.ShapeType["FACE"])) + + groups = [] + groups.append(inlet) + groups.append(outlet) + groups.extend(symetry) + wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) + + return shape, groups + diff --git a/samples/test.py b/samples/simple.py similarity index 85% rename from samples/test.py rename to samples/simple.py index 2cfef0f..2e98b56 100644 --- a/samples/test.py +++ b/samples/simple.py @@ -7,7 +7,7 @@ geompy = geomBuilder.New() from math import pi, sqrt -def simpleCubic(theta = 0.01, fillet = False): +def simpleCubic(theta = 0.01, fillet = False, direction = [1, 0, 0]): ### # Parameters @@ -37,17 +37,33 @@ def simpleCubic(theta = 0.01, fillet = False): ### # Bounding box ## - sk = geompy.Sketcher3D() - sk.addPointsAbsolute(xl, 0, 0) - sk.addPointsAbsolute(0, yw, 0) - sk.addPointsAbsolute(0, yw, zh) - sk.addPointsAbsolute(xl, 0, zh) - sk.addPointsAbsolute(xl, 0, 0) - - inletface = geompy.MakeFaceWires([sk.wire()], True) - vecflow = geompy.GetNormal(inletface) - cubic = geompy.MakePrismVecH(inletface, vecflow, width) + if direction == [1, 0, 0]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(xl, 0, 0) + sk.addPointsAbsolute(0, yw, 0) + sk.addPointsAbsolute(0, yw, zh) + sk.addPointsAbsolute(xl, 0, zh) + sk.addPointsAbsolute(xl, 0, 0) + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, width) + + elif direction == [0, 0, 1]: + sk = geompy.Sketcher3D() + sk.addPointsAbsolute(0, yw, 0) + sk.addPointsAbsolute(xl, 0, 0) + sk.addPointsAbsolute(2 * xl, yw, 0) + sk.addPointsAbsolute(xl, 2 * yw, 0) + sk.addPointsAbsolute(0, yw, 0) + + inletface = geompy.MakeFaceWires([sk.wire()], True) + vecflow = geompy.GetNormal(inletface) + cubic = geompy.MakePrismVecH(inletface, vecflow, height) + + else: + raise Exception("The direction is not implemented") + inletface = geompy.MakeScaleTransform(inletface, oo, scale) cubic = geompy.MakeScaleTransform(cubic, oo, scale) @@ -119,11 +135,12 @@ def simpleCubic(theta = 0.01, fillet = False): groups.append(outlet) groups.extend(symetry) wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) - return shape + return shape, groups -def simpleHexagonalPrism(theta = 0.01, fillet = False): +def simpleHexagonalPrism(theta = 0.01, fillet = False, direction = [1, 1, 1]): ### # Parameters @@ -239,6 +256,7 @@ def simpleHexagonalPrism(theta = 0.01, fillet = False): groups.append(outlet) groups.extend(symetry) wall = geompy.CutListOfGroups([sall], groups, theName = "wall") + groups.append(wall) + + return shape, groups - return shape - diff --git a/src/mesh_utils.py b/src/mesh_utils.py index ad5f073..ca23bff 100644 --- a/src/mesh_utils.py +++ b/src/mesh_utils.py @@ -3,12 +3,38 @@ from salome.smesh import smeshBuilder smesh = smeshBuilder.New() import logging +from collections import namedtuple + +Parameters = namedtuple("Parameters", [ + "minSize", + "maxSize", + "growthRate", + "nbSegPerEdge", + "nbSegPerRadius", + "chordalErrorEnabled", + "chordalError", + "secondOrder", + "optimize", + "quadAllowed", + "useSurfaceCurvature", + "fuseEdges", + "checkChartBoundary" +]) + +ViscousLayers = namedtuple("ViscousLayers", [ + "thickness", + "numberOfLayers", + "stretchFactor", + "isFacesToIgnore", + "facesToIgnore", + "extrusionMethod" +]) def getSmesh(): return smesh -def meshCreate(gobj, boundary, fineness, viscousLayers=None): +def meshCreate(shape, groups, fineness, parameters, viscousLayers = None): """ Creates a mesh from a geometry. @@ -21,19 +47,13 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None): 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. """ + ### + # Netgen + ## Fineness = { 0: "Very coarse", 1: "Coarse", @@ -43,31 +63,34 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None): 5: "Custom" }[fineness] - mesh = smesh.Mesh(gobj) + # Mesh + mesh = smesh.Mesh(shape) netgen = mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D) + # Parameters 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(parameters.minSize) + param.SetMaxSize(parameters.maxSize) param.SetFineness(fineness) - # TODO: add customization if fineness == 5: - param.SetGrowthRate( 0.1 ) - param.SetNbSegPerEdge( 5 ) - param.SetNbSegPerRadius( 10 ) + param.SetGrowthRate(parameters.growthRate) + param.SetNbSegPerEdge(parameters.nbSegPerEdge) + param.SetNbSegPerRadius(parameters.nbSegPerRadius) + param.SetChordalErrorEnabled(parameters.chordalErrorEnabled) + param.SetChordalError(parameters.chordalError) + + param.SetSecondOrder(parameters.secondOrder) + param.SetOptimize(parameters.optimize) + param.SetQuadAllowed(parameters.quadAllowed) + + param.SetUseSurfaceCurvature(parameters.useSurfaceCurvature) + param.SetFuseEdges(parameters.fuseEdges) + param.SetCheckChartBoundary(parameters.checkChartBoundary) + + logging.info("""meshCreate: fineness:\t{} min size:\t{} @@ -86,29 +109,41 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None): True if param.GetSecondOrder() else False, True if param.GetOptimize() else False)) + + ### + # Groups + ## + for group in groups: + mesh.GroupOnGeom(group, "{}_".format(group.GetName()), SMESH.FACE) - + ### + # Viscous layers + ## if not viscousLayers is None: + vlayer = netgen.ViscousLayers( + viscousLayers.thickness, + viscousLayers.numberOfLayers, + viscousLayers.stretchFactor, + viscousLayers.facesToIgnore, + viscousLayers.isFacesToIgnore, + viscousLayers.extrusionMethod + ) + logging.info("""meshCreate: viscous layers: thickness:\t{} number:\t{} stretch factor:\t{}""".format( - viscousLayers["thickness"], viscousLayers["number"], viscousLayers["stretch"])) + viscousLayers.thickness, + viscousLayers.numberOfLayers, + viscousLayers.stretchFactor)) - vlayer = netgen.ViscousLayers( - viscousLayers["thickness"], - viscousLayers["number"], - viscousLayers["stretch"], - [boundary["inlet"], boundary["outlet"]], - 1, smeshBuilder.NODE_OFFSET) + else: logging.info("""meshCreate: viscous layers: disabled""") - for name, b in boundary.items(): - mesh.GroupOnGeom(b, "{}_".format(name), SMESH.FACE) return mesh @@ -175,7 +210,7 @@ def meshExport(mobj, path): mobj.ExportUNV(path) logging.info("""meshExport: - format:\t{}""".format("unv")) + format:\t{}""".format("unv")) except: logging.error("""meshExport: Cannot export.""") diff --git a/temp/samples.old.2/__init__.py b/temp/samples.old.2/__init__.py new file mode 100644 index 0000000..09764e6 --- /dev/null +++ b/temp/samples.old.2/__init__.py @@ -0,0 +1,100 @@ +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) + diff --git a/samples/bodyCenteredCubic.py b/temp/samples.old.2/bodyCenteredCubic.py similarity index 100% rename from samples/bodyCenteredCubic.py rename to temp/samples.old.2/bodyCenteredCubic.py diff --git a/samples/faceCenteredCubic.py b/temp/samples.old.2/faceCenteredCubic.py similarity index 100% rename from samples/faceCenteredCubic.py rename to temp/samples.old.2/faceCenteredCubic.py diff --git a/samples/simpleCubic.py b/temp/samples.old.2/simpleCubic.py similarity index 100% rename from samples/simpleCubic.py rename to temp/samples.old.2/simpleCubic.py