Add helper function for creating midsurface (a.k.a VLM) meshes

doc/geo_utils.rst

@@ -5,6 +5,11 @@ File IO
 .. automodule:: pygeo.geo_utils.file_io
     :members:
 
+Mesh Generation
+~~~~~~~~~~~~~~
+.. automodule:: pygeo.geo_utils.mesh_generation
+    :members:
+
 FFD Generation
 ~~~~~~~~~~~~~~
 .. automodule:: pygeo.geo_utils.ffd_generation

pygeo/geo_utils/__init__.py

@@ -8,6 +8,7 @@
 from .dcel import *  # noqa: F401, F403
 from .file_io import *  # noqa: F401, F403
 from .ffd_generation import *  # noqa: F401, F403
+from .mesh_generation import *  # noqa: F401, F403
 from .index_position import *  # noqa: F401, F403
 from .knotvector import *  # noqa: F401, F403
 from .misc import *  # noqa: F401, F403

pygeo/geo_utils/mesh_generation.py

@@ -0,0 +1,139 @@
+import numpy as np
+from baseclasses.utils import Error
+
+from .projection import projectNode
+
+
+def createMidsurfaceMesh(
+    surf, surfFormat, leList, teList, nSpan, nChord, liftIndex, chordCosSpacing=0.0, rootOffset=1e-5
+):
+    """Create a cambered midsurface mesh (e.g to be used with a VLM) by projecting points to the upper and lower wing
+    surfaces
+
+    Parameters
+    ----------
+    surf : pyGeo object or list or str
+        The surface around which the FFD will be created.
+        See the documentation for :meth:`pygeo.constraints.DVCon.DVConstraints.setSurface` for details.
+
+    surfFormat : str
+        The surface format.
+        See the documentation for :meth:`pygeo.constraints.DVCon.DVConstraints.setSurface` for details.
+
+    leList : list or array
+        List or array of points (of size Nx3 where N is at least 2) defining the 'leading edge'.
+
+    teList : list or array
+        Same as leList but for the trailing edge.
+
+    nSpan : int or list of int
+        Number of spanwise points in the mesh.
+        Use a list of length N-1 to specify the number for each segment defined by leList and teList
+        and to precisely match intermediate locations.
+
+    nChord : int
+        Number of chordwise points in the mesh.
+
+    liftIndex : int
+        Index specifying which direction lift is in (same as the ADflow option).
+        Either 2 for the y-axis or 3 for the z-axis.
+        This is used to determine the wing's spanwise direction.
+
+    chord_cos_spacing : float
+        Cosine spacing factor for the chordwise points. 0.0 is uniform, 1.0 is cosine spacing.
+
+    rootOffset : float
+        Point projection often fails if attempted at the exact wing root, this offset shifts the points at the wing root
+        inward slightly to avoid this, before shifting them back to the correct location after projection. The value of
+        the offset is how far to shift the points as a fraction of the distance between the first and second
+        leading/trailing edge points
+
+    Returns
+    -------
+    meshCoords : ndarray
+        3D array of size (nChord, nSpanTotal, 3) containing the coordinates of the mesh points.
+        The first index is the chordwise direction (LE to TE), the second is the spanwise direction (Root to tip), and the third is the x, y, z
+
+    Notes
+    -----
+        If you plan to use the mesh generated by this function with OpenAeroStruct, you will need to reverse the order
+        of the spanwise points so that they go from the tip to the root. This can be done with
+        `meshCoords = np.flip(meshCoords, 1)`
+    """
+
+    # Import inside this function to avoid circular imports
+    # First party modules
+    from pygeo import DVConstraints
+
+    # Set the triangulated surface in DVCon
+    DVCon = DVConstraints()
+    DVCon.setSurface(surf, surfFormat=surfFormat)
+
+    p0, p1, p2 = DVCon._getSurfaceVertices(surfaceName="default")
+
+    # Create 2D grid of points between leading and trailing edge
+    cosineSpacing = 0.5 * (1 - np.cos(np.linspace(0, np.pi, nChord)))
+    uniformSpacing = np.linspace(0, 1, nChord)
+    chordwiseSpacing = chordCosSpacing * cosineSpacing + (1 - chordCosSpacing) * uniformSpacing
+    if isinstance(nSpan, int):
+        nSpan = [nSpan] * (len(leList) - 1)
+    else:
+        if len(nSpan) != len(leList) - 1:
+            raise Error(
+                f"nSpan must be an integer or a list of length len(leList) - 1. Got {len(nSpan)=} and {len(leList)=}."
+            )
+    nSpanTotal = np.sum(nSpan) - (len(nSpan) - 1)
+    meshCoords = np.zeros((nChord, nSpanTotal, 3))
+
+    # Create a flat mesh by linearly interpolating between the leading and trailing edge
+    nSegments = len(nSpan)
+    spanEnd = 1
+    for ii in range(nSegments):
+        nSpanSegment = nSpan[ii]
+        spanStart = spanEnd - 1
+        spanEnd = spanStart + nSpanSegment
+
+        leCoords = np.linspace(leList[ii], leList[ii + 1], nSpanSegment)
+        teCoords = np.linspace(teList[ii], teList[ii + 1], nSpanSegment)
+        for jj, chordFraction in enumerate(chordwiseSpacing):
+            meshCoords[jj, spanStart:spanEnd, :] = (1 - chordFraction) * leCoords + chordFraction * teCoords
+
+    # Shift the root points inward slightly to avoid projection errors, we want to shift the point at the root LE along
+    # the LE direction and the point at the root TE along the TE direction
+    rootShiftVec = np.outer((1 - chordwiseSpacing), (leCoords[1] - leCoords[0])) + np.outer(
+        chordwiseSpacing, (teCoords[1] - teCoords[0])
+    )
+    meshCoords[:, 0, :] += rootShiftVec * rootOffset
+
+    # Now loop through all but the leading and trailing edge points and project to the surface, then take the midpoint
+    # and set that as the midsurface
+    for jj in range(nSpanTotal):
+        # Project in the component of the lift direction normal to the chord line between the leading and trailing edge
+        projectDir = np.zeros(3)
+        projectDir[liftIndex - 1] = 1
+        chordDir = meshCoords[0, jj] - meshCoords[-1, jj]
+        chordDir /= np.linalg.norm(chordDir)
+        projectDir -= np.dot(projectDir, chordDir) * chordDir
+        projectDir /= np.linalg.norm(projectDir)
+        for ii in range(1, nChord - 1):
+            up, down, fail = projectNode(meshCoords[ii, jj], projectDir, p0, p1 - p0, p2 - p0)
+
+            if fail in [0, -1]:
+                meshCoords[ii, jj] = 0.5 * (up + down)
+            else:
+                raise Error(
+                    "There was an error projecting a node at (%f, %f, %f) with normal (%f, %f, %f)."
+                    % (
+                        meshCoords[ii, jj, 0],
+                        meshCoords[ii, jj, 1],
+                        meshCoords[ii, jj, 2],
+                        projectDir[0],
+                        projectDir[1],
+                        projectDir[2],
+                    )
+                )
+
+    # Shift the root points back
+    meshCoords[:, 0, :] -= rootShiftVec * rootOffset
+
+    return meshCoords

