#!/usr/bin/env python # -*- coding: utf-8 -*- import salome, GEOM, SMESH, SALOMEDS from salome.geom import geomBuilder from salome.smesh import smeshBuilder import math import os, sys import logging import time from datetime import timedelta class simpleCubic: def __init__(self, name = None): self.name = name if type(name) != None else "simpleCubic" self.geometry = None self.geometrybbox = None self.mesh = None self.boundary = None self.rombus = None self.rombusbbox = None 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() # R_0 = 1 R = R_0 / (1 - alpha) R_fillet = 0.05 # 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], R) sphere = geompy.MakeMultiTranslation2D(sphere, None, 2, 3, None, 2, 3) sphere = geompy.MakeTranslation(sphere, -2, 0, 0) sphere2 = geompy.MakeTranslation(sphere, 0, 0, 2) sphere3 = geompy.MakeTranslation(sphere2, 0, 0, 2) sphere = geompy.ExtractShapes(sphere, geompy.ShapeType["SOLID"], True) sphere2 = geompy.ExtractShapes(sphere2, geompy.ShapeType["SOLID"], True) sphere3 = geompy.ExtractShapes(sphere3, geompy.ShapeType["SOLID"], True) sphere = geompy.MakeFuseList(sphere + sphere2 + sphere3, True, True) sphere = geompy.MakeFilletAll(sphere, R_fillet) self.geometry = geompy.MakeCutList(box, [sphere], True) self.geometrybbox = box geompy.addToStudy(self.geometry, self.name) # 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 geompy.addToStudy(self.rombus, "rombus") return self.geometry def boundaryCreate(self, direction): """ 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": } """ # # _____ 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 if direction == "001": center = geompy.MakeVertex(2, 2, 1) norm = geompy.MakeVector(center, geompy.MakeVertexWithRef(center, 0, 0, 1)) vstep = 1 hstep = math.sqrt(2) elif direction == "100": center = geompy.MakeVertex(2, 2, 1) 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), 0)) vstep = math.sqrt(2) hstep = 1 elif direction == "111": center = geompy.MakeVertex(2, 2, 2) 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)) vstep = math.sqrt(2) hstep = 1 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 else: box = self.geometrybbox planes = geompy.ExtractShapes(box, geompy.ShapeType["FACE"], True) vplanes = [] hplanes = [] n = 0 for plane in planes: planeNorm = geompy.GetNormal(plane) n += 1 geompy.addToStudy(planeNorm, "normalplane-{}".format(n)) angle = int(abs(geompy.GetAngle(planeNorm, norm))) logging.info("angle = {}".format(angle)) if angle == 0 or angle == 180: vplanes.append(plane) else: hplanes.append(plane) if salome.sg.hasDesktop(): salome.sg.updateObjBrowser() logging.info(len(vplanes)) logging.info(len(hplanes)) if direction == "001": z1 = geompy.GetPosition(vplanes[0])[3] z2 = geompy.GetPosition(vplanes[1])[3] if z1 > z2: inletplane = vplanes[0] outletplane = vplanes[1] else: inletplane = vplanes[1] outletplane = vplanes[0] elif direction == "100" or direction == "111": x1 = geompy.GetPosition(vplanes[0])[1] x2 = geompy.GetPosition(vplanes[1])[1] logging.info("x1 = {}, x2 = {}".format(x1, x2)) if x1 > x2: inletplane = vplanes[0] outletplane = vplanes[1] else: inletplane = vplanes[1] outletplane = vplanes[0] # inlet and outlet common1 = geompy.MakeCommonList([self.geometry, inletplane], True) inlet = createGroup(common1, "inlet") 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) faces = geompy.SubShapeAll(gip, geompy.ShapeType["FACE"]) geompy.UnionList(symetryPlane, faces) # wall allgroup = geompy.CreateGroup(self.geometry, geompy.ShapeType["FACE"]) faces = geompy.SubShapeAllIDs(self.geometry, geompy.ShapeType["FACE"]) geompy.UnionIDs(allgroup, faces) wall = geompy.CutListOfGroups([allgroup], [inlet, outlet, symetryPlane], "wall") self.boundary = { "inlet": inlet, "outlet": outlet, "symetryPlane": symetryPlane, "wall": wall } return self.boundary 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(self.geometry) 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(fineness) #param.SetGrowthRate( 0.1 ) #param.SetNbSegPerEdge( 5 ) #param.SetNbSegPerRadius( 10 ) param.SetQuadAllowed( 0 ) 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: logging.info("Mesh succesfully computed.") else: logging.warning("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: self.mesh.ExportUNV(exportpath) except: logging.error("Cannot export mesh to '{}'".format(exportpath)) if __name__ == "__main__": # Arguments buildpath = str(sys.argv[1]) alpha = float(sys.argv[2]) direction = str(sys.argv[3]) name = "simpleCubic-{}-{}".format(direction, alpha) # Logger logging.basicConfig( level=logging.INFO, format="%(levelname)s: %(message)s", handlers = [ logging.StreamHandler(), logging.FileHandler(os.path.join(buildpath, "{}.log".format(name))) ]) start_time = time.monotonic() # Simple cubic sc = simpleCubic(name) logging.info("Creating the geometry ...") sc.geometryCreate(alpha) logging.info("Extracting boundaries ...") sc.boundaryCreate(direction) logging.info("Creating the mesh ...") sc.meshCreate(2, { "thickness": 0.02, "number": 2, "stretch": 1.1 }) sc.meshCompute() logging.info("Exporting the mesh ...") sc.meshExport(buildpath) end_time = time.monotonic() logging.info("Elapsed time: {}".format(timedelta(seconds=end_time - start_time))) logging.info("Done.") if salome.sg.hasDesktop(): salome.sg.updateObjBrowser()