# External modules
from baseclasses.utils import Error
import numpy as np
# Local modules
from .projection import projectNode
[docs]
def createMidsurfaceMesh(
surf, surfFormat, leList, teList, nSpan, nChord, liftIndex, chordCosSpacing=0.0, spanOffset=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). 1/2/3 for x/y/z-axis. This is used to
determine the projection direction.
chord_cos_spacing : float, optional
Cosine spacing factor for the chordwise points. 0.0 is uniform, 1.0 is cosine spacing, by default 0.0.
spanOffset : float, optional
Point projection often fails if attempted at the exact wing root and tip, this offset shifts the points at the
wing root/tip 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 at each end of the wing. By default, 1e-5.
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, and the third is
the x, y, z. The spanwise ordering of the points will match the order of the leList and teList.
Notes
-----
If you plan to use the mesh generated by this function with OpenAeroStruct, and your leList/teList are ordered
from root to tip, then 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 end points inward slightly to avoid projection errors, we want to shift the points at the end of the LE
# along the LE direction and the point at the ends of the TE along the TE direction
shiftVec1 = np.outer((1 - chordwiseSpacing), (leCoords[1] - leCoords[0])) + np.outer(
chordwiseSpacing, (teCoords[1] - teCoords[0])
)
meshCoords[:, 0, :] += shiftVec1 * spanOffset
shiftVec2 = np.outer((1 - chordwiseSpacing), (leCoords[-1] - leCoords[-2])) + np.outer(
chordwiseSpacing, (teCoords[-1] - teCoords[-2])
)
meshCoords[:, -1, :] += shiftVec2 * spanOffset
# 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, :] -= shiftVec1 * spanOffset
meshCoords[:, -1, :] -= shiftVec2 * spanOffset
return meshCoords