tests/reg_tests/test_ffdGeneration.py

@@ -8,14 +8,25 @@
 
 # First party modules
 from pygeo import DVGeometry
-from pygeo.geo_utils import createFittedWingFFD, write_wing_FFD_file
+from pygeo.geo_utils import createFittedWingFFD, write_wing_FFD_file, createMidsurfaceMesh
 
 baseDir = os.path.dirname(os.path.abspath(__file__))
 
 
 class TestFFDGeneration(unittest.TestCase):
     N_PROCS = 1
 
+    def setUp(self):
+        # Get the surface definition from the STL file
+        meshFile = os.path.join(baseDir, "../../input_files/c172.stl")
+        stlMesh = Mesh.from_file(meshFile)
+        self.p0 = stlMesh.vectors[:, 0, :] * 1e-3
+        self.p1 = stlMesh.vectors[:, 1, :] * 1e-3
+        self.p2 = stlMesh.vectors[:, 2, :] * 1e-3
+        self.surf = [self.p0, self.p1, self.p2]
+        self.surfFormat = "point-point"
+        self.liftIndex = 2
+
     def test_box_ffd(self, train=False, refDeriv=False):
         # Write duplicate of outerBoxFFD
         axes = ["i", "k", "j"]
@@ -50,24 +61,25 @@ def test_c172_fitted(self):
         leList = np.array([[0.0, 0.0, 0.1], [10.0, 0.0, 2500.0], [160.0, 0.0, 5280.0]]) * 1e-3
         teList = np.array([[1600.0, 0.0, 0.1], [1650.0, 0.0, 2500.0], [1320.0, 0.0, 5280.0]]) * 1e-3
 
