Source code for idwarp.UnstructuredMesh

"""Unstructured Mesh

The UnstructuredMesh module is used for interacting with an
unstructured (or structured!)  mesh - typically used in a 3D CFD
program.

It contains the following classes:

USMesh: General class for working with unstructured meshes

Copyright (c) 2014 by C.A.(Sandy) Mader
All rights reserved. Not to be used for commercial purposes.
Revision: 1.0   $Date: 09/07/2014$

Developers:
-----------
- C.A.(Sandy) Mader (CAM)
- Gaetan. K. W. Kenway (GKK)
History
-------
    v. 1.0 - Initial Class Creation (CAM, 2014)
"""
# =============================================================================
# Imports
# =============================================================================
import os
import shutil
import warnings
import numpy as np
from mpi4py import MPI
from .MExt import MExt
from baseclasses import BaseSolver
from baseclasses.utils import Error


# =============================================================================
# UnstructuredMesh class
# =============================================================================


[docs]class USMesh(BaseSolver): """ This is the main Unstructured Mesh. This mesh object is designed to interact with an structured or unstructured CFD solver though a variety of interface functions. """ def __init__(self, options=None, comm=None, debug=False): """ Create the USMesh object. Parameters ---------- options : dictionary A dictionary containing the options for the mesh movement strategies. THIS IS NOT OPTIONAL. comm : MPI_INTRA_COMM MPI communication (as obtained from mpi4py) on which to create the USMesh object. If not provided, MPI_COMM_WORLD is used. debug : bool Flag specifying if the MExt import is automatically deleted. This needs to be true ONLY when a symbolic debugger is used. """ name = "IDWarp" category = "Volume mesh warping" if comm is None: comm = MPI.COMM_WORLD # Default options for mesh warping defOpts = self._getDefaultOptions() if options is None: raise Error( "The 'options' keyword argument is *NOT* " "optional. An options dictionary must be passed upon " "creation of this object" ) # Initialize the inherited BaseSolver super().__init__(name, category, defaultOptions=defOpts, options=options, comm=comm) self.printOptions() # Check if warp has already been set if this has been # inherited to complex version try: self.warp except AttributeError: curDir = os.path.basename(os.path.dirname(os.path.realpath(__file__))) self.warp = MExt("libidwarp", curDir, debug=debug)._module # Initialize PETSc if not done so self.warp.initpetsc(self.comm.py2f()) # Set realtype of 'd'. 'D' is used in Complex and set in # UnstructuredMesh_C.py self.dtype = "d" # Set Fortran options values self._setMeshOptions() # Initialize various bits of stored information self.OFData = {} self.warpInitialized = False self.faceSizes = None self.fileType = self.getOption("fileType") fileName = self.getOption("gridFile") # Determine how to read if self.fileType == "CGNS": # Determine type of CGNS mesh we have self.warp.readcgns(fileName) elif self.fileType == "OpenFOAM": self._readOFGrid(fileName) elif self.fileType == "PLOT3D": self.warp.readplot3d(fileName) @staticmethod def _getDefaultOptions(): defOpts = { "gridFile": [(str, type(None)), None], "fileType": [str, ["CGNS", "OpenFOAM", "PLOT3D"]], "specifiedSurfaces": [(list, type(None)), None], "symmetrySurfaces": [(list, type(None)), None], "symmetryPlanes": [(list, type(None)), None], "aExp": [float, 3.0], "bExp": [float, 5.0], "LdefFact": [float, 1.0], "alpha": [float, 0.25], "errTol": [float, 0.0005], "evalMode": [str, ["fast", "exact"]], "symmTol": [float, 1e-6], "useRotations": [bool, True], "zeroCornerRotations": [bool, True], "cornerAngle": [float, 30.0], "restartFile": [(str, type(None)), None], "bucketSize": [int, 8], } return defOpts
[docs] def getSurfaceCoordinates(self): """Returns all defined surface coordinates on this processor Returns ------- coords : numpy array size (N,3) Specified surface coordinates residing on this processor. This may be empty array, size (0,3) """ self._setInternalSurface() coords = np.zeros((self.nSurf, 3), self.dtype) self.warp.getsurfacecoordinates(np.ravel(coords)) return coords
[docs] def setSurfaceCoordinates(self, coordinates): """Sets all surface coordinates on this processor Parameters ---------- coordinates : numpy array, size(N, 3) The coordinate to set. This MUST be exactly the same size as the array obtained from getSurfaceCoordinates() """ if len(coordinates) != self.nSurf: raise Error( "Incorrect length of coordinates supplied to " "setSurfaceCoordinates on proc %d. Expected " "aray of length %d, received an array of length " "%d." % (self.comm.rank, self.nSurf, len(coordinates)) ) self.warp.setsurfacecoordinates(np.ravel(coordinates))
[docs] def setExternalMeshIndices(self, ind): """ Set the indicies defining the transformation of an external solver grid to the original CGNS grid. This is required to use USMesh functions that involve the word "Solver" and warpDeriv. The indices must be zero-based. Parameters ---------- ind : numpy integer array The list of indicies this processor needs from the common mesh file """ self.warp.setexternalmeshindices(ind)
[docs] def setSurfaceDefinition(self, pts, conn, faceSizes, cgnsBlockID=None): """This is the master function that determines the definition of the surface to be used for the mesh movement. This surface may be supplied from an external solver (such as SUmb) or it may be generated by IDWarp internally. Parameters ---------- pts : array, size (M, 3) Nodes on this processor conn : int array, size (sum(faceSizes)) Connectivity of the nodes on this processor faceSizes : int Array size (N) Treat the conn array as a flat list with the faceSizes giving connectivity offset for each element. cgnsBlockID : dummy argument. This argument is not used at all. It is here just to have the same API as the IDWarpMulti class. """ conn = np.array(conn).flatten() # Call the fortran initialze warping command with the # coordinates and the patch connectivity given to us. if self.getOption("restartFile") is None: restartFile = "" else: restartFile = self.getOption("restartFile") # The symmetry conditions but be set before initializing the # warping. self._setSymmetryConditions() allPts = np.vstack(self.comm.allgather(pts)) allFaceSizes = np.hstack(self.comm.allgather(faceSizes)) # Be careful with conn...need to increment by the offset from # the points. ptSizes = self.comm.allgather(len(pts)) offsets = np.zeros(len(ptSizes), "intc") offsets[1:] = np.cumsum(ptSizes)[:-1] allConn = np.hstack(self.comm.allgather(conn + offsets[self.comm.rank])) # Next point reduce all nodes: uniquePts, link, nUnique = self.warp.pointreduce(allPts.T, 1e-12) uniquePts = uniquePts.T[0:nUnique] # Now sort them by z, y and x to ensure that the tree is # always created the same way. It is currently sensitivte to # the ordering of the original nodes. Ie, if you reorder the # nodes you get a differnet tree which gives slightly # differnet warp values/derivatives. ind = np.lexsort(uniquePts.T) uniquePts = uniquePts[ind] # Update the link array such that it points back to the original list. inv = np.zeros_like(ind, "intc") rng = np.arange(len(ind)).astype("intc") inv[ind[rng]] = rng link = inv[link - 1] + 1 self.warp.initializewarping(np.ravel(pts), uniquePts.T, link, allFaceSizes, allConn, restartFile) self.nSurf = len(pts) self.warpInitialized = True # Clean-up patch data no longer necessary. self.warp.deallocatepatchio()
[docs] def setSurfaceDefinitionFromFile(self, surfFile): """Set the surface definition for the warping from a multiblock PLOT3D surface file Parameters ---------- surfFile : filename of multiblock PLOT3D surface file. """ # Read the PLOT3D surface file self.warp.readplot3dsurface(surfFile) if self.comm.rank == 0: # Extract out the data we want from the module. Note that # becuase we are in python, convert conn to 0 based. pts = self.warp.plot3dsurface.pts.T.copy() conn = self.warp.plot3dsurface.conn.T.copy() - 1 faceSizes = 4 * np.ones(len(conn)) else: pts = np.zeros((0, 3)) conn = np.zeros((0, 4)) faceSizes = 4 * np.ones(0) # Run the regular setSurfaceDefinition routine. Note that only # the root proc has non-zero length arrays. That's ok. self.setSurfaceDefinition(pts, conn, faceSizes) # Cleanup the fortran memory we used to read the surface self.warp.plot3dsurface.pts = None self.warp.plot3dsurface.conn = None
[docs] def setSurfaceFromFile(self, surfFile): """Update the internal surface surface coordinates using an external PLOT3D surface file. This can be used in an analogous way to setSurfaceDefinitionFromFile. The 'sense' of the file must be same as the file used with setSurfaceDefinitionFromFile. That means, the same number of blocks, with the same sizes, in the same order. Parameters ---------- surfFile: filename of multiblock PLOT3D surface file' """ # Read the PLOT3D surface file self.warp.readplot3dsurface(surfFile) if self.comm.rank == 0: # Extract out the data we want from the module pts = self.warp.plot3dsurface.pts.T.copy() else: pts = np.zeros((0, 3)) # Do the regular update self.setSurfaceCoordinates(pts) # Cleanup the memory we used to read the surface self.warp.plot3dsurface.pts = None self.warp.plot3dsurface.conn = None
[docs] def getSolverGrid(self): """Return the current grid in the order specified by setExternalMeshIndices(). This is the main routine for returning the defomed mesh to the external CFD solver. Returns ------- solverGrid, numpy array, real: The resulting grid. The output is returned in flatted 1D coordinate format. The len of the array is 3*len(indices) as set by setExternalMeshIndices() """ solverGrid = np.zeros(self.warp.griddata.solvermeshdof, self.dtype) warpGrid = self.getWarpGrid() self.warp.warp_to_solver_grid(warpGrid, solverGrid) return solverGrid
[docs] def getWarpGrid(self): """ Return the current grid. This funtion is typically unused. See getSolverGrid for the more useful interface functionality. Returns ------- warpGrid, numpy array, real: The resulting grid. The output is returned in flatted 1D coordinate format. """ return self.warp.getvolumecoordinates(self.warp.griddata.warpmeshdof)
[docs] def getCommonGrid(self): """Return the grid in the original ordering. This is required for the OpenFOAM tecplot writer since the connectivity is only known in this ordering. """ return self.warp.getcommonvolumecoordinates(self.warp.griddata.commonmeshdof)
[docs] def getdXs(self): """Return the current values in dXs. This is the result from a mesh-warp derivative computation. Returns ------- dXs : numpy array The specific components of dXs. size(N,3). This the same size as the array obtained with getSurfaceCoordiantes(). N may be zero if this processor does not have any surface coordinates. """ dXs = np.zeros((self.nSurf, 3), self.dtype) self.warp.getdxs(np.ravel(dXs)) return dXs
[docs] def warpMesh(self): """ Run the applicable mesh warping strategy. This will update the volume coordinates to match surface coordinates set with setSurfaceCoordinates() """ # Initialize internal surface if one isn't already set. self._setInternalSurface() # Set Options self._setMeshOptions() # Now run mesh warp command self.warp.warpmesh()
[docs] def warpDeriv(self, dXv, solverVec=True): """Compute the warping derivative (dXv/dXs^T)*Vec (where vec is the dXv argument to this function. This is the main routine to compute the mesh warping derivative. Parameters ---------- dXv : numpy array Vector of size external solver_grid. This is typically obtained from the external solver's dRdx^T * psi calculation. solverVec : logical Flag to indicate that the dXv vector is in the solver ordering and must be converted to the warp ordering first. This is the usual approach and thus defaults to True. Returns ------- None. The resulting calculation is available from the getdXs() function. """ if solverVec: dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, self.dtype) self.warp.solver_to_warp_grid(dXv, dXvWarp) else: dXvWarp = dXv self.warp.warpderiv(dXvWarp)
[docs] def warpDerivFwd(self, dXs, solverVec=True): """Compute the forward mode warping derivative This routine is not used for "regular" optimization; it is used for matrix-free type optimization. dXs is assumed to be the the peturbation on all the surface nodes. Parameters ---------- dXs : array, size Nsx3 This is the forward mode peturbation seed. Same size as the surface mesh from getSurfaceCoordinates(). solverVec : bool Whether or not to convert to the solver ordering. Returns ------- dXv : array The peturbation on the volume meshes. It may be in warp ordering or solver ordering depending on the solverVec flag. """ dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, dtype=self.dtype) self.warpMesh() self.warp.warpderivfwd(np.ravel(dXs), dXvWarp) if solverVec: dXv = np.zeros(self.warp.griddata.solvermeshdof, dtype=self.dtype) self.warp.warp_to_solver_grid(dXvWarp, dXv) return dXv else: return dXvWarp
[docs] def verifyWarpDeriv(self, dXv=None, solverVec=True, dofStart=0, dofEnd=10, h=1e-6, randomSeed=314): """Run an internal verification of the solid warping derivatives""" if dXv is None: np.random.seed(randomSeed) # 'Random' seed to ensure runs are same dXvWarp = np.random.random(self.warp.griddata.warpmeshdof, dtype=self.dtype) else: if solverVec: dXvWarp = np.zeros(self.warp.griddata.warpmeshdof, self.dtype) self.warp.solver_to_warp_grid(dXv, dXvWarp) else: dXvWarp = dXv self.warp.verifywarpderiv(dXvWarp, dofStart, dofEnd, h)
# ========================================================================== # Output Functionality # ==========================================================================
[docs] def writeGrid(self, fileName=None): """ Write the current grid to the correct format Parameters ---------- fileName : str or None Filename for grid. Should end in .cgns for CGNS files. For PLOT3D whatever you want. It is not optional for CGNS/PLOT3D. It is not required for OpenFOAM meshes. This call will update the 'points' file. """ if self.fileType in ["CGNS", "PLOT3D"]: if fileName is None: raise Error("fileName is not optional for writeGrid with " "gridType of CGNS or PLOT3D") if self.fileType == "CGNS": # Copy the default and then write if self.comm.rank == 0: shutil.copy(self.getOption("gridFile"), fileName) self.warp.writecgns(fileName) elif self.fileType == "PLOT3D": self.warp.writeplot3d(fileName) elif self.fileType == "OpenFOAM": self._writeOpenFOAMVolumePoints(self.getCommonGrid())
[docs] def writeOFGridTecplot(self, fileName): """ Write the current OpenFOAM grid to a Tecplot FE polyhedron file. This is generally used for debugging/visualization purposes. Parameters ---------- fileName : str Filename to use. Should end in .dat for tecplot ascii file. """ if not self.OFData: raise ValueError("Cannot write OpenFOAM tecplot file since there is " "no OpenFOAM data present") if self.comm.size == 1: f = open(fileName, "w") else: fname, fext = os.path.splitext(fileName) fileName = fname + "%d" % self.comm.rank + fext f = open(fileName, "w") # Extract the data we need from the OFDict to make the code a # little easier to read: nodes = self.getCommonGrid() nodes = nodes.reshape((len(nodes) // 3, 3)) nPoints = len(nodes) faces = self.OFData["faces"] nFaces = len(faces) owner = self.OFData["owner"] nCells = np.max(owner) + 1 neighbour = self.OFData["neighbour"] nNeighbours = len(neighbour) # write the node numbers for each face faceNodeSum = 0 for i in range(nFaces): faceNodeSum += len(faces[i]) # Tecplot Header Info: f.write('TITLE = "Example Grid File"\n') f.write("FILETYPE = GRID\n") f.write('VARIABLES = "X" "Y" "Z"\n') f.write("Zone\n") f.write("ZoneType=FEPOLYHEDRON\n") f.write("NODES=%d\n" % nPoints) f.write("FACES=%d\n" % nFaces) f.write("ELEMENTS=%d\n" % nCells) f.write("TotalNumFaceNodes=%d\n" % faceNodeSum) f.write("NumConnectedBoundaryFaces=%d\n" % 0) f.write("TotalNumBoundaryConnections=%d\n" % 0) # Write the points to file in block data format for idim in range(3): for j in range(nPoints): f.write("%g\n" % nodes[j, idim]) # Write the number of nodes in each face for i in range(nFaces): f.write("%d " % len(faces[i])) if i % 300 == 0: f.write("\n") f.write("\n") # write the node numbers for each face (+1 for tecplot 1-based # ordering) for i in range(nFaces): nPointsFace = len(faces[i]) for j in range(nPointsFace): f.write("%d " % (faces[i][j] + 1)) f.write("\n") # write left elements (+1 for 1 based ordering) for i in range(nFaces): f.write("%d\n" % (owner[i] + 1)) # write right elements (+1 for 1 based ordering) for i in range(nFaces): if i < nNeighbours: f.write("%d\n" % (neighbour[i] + 1)) else: f.write("%d\n" % (0)) f.close()
# ========================================================================= # Internal Private Functions # ========================================================================= def _setInternalSurface(self): """This function is used by default if setSurfaceDefinition() is not set BEFORE an operation is requested that requires this information. """ if self.warpInitialized: return if self.comm.rank == 0: warnings.warn( "Using internally generated IDWarp surfaces. If " "this mesh object is to be used with an " "external solver, ensure the mesh object is " "passed to the solver immediatedly after it is created. " "The external solver must then call " "'setExternalMeshIndices()' and 'setSurfaceDefinition()' " "routines.", stacklevel=2, ) conn = [] pts = [] faceSizes = [] if self.fileType == "CGNS": if self.comm.rank == 0: # Do the necessary fortran preprocessing if self.warp.cgnsgrid.cgnsstructured: self.warp.processstructuredpatches() else: self.warp.processunstructuredpatches() fullConn = self.warp.cgnsgrid.surfaceconn - 1 fullPts = self.warp.cgnsgrid.surfacepoints nPatch = self.warp.cgnsgrid.getnpatch() fullPatchNames = [] for i in range(nPatch): fullPatchNames.append(self.warp.cgnsgrid.getsurf(i + 1).strip().lower()) # We now have all surfaces belonging to # boundary conditions. We need to decide which # ones to use depending on what the user has # told us. surfaceFamilies = set() if self.getOption("specifiedSurfaces") is None: # Use all wall surfaces: for i in range(len(self.warp.cgnsgrid.surfaceiswall)): if self.warp.cgnsgrid.surfaceiswall[i]: surfaceFamilies.add(fullPatchNames[i].lower()) else: # The user has supplied a list of surface families for name in self.getOption("specifiedSurfaces"): surfaceFamilies.add(name.lower()) usedFams = set() if self.warp.cgnsgrid.cgnsstructured: if self.comm.rank == 0: # Pull out data and convert to 0-based ordering fullPatchSizes = self.warp.cgnsgrid.surfacesizes.T # Now we loop back over the "full" versions of # thing and just take the ones that correspond # to the families we are using. curNodeIndex = 0 curCellIndex = 0 curOffset = 0 for i in range(len(fullPatchNames)): curNodeSize = fullPatchSizes[i][0] * fullPatchSizes[i][1] curCellSize = (fullPatchSizes[i][0] - 1) * (fullPatchSizes[i][1] - 1) if fullPatchNames[i] in surfaceFamilies: # Keep track of the families we've actually used usedFams.add(fullPatchNames[i]) pts.extend(fullPts[curNodeIndex : curNodeIndex + curNodeSize * 3]) conn.extend(fullConn[curCellIndex : curCellIndex + curCellSize * 4] - curOffset) else: # If we skipped, we increment the offset curOffset += curNodeSize curNodeIndex += curNodeSize * 3 curCellIndex += curCellSize * 4 # end for (root proc) # Run the common surface definition routine pts = np.array(pts).reshape((len(pts) // 3, 3)) faceSizes = 4 * np.ones(len(conn) // 4, "intc") self.setSurfaceDefinition(pts=pts, conn=np.array(conn, "intc"), faceSizes=faceSizes) else: # unstructured if self.comm.rank == 0: # Pull out data and convert to 0-based ordering fullPtr = self.warp.cgnsgrid.surfaceptr - 1 fullPatchPtr = self.warp.cgnsgrid.surfacepatchptr - 1 fullFaceSizes = fullPtr[1:-1] - fullPtr[0:-2] # Now we loop back over the "full" versions of # thing and just take the ones that correspond # to the families we are using. curOffset = 0 for i in range(len(fullPatchNames)): # Start/end indices into fullPtr array iStart = fullPatchPtr[i] iEnd = fullPatchPtr[i + 1] # Start/end indices into the fullConn/fullPts array iStart2 = fullPtr[iStart] iEnd2 = fullPtr[iEnd] if fullPatchNames[i] in surfaceFamilies: # Keep track of the families we've actually used usedFams.add(fullPatchNames[i]) faceSizes.extend(fullFaceSizes[iStart:iEnd]) conn.extend(fullConn[iStart2:iEnd2] - curOffset) pts.extend(fullPts[iStart2 * 3 : iEnd2 * 3]) else: # If we skipped, we increment the offset curOffset += iEnd2 - iStart2 pts = np.array(pts).reshape((len(pts) // 3, 3)) # Run the common surface definition routine self.setSurfaceDefinition(pts=pts, conn=np.array(conn, "intc"), faceSizes=faceSizes) # Check if all supplied family names were actually # used. The user probably wants to know if a family # name was specified incorrectly. if self.comm.rank == 0: if usedFams < surfaceFamilies: missing = list(surfaceFamilies - usedFams) warnings.warn( "Not all specified surface families that " "were given were found the CGNS file. " "The families not found are %s." % (repr(missing)), stacklevel=2, ) if len(usedFams) == 0: raise Error( "No surfaces were selected. Check the names " "given in the 'specifiedSurface' option. The " "list of families is %s." % (repr(list(fullPatchNames))) ) elif self.fileType == "OpenFOAM": faceSizes, conn, pts = self._computeOFConn() # Run the "external" command self.setSurfaceDefinition(pts=pts, conn=conn, faceSizes=faceSizes) def _setSymmetryConditions(self): """This function determines the symmetry planes used for the computation. It has a similar structure to setInternalSurface. Symmetry planes are specified throught 'symmetryPlanes' option, which has the form 'symmetryPlanes':[[[x1,y1, z1],[dir_x1, dir_y1, dir_z1]],[[x2,y2, z2],[dir_x2, dir_y2, dir_z2]],...] Examples -------- meshOptions = {'symmetryPlanes':[[[0.,0., 0.],[0., 1., 0.]]]} mesh = USMesh(options=meshOptions,comm=gcomm) """ symmList = self.getOption("symmetryPlanes") if symmList is not None: # The user has explictly supplied symmetry planes. Use those pts = [] normals = [] for i in range(len(symmList)): pts.append(symmList[i][0]) normals.append(symmList[i][1]) else: # Otherwise generate from the geometry. planes = [] if self.fileType == "CGNS": if self.comm.rank == 0: # Do the necessary fortran preprocessing if self.warp.cgnsgrid.cgnsstructured: self.warp.processstructuredpatches() else: self.warp.processunstructuredpatches() fullConn = self.warp.cgnsgrid.surfaceconn - 1 fullPts = self.warp.cgnsgrid.surfacepoints nPatch = self.warp.cgnsgrid.getnpatch() fullPatchNames = [] for i in range(nPatch): fullPatchNames.append(self.warp.cgnsgrid.getsurf(i + 1).strip().lower()) symmFamilies = set() if self.getOption("symmetrySurfaces") is None: # Use all symmetry surfaces: for i in range(len(self.warp.cgnsgrid.surfaceissymm)): if self.warp.cgnsgrid.surfaceissymm[i]: symmFamilies.add(fullPatchNames[i].lower()) else: # The user has supplied a list of surface families for name in self.getOption("symmetrySurfaces"): symmFamilies.add(name.lower()) usedFams = set() if self.warp.cgnsgrid.cgnsstructured: if self.comm.rank == 0: # Pull out data and convert to 0-based ordering fullPatchSizes = self.warp.cgnsgrid.surfacesizes.T # Now we loop back over the "full" versions of # thing and just take the ones that correspond # to the families we are using. curNodeIndex = 0 curCellIndex = 0 curOffset = 0 for i in range(len(fullPatchNames)): curNodeSize = fullPatchSizes[i][0] * fullPatchSizes[i][1] curCellSize = (fullPatchSizes[i][0] - 1) * (fullPatchSizes[i][1] - 1) if fullPatchNames[i] in symmFamilies: # Keep track of the families we've actually used usedFams.add(fullPatchNames[i]) # Determine the average normal for these points: conn = ( fullConn[curCellIndex : curCellIndex + curCellSize * 4] - fullConn[curCellIndex] + 1 ) pts = fullPts[curNodeIndex : curNodeIndex + curNodeSize * 3] avgNorm = self.warp.averagenormal(pts, conn, 4 * np.ones(curCellSize, "intc")) planes.append([pts[0:3], avgNorm]) else: # If we skipped, we increment the offset curOffset += curNodeSize curNodeIndex += curNodeSize * 3 curCellIndex += curCellSize * 4 # end for (root proc) else: # unstructured # We won't do this in general. The issue is that # each of the elements needs to be checked # individually since one sym BC may have multiple # actual symmetry planes. This could be done, but # since plane elemination code is in python and # slow, we'll just defer this and make the user # supply the symmetry planes. raise Error( "Automatic determine of symmetry surfaces is " "not supported for unstructured CGNS. Please " "specify the symmetry planes using the " "'symmetryPlanes' option. See the _setSymmetryConditions()" " documentation string for the option prototype." ) # Check if all supplied family names were actually # used. The user probably wants to know if a family # name was specified incorrectly. if self.comm.rank == 0: if usedFams < symmFamilies: missing = list(symmFamilies - usedFams) warnings.warn( "Not all specified symm families that " "were given were found the CGNS file. " "The families not found are %s." % (repr(missing)), stacklevel=2, ) elif self.fileType in ["OpenFOAM", "PLOT3D"]: # We could probably implement this at some point, but # it is not critical raise Error( "Automatic determine of symmetry surfaces is " "not supported for OpenFOAM or PLOT3D meshes. Please " "specify the symmetry planes using the " "'symmetryPlanes' option. See the _setSymmetryConditions()" " documentation string for the option prototype." ) # Now we have a list of planes. We have to reduce them to the # set of independent planes. This is tricky since you can have # have to different normals belonging to the same physical # plane. Since we don't have that many, we just use a dumb # double loop. def checkPlane(p1, n1, p2, n2): # Determine if two planes defined by (pt, normal) are # actually the same up to a normal sign. # First check the normal...if these are not the same, # cannot be the same plane if abs(np.dot(n1, n2)) < 0.99: return False # Normals are the same direction. Check if p2 is on the # first plane up to a tolerance. d = p2 - p1 d1 = p2 - np.dot(d, n1) * n1 if np.linalg.norm(d1 - p2) / (np.linalg.norm(d) + 1e-30) > 1e-8: return False return True uniquePlanes = [] flagged = np.zeros(len(planes), "intc") for i in range(len(planes)): if not flagged[i]: uniquePlanes.append(planes[i]) curPlane = planes[i] # Loop over remainder to check: for j in range(i + 1, len(planes)): if checkPlane(curPlane[0], curPlane[1], planes[j][0], planes[j][1]): flagged[j] = 1 # Before we return, reset the point for each plane to be as # close to the origin as possible. This will help slightly # with the numerics. pts = [] normals = [] for i in range(len(uniquePlanes)): p = uniquePlanes[i][0] n = uniquePlanes[i][1] p2 = np.zeros(3) d = p2 - p pstar = p2 - np.dot(d, n) * n normals.append(n) pts.append(pstar) normals = self.comm.bcast(np.array(normals)) pts = self.comm.bcast(np.array(pts)) # Let the user know what they are: if self.comm.rank == 0: print("+-------------------- Symmetry Planes -------------------+") print("| Point Normal |") for i in range(len(pts)): print( "| (%7.3f %7.3f %7.3f) (%7.3f %7.3f %7.3f) |" % ( np.real(pts[i][0]), np.real(pts[i][1]), np.real(pts[i][2]), np.real(normals[i][0]), np.real(normals[i][1]), np.real(normals[i][2]), ) ) print("+--------------------------------------------------------+") # Now set the data into fortran. self.warp.setsymmetryplanes(pts.T, normals.T) # ========================================================================= # Local OpenFOAM Functions # ========================================================================= def _computeOFConn(self): """ The user has specified an OpenFOAM mesh. Loop through the mesh data and create the arrays necessary to initialize the warping. """ # Finally create the underlying data structure: faceSizes = [] faceConn = [] pts = [] x0 = self.OFData["x0"] faces = self.OFData["faces"] connCount = 0 for bName in self.OFData["boundaries"]: bType = self.OFData["boundaries"][bName]["type"].lower() if bType in ["patch", "wall", "slip"]: # Apparently OpenFOAM will list boundaries with zero # faces on them: nFace = len(self.OFData["boundaries"][bName]["faces"]) if nFace > 0: for iFace in self.OFData["boundaries"][bName]["faces"]: face = faces[iFace] nNodes = len(face) # Pull out the 'len(face)' nodes from x0. for j in range(nNodes): pts.append(x0[face[j]]) faceSizes.append(nNodes) faceConn.extend(range(connCount, connCount + nNodes)) connCount += nNodes pts = np.array(pts) faceConn = np.array(faceConn) return faceSizes, faceConn, pts def _readOFGrid(self, caseDir): """ Read in the mesh points and connectivity from the polyMesh directory. NOTE: To ensure a consistant starting point in DAFoam. One needs to copy points_orig to points in pyDAFoam before initializing USMesh. This function will no longer perform this copying! Parameters ---------- caseDir : str The directory containing the OpenFOAM mesh files """ # import the pyOFM repo from pyofm import PYOFM ofm = PYOFM(comm=self.comm) # Copy the reference points file to points to ensure # consistant starting point self.OFData = ofm.getFileNames(caseDir, comm=self.comm) # Read in the volume points self.OFData["x0"] = ofm.readVolumeMeshPoints() # Read the face info for the mesh self.OFData["faces"] = ofm.readFaceInfo() # Read the boundary info self.OFData["boundaries"] = ofm.readBoundaryInfo(self.OFData["faces"]) # Read the cell info for the mesh self.OFData["owner"], self.OFData["neighbour"] = ofm.readCellInfo() # Create the internal structures for volume mesh x0 = self.OFData["x0"] self.warp.createcommongrid(x0.T) ofm = None # ========================================================================= # Utility functions # ========================================================================= def _setMeshOptions(self): """Private function to set the options currently in self.options to the corresponding place in Fortran""" self.warp.gridinput.ldeffact = self.getOption("LdefFact") self.warp.gridinput.alpha = self.getOption("alpha") self.warp.gridinput.aexp = self.getOption("aExp") self.warp.gridinput.bexp = self.getOption("bExp") self.warp.gridinput.symmtol = self.getOption("symmTol") self.warp.gridinput.userotations = self.getOption("useRotations") self.warp.gridinput.zerocornerrotations = self.getOption("zeroCornerRotations") self.warp.gridinput.cornerangle = self.getOption("cornerAngle") self.warp.gridinput.errtol = self.getOption("errTol") self.warp.kd_tree.bucket_size = self.getOption("bucketSize") if self.getOption("evalMode") == "fast": self.warp.gridinput.evalmode = self.warp.gridinput.eval_fast elif self.getOption("evalMode") == "exact": self.warp.gridinput.evalmode = self.warp.gridinput.eval_exact def _processFortranStringArray(self, strArray): """Getting arrays of strings out of Fortran can be kinda nasty. This takes the array and returns a nice python list of strings""" shp = strArray.shape arr = strArray.reshape((shp[1], shp[0]), order="F") tmp = [] for i in range(arr.shape[1]): tmp.append("") for j in range(arr.shape[0]): tmp[-1] += arr[j, i].decode() tmp[-1] = tmp[-1].strip().lower() return tmp def __del__(self): """Release all the mesh warping memory. This should be called automatically when the object is garbage collected.""" self.warp.releasememory()