Skip to content

Commit

Permalink
Added coordinate transformations and a few fixes (#157)
Browse files Browse the repository at this point in the history 
* checkpoint

* fix typo in docstring

* several fixes related to use cases with a symmetry plane

* fixed some integer divisions

* added the optional coordinate transformation routines

* bug fix in the section local dv class

* editing the symmetry fix a bit

* clarified comments

* too many typo fixes. added test for the coord xfer feature

* format and lint. changed "dir" argument to "mode"

* added pyspline solver test. this should fail until the pyspline pr #47 is merged. Updated the input file archive to include sabet's fuselage FFD

* flake fix

* minor tolerance adjustment on ffd generation test. this fails on my laptop with 1e-15 rel tolerance. the error is at 1.1e-15. Moved the tolerance to 1e-14, which should be plenty accurate

* addressing todos

* minor rewording

* addressing comments
  • Loading branch information
@anilyil
anilyil committed Oct 5, 2022
1 parent 4641625 commit 022e6b2
Show file tree
Hide file tree
Showing 5 changed files with 301 additions and 16 deletions.
190 changes: 177 additions & 13 deletions pygeo/parameterization/DVGeo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,9 @@ def __init__(self, fileName, *args, isComplex=False, child=False, faceFreeze=Non
self.JT = {}
self.nPts = {}

# dictionary to save any coordinate transformations we are given
self.coord_xfer = {}

# Derivatives of Xref and Coef provided by the parent to the
# children
self.dXrefdXdvg = None
Expand Down Expand Up @@ -428,10 +431,13 @@ def addRefAxis(
for volume in volumes:
volumesSymm.append(volume + self.FFD.nVol / 2)

curveSymm = copy.deepcopy(curve)
curveSymm.reverse()
for _coef in curveSymm.coef:
curveSymm.coef[:, index] = -curveSymm.coef[:, index]
# We want to create a curve that is symmetric of the current one
symm_curve_X = curve.X.copy()

# flip the coefs
symm_curve_X[:, index] = -symm_curve_X[:, index]
curveSymm = Curve(k=curve.k, X=symm_curve_X)

self.axis[name] = {
"curve": curve,
"volumes": volumes,
Expand Down Expand Up @@ -600,9 +606,39 @@ def addRefAxis(
# Add the raySize multiplication factor for this axis
self.axis[name]["raySize"] = raySize

# do the same for the other half if we have a symmetry plane
if self.FFD.symmPlane is not None:
# we need to figure out the correct indices to ignore for the mirrored FFDs

# first get the matching indices between the current and mirroring FFDs.
# we want to include the the nodes on the symmetry plane.
# these will appear as the same indices on FFDs on both sides
indSetA, indSetB = self.getSymmetricCoefList(getSymmPlane=True)

# loop over the inds_to_ignore list and find the corresponding symmetries
ignoreIndSymm = []
for ind in ignoreInd:
try:
tmp = indSetA.index(ind)
except ValueError:
raise Error(
f"""The index {ind} is not in indSetA. This is likely due to a weird
issue caused by the point reduction routines during initialization.
Reduce the offset of the FFD control points from the symmetry plane
to avoid it. The max deviation from the symmetry plane needs to be
less than around 1e-5 if rest of the default tolerances in pygeo is used."""
)
ind_mirror = indSetB[tmp]
ignoreIndSymm.append(ind_mirror)

self.axis[name + "Symm"]["ignoreInd"] = ignoreIndSymm

# we just take the same raySize as the original curve
self.axis[name + "Symm"]["raySize"] = raySize

return nAxis

def addPointSet(self, points, ptName, origConfig=True, **kwargs):
def addPointSet(self, points, ptName, origConfig=True, coord_xfer=None, **kwargs):
"""
Add a set of coordinates to DVGeometry
Expand All @@ -624,6 +660,72 @@ def addPointSet(self, points, ptName, origConfig=True, **kwargs):
undeformed or deformed configuration. This should almost
always be True except in circumstances when the user knows
exactly what they are doing.
coord_xfer : function
A callback function that performs a coordinate transformation
between the DVGeo reference frame and any other reference
frame. The DVGeo object uses this routine to apply the coordinate
transformation in "forward" and "reverse" directions to go between
the two reference frames. Derivatives are also adjusted since they
are vectors coming into DVGeo (in the reverse AD mode)
and need to be rotated. We have a callback function here that lets
the user to do whatever they want with the coordinate transformation.
The function must have the first positional argument as the array that is
(npt, 3) and the two keyword arguments that must be available are "mode"
("fwd" or "bwd") and "apply_displacement" (True or False). This function
can then be passed to DVGeo through something like ADflow, where the
set DVGeo call can be modified as:
CFDSolver.setDVGeo(DVGeo, pointSetKwargs={"coord_xfer": coord_xfer})
An example function is as follows:
.. code-block:: python
def coord_xfer(coords, mode="fwd", apply_displacement=True, **kwargs):
# given the (npt by 3) array "coords" apply the coordinate transformation.
# The "fwd" mode implies we go from DVGeo reference frame to the
# application, e.g. CFD, the "bwd" mode is the opposite;
# goes from the CFD reference frame back to the DVGeo reference frame.
# the apply_displacement flag needs to be correctly implemented
# by the user; the derivatives are also passed through this routine
# and they only need to be rotated when going between reference frames,
# and they should NOT be displaced.
# In summary, all the displacements MUST be within the if apply_displacement == True
# checks, otherwise the derivatives will be wrong.
# Example transfer: The CFD mesh
# is rotated about the x-axis by 90 degrees with the right hand rule
# and moved 5 units below (in z) the DVGeo reference.
# Note that the order of these operations is important.
# a different rotation matrix can be created during the creation of
# this function. This is a simple rotation about x-axis.
# Multiple rotation matrices can be used; the user is completely free
# with whatever transformations they want to apply here.
rot_mat = np.array([
[1, 0, 0],
[0, 0, -1],
[0, 1, 0],
])
if mode == "fwd":
# apply the rotation first
coords_new = np.dot(coords, rot_mat)
# then the translation
if apply_displacement:
coords_new[:, 2] -= 5
elif mode == "bwd":
# apply the operations in reverse
coords_new = coords.copy()
if apply_displacement:
coords_new[:, 2] += 5
# and the rotation. note the rotation matrix is transposed
# for switching the direction of rotation
coords_new = np.dot(coords_new, rot_mat.T)
return coords_new
"""

Expand All @@ -636,6 +738,17 @@ def addPointSet(self, points, ptName, origConfig=True, **kwargs):
self.nPts[ptName] = None

points = np.array(points).real.astype("d")

# save the coordinate transformation info
if coord_xfer is not None:
self.coord_xfer[ptName] = coord_xfer

# Also apply the first coordinate transformation while adding this ptset.
# The child FFDs only interact with their parent FFD, and therefore,
# do not need to access the coordinate transformation routine; i.e.
# all transformations are applied once during the highest level DVGeo object.
points = self.coord_xfer[ptName](points, mode="bwd", apply_displacement=True)

self.points[ptName] = points

# Ensure we project into the undeformed geometry
Expand Down Expand Up @@ -834,7 +947,7 @@ def addLocalDV(
for vol in volList:
volListTmp.append(vol)
for vol in volList:
volListTmp.append(vol + self.FFD.nVol / 2)
volListTmp.append(vol + self.FFD.nVol // 2)
volList = volListTmp

volList = np.atleast_1d(volList).astype("int")
Expand Down Expand Up @@ -952,7 +1065,7 @@ def addSpanwiseLocalDV(
for vol in volList:
volListTmp.append(vol)
for vol in volList:
volListTmp.append(vol + self.FFD.nVol / 2)
volListTmp.append(vol + self.FFD.nVol // 2)
volList = volListTmp

volList = np.atleast_1d(volList).astype("int")
Expand Down Expand Up @@ -1186,7 +1299,7 @@ class in geo_utils. Using pointSelect discards everything in volList.
for vol in volList:
volListTmp.append(vol)
for vol in volList:
volListTmp.append(vol + self.FFD.nVol / 2)
volListTmp.append(vol + self.FFD.nVol // 2)
volList = volListTmp

volList = np.atleast_1d(volList).astype("int")
Expand Down Expand Up @@ -1284,7 +1397,7 @@ def addCompositeDV(self, dvName, ptSetName=None, u=None, scale=None):
self.DVComposite = geoDVComposite(dvName, values, NDV, u, scale=scale, s=s)
self.useComposite = True

def getSymmetricCoefList(self, volList=None, pointSelect=None, tol=1e-8):
def getSymmetricCoefList(self, volList=None, pointSelect=None, tol=1e-8, getSymmPlane=False):
"""
Determine the pairs of coefs that need to be constrained for symmetry.
Expand All @@ -1300,6 +1413,13 @@ def getSymmetricCoefList(self, volList=None, pointSelect=None, tol=1e-8):
tol : float
Tolerance for ignoring nodes around the symmetry plane. These should be
merged by the network/connectivity anyway
getSymmPlane : bool
If this flag is set to True, we also return the points on the symmetry plane
for all volumes. e.g. a reduced point on the symmetry plane with the same
indices on both volumes will show up as the same value in both arrays. This
is useful when determining the indices to ignore when adding pointsets. The
default behavior will not include the points exactly on the symmetry plane.
this is more useful for adding them as linear constraints
Returns
-------
Expand Down Expand Up @@ -1335,7 +1455,7 @@ def getSymmetricCoefList(self, volList=None, pointSelect=None, tol=1e-8):
for vol in volList:
volListTmp.append(vol)
for vol in volList:
volListTmp.append(vol + self.FFD.nVol / 2)
volListTmp.append(vol + self.FFD.nVol // 2)
volList = volListTmp

volList = np.atleast_1d(volList).astype("int")
Expand Down Expand Up @@ -1368,11 +1488,28 @@ def getSymmetricCoefList(self, volList=None, pointSelect=None, tol=1e-8):
# Now find any matching nodes within tol. there should be 2 and
# only 2 if the mesh is symmetric
Ind = tree.query_ball_point(pt, tol) # should this be a separate tol
if not (len(Ind) == 2):
raise Error("more than 2 coefs found that match pt")
if len(Ind) == 2:
# check which point is on which side
if pts[Ind[0], index] > 0:
# first one is on the primary side
indSetA.append(Ind[0])
indSetB.append(Ind[1])
else:
# flip the order
indSetA.append(Ind[1])
indSetB.append(Ind[0])
else:
raise Error("more than 2 coefs found that match pt")

elif (abs(pt[index]) < tol) and getSymmPlane:
# this point is on the symmetry plane
# if everything went right so far, this should return only one point
Ind = tree.query_ball_point(pt, tol)
if len(Ind) == 1:
indSetA.append(Ind[0])
indSetB.append(Ind[1])
indSetB.append(Ind[0])
else:
raise Error("more than 1 coefs found that match pt on symmetry plane")

return indSetA, indSetB

Expand Down Expand Up @@ -1826,6 +1963,11 @@ def update(self, ptSetName, childDelta=True, config=None):
if self.isChild and childDelta:
return Xfinal - Xstart
else:
# we only check if we need to apply the coordinate transformation
# and move the pointset to the reference frame of the application,
# if this is the last pygeo in the chain
if ptSetName in self.coord_xfer:
Xfinal = self.coord_xfer[ptSetName](Xfinal, mode="fwd", apply_displacement=True)
return Xfinal

def applyToChild(self, iChild):
Expand Down Expand Up @@ -2061,6 +2203,14 @@ def totalSensitivity(self, dIdpt, ptSetName, comm=None, config=None):
dIdpt = np.array([dIdpt])
N = dIdpt.shape[0]

# apply the coordinate transformation on dIdpt if this pointset has it.
if ptSetName in self.coord_xfer:
# loop over functions
for ifunc in range(N):
# its important to remember that dIdpt are vector-like values,
# so we don't apply the transformations and only the rotations!
dIdpt[ifunc] = self.coord_xfer[ptSetName](dIdpt[ifunc], mode="bwd", apply_displacement=False)

# generate the total Jacobian self.JT
self.computeTotalJacobian(ptSetName, config=config)

Expand Down Expand Up @@ -2161,6 +2311,13 @@ def totalSensitivityProd(self, vec, ptSetName, config=None):
else:
xsdot = self.JT[ptSetName].T.dot(newvec)
xsdot.reshape(len(xsdot) // 3, 3)

# check if we have a coordinate transformation on this ptset
if ptSetName in self.coord_xfer:
# its important to remember that dIdpt are vector-like values,
# so we don't apply the transformations and only the rotations!
xsdot = self.coord_xfer[ptSetName](xsdot, mode="fwd", apply_displacement=False)

# Maybe this should be:
# xsdot = xsdot.reshape(len(xsdot)//3, 3)

Expand Down Expand Up @@ -2216,6 +2373,13 @@ def totalSensitivityTransProd(self, vec, ptSetName, config=None):
if self.JT[ptSetName] is None:
xsdot = np.zeros((0, 3))
else:

# check if we have a coordinate transformation on this ptset
if ptSetName in self.coord_xfer:
# its important to remember that dIdpt are vector-like values,
# so we don't apply the transformations and only the rotations!
vec = self.coord_xfer[ptSetName](vec, mode="bwd", apply_displacement=False)

xsdot = self.JT[ptSetName].dot(np.ravel(vec))

# Pack result into dictionary
Expand Down
2 changes: 1 addition & 1 deletion pygeo/parameterization/designVars.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ def __init__(self, name, lower, upper, scale, axis, coefListIn, mask, config, se
self.coefList.append(coefListIn[i])

nVal = len(self.coefList)
super().__init__(name=name, value=np.zeros(nVal, "D"), nVal=nVal, lower=None, upper=None, scale=scale)
super().__init__(name=name, value=np.zeros(nVal, "D"), nVal=nVal, lower=lower, upper=upper, scale=scale)

self.config = config

Expand Down
2 changes: 1 addition & 1 deletion pygeo/pyBlock.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ class pyBlock:
uniform (and symmetric) knot vectors are assumed
everywhere. This ensures a seamless FFD.
symPlane : {"x", "y", or "z"}
symmPlane : {"x", "y", or "z"}
if a coordinate direciton is provided, the code will duplicate
the FFD in the mirroring direction.
Expand Down

0 comments on commit 022e6b2

Please sign in to comment.