-        # Get the surface definition from the STL file
-        meshFile = os.path.join(baseDir, "../../input_files/c172.stl")
-        stlMesh = Mesh.from_file(meshFile)
-        p0 = stlMesh.vectors[:, 0, :] * 1e-3
-        p1 = stlMesh.vectors[:, 1, :] * 1e-3
-        p2 = stlMesh.vectors[:, 2, :] * 1e-3
-        surf = [p0, p1, p2]
-        surfFormat = "point-point"
-
         # Set the other FFD generation inputs
         outFile = "wing_ffd.xyz"
         nSpan = [4, 4]
         nChord = 8
         relMargins = [0.02, 0.01, 0.01]
         absMargins = [0.04, 0.01, 0.02]
-        liftIndex = 2
 
-        createFittedWingFFD(surf, surfFormat, outFile, leList, teList, nSpan, nChord, absMargins, relMargins, liftIndex)
+        createFittedWingFFD(
+            self.surf,
+            self.surfFormat,
+            outFile,
+            leList,
+            teList,
+            nSpan,
+            nChord,
+            absMargins,
+            relMargins,
+            self.liftIndex,
+        )
 
         # Check that the generated FFD file matches the reference
         referenceFFD = DVGeometry(os.path.join(baseDir, "../../input_files/c172_fitted.xyz"))
@@ -77,13 +89,42 @@ def test_c172_fitted(self):
         # Check that the embedding works
         # This is not an actual test because no errors are raised if the projection does not work
         DVGeo = DVGeometry(outFile)
-        DVGeo.addPointSet(p0, "wing_p0")
-        DVGeo.addPointSet(p1, "wing_p1")
-        DVGeo.addPointSet(p2, "wing_p2")
+        DVGeo.addPointSet(self.p0, "wing_p0")
+        DVGeo.addPointSet(self.p1, "wing_p1")
+        DVGeo.addPointSet(self.p2, "wing_p2")
 
         # Delete the generated FFD file
         os.remove(outFile)
 
+    def test_c172_midsurface_mesh(self):
+        # These LE and TE coordinates account for the twist in the wing unlike those used in test_c172_fitted
+        x = [0, 0, 0.125 * 1.18]
+        y = [0, 0, 0]
+        z = [0.0, 2.5, 10.58 / 2]
+        chord = [1.67, 1.67, 1.18]
+        rot_z = [0, 0, 2]
+        leList = np.array([x, y, z]).T
+        teList = np.array([x, y, z]).T
+        for ii in range(len(chord)):
+            teList[ii] = leList[ii] + np.array(
+                [chord[ii] * np.cos(np.deg2rad(rot_z[ii])), chord[ii] * np.sin(np.deg2rad(rot_z[ii])), 0]
+            )
+        nSpan = 21
+        nChord = 21
+        mesh = createMidsurfaceMesh(
+            self.surf, self.surfFormat, leList, teList, nSpan, nChord, self.liftIndex, chordCosSpacing=0.75
+        )
+
+        # Check that the generated mesh matches the reference
+        flattenedMesh = mesh.reshape((-1, 3))
+        refMesh = np.loadtxt(os.path.join(baseDir, "../../input_files/c172-midsurfaceMesh.txt"))
+        np.testing.assert_allclose(flattenedMesh, refMesh, rtol=1e-13)
+
+        # Check that the generated mesh fits inside the FFD
+        # This is not an actual test because no errors are raised if the projection does not work
+        DVGeo = DVGeometry(os.path.join(baseDir, "../../input_files/c172_fitted.xyz"))
+        DVGeo.addPointSet(flattenedMesh, "midsurfaceMesh")
+
 
 if __name__ == "__main__":
     unittest.main()