Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bug fix for non-aligned FFD blocks for rotType=0 and generalization of xFraction #67

Merged
merged 23 commits into from
Mar 24, 2021
Merged

Bug fix for non-aligned FFD blocks for rotType=0 and generalization of xFraction #67

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
fe0f6ee
added new vars for rotType0 scaling handling
marcomangano Feb 8, 2021
027d285
applying axis scaling to rot0 - Mads way
marcomangano Feb 8, 2021
f49f010
extended nonaligned scaling to all rotTypes
marcomangano Feb 8, 2021
86cf8d1
added WARNING - nonzero rotType have inaccurate sensitivities
marcomangano Feb 8, 2021
24d56c0
bugfix on single axis scaling
marcomangano Feb 11, 2021
36ed707
generalized axis nodes xyz location
marcomangano Feb 15, 2021
42b7d61
implemented Mads-Sandy fix for rotAxisnodes when ffds are not aligned…
marcomangano Feb 15, 2021
0f9890f
Merge remote-tracking branch 'origin' into rot0-fix
marcomangano Feb 24, 2021
1e0ee66
added test for generalized xyz fraction
marcomangano Feb 25, 2021
f96758c
added functions to commonUtils
marcomangano Feb 25, 2021
e375a8f
WIP: added test for non-aligned FFDs - FD-trained test is failing
marcomangano Feb 25, 2021
ff2fa29
updated test 24 ref
marcomangano Feb 25, 2021
a0dde96
Merge remote-tracking branch 'origin' into rot0-fix
marcomangano Feb 26, 2021
5d5982a
Merge remote-tracking branch 'origin' into rot0-fix
marcomangano Mar 1, 2021
f063378
added tests explanation
marcomangano Mar 3, 2021
94812eb
renamed vars to singular
marcomangano Mar 3, 2021
b2c4a8a
cleaner coordinates array reshaping - thx Josh
marcomangano Mar 3, 2021
0df4ddc
less intrusive xyz fraction behavior
marcomangano Mar 3, 2021
7fda92b
bugfix on a check that did not make sense
marcomangano Mar 3, 2021
991d5d7
updated ang checks
marcomangano Mar 16, 2021
71798d4
limited rotation to rotType=0 only
marcomangano Mar 16, 2021
a598be5
Merge branch 'master' into rot0-fix
ewu63 Mar 16, 2021
cd3db0d
Merge branch 'master' into rot0-fix
marcomangano Mar 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
131 changes: 115 additions & 16 deletions pygeo/DVGeometry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,8 +196,9 @@ def __init__(self, fileName, complex=False, child=False, faceFreeze=None, name=N
tmp[ind] = True
self.masks = tmp

def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
def addRefAxis(self, name, curve=None, xFraction=None, yFraction=None, zFraction=None, volumes=None,
rotType=5, axis='x', alignIndex=None, rotAxisVar=None,
rot0ang=None, rot0axis=[1, 0, 0],
xFractionOrder=2, includeVols=[], ignoreInd=[],
raySize=1.5):
"""
Expand Down Expand Up @@ -268,7 +269,20 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
variable which should be used to compute the orientation of the theta
rotation.

xFractionOrder : int
rot0ang: float
If rotType == 0, defines the offset angle of the (child) FFD with respect
to the main system of reference. This is necessary to use the scaling functions
`scale_x`, `scale_y`, and `scale_z` with rotType == 0. The axis of rotation is
defined by `rot0axis`.

rot0axis: list
If rotType == 0, defines the rotation axis for the rotation offset of the
FFD grid given by `rot0ang`. The variable has to be a list of 3 floats
defining the [x,y,z] components of the axis direction.
This is necessary to use the scaling functions `scale_x`, `scale_y`,
and `scale_z` with rotType == 0.

xFractionOrder : int (NOT USED?)
Order of spline used for refaxis curve.

includeVols : list
Expand Down Expand Up @@ -325,7 +339,7 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
if volumes is None:
volumes = numpy.arange(self.FFD.nVol)
self.axis[name] = {'curve':curve, 'volumes':volumes,
'rotType':rotType, 'axis':axis}
'rotType':rotType, 'axis':axis, 'rot0ang':rot0ang, 'rot0axis':rot0axis}

else:
# get the direction of the symmetry plane
Expand All @@ -349,11 +363,11 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
for coef in curveSymm.coef:
curveSymm.coef[:,index]=-curveSymm.coef[:,index]
self.axis[name] = {'curve':curve, 'volumes':volumes,
'rotType':rotType, 'axis':axis}
'rotType':rotType, 'axis':axis,'rot0ang':rot0ang, 'rot0axis':rot0axis}
self.axis[name+'Symm'] = {'curve':curveSymm, 'volumes':volumesSymm,
'rotType':rotType, 'axis':axis}
'rotType':rotType, 'axis':axis, 'rot0ang':rot0ang, 'rot0axis':rot0axis}
nAxis = len(curve.coef)
elif xFraction is not None:
elif xFraction or yFraction or zFraction:
# Some assumptions
# - FFD should be a close approximation of geometry surface so that
# xFraction roughly corresponds to airfoil LE, TE, or 1/4 chord
Expand All @@ -363,6 +377,14 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
# included
# - 'x' is streamwise direction

# Default to "mean" ref axis location along non-user specified direction
if xFraction is None:
xFraction = 0.5
if yFraction is None:
yFraction = 0.5
if zFraction is None:
zFraction = 0.5

# This is the block direction along which the reference axis will lie
# alignIndex = 'k'
if alignIndex is None:
Expand Down Expand Up @@ -416,16 +438,51 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
# Loop through sections and compute node location
place = 0
for j, vol in enumerate(volOrd):
# sectionArr: indices of FFD points grouped by section
sectionArr = numpy.rollaxis(lIndex[vol], alignIndex, 0)
skip = 0
if j > 0:
skip = 1
for i in range(nSections[j]):
LE = numpy.min(self.FFD.coef[sectionArr[i+skip,:,:],0])
TE = numpy.max(self.FFD.coef[sectionArr[i+skip,:,:],0])
refaxisNodes[place+i,0] = xFraction*(TE - LE) + LE
refaxisNodes[place+i,1] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],1])
refaxisNodes[place+i,2] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],2])
# getting all the section control points coordinates
pts_tens = self.FFD.coef[sectionArr[i + skip, :, :], :] # shape=(A,B,3)
marcomangano marked this conversation as resolved.
Show resolved Hide resolved
# reshaping into vector to allow rotation (if needed)
pts_vec = numpy.copy(pts_tens.reshape(numpy.shape(pts_tens)[0] * numpy.shape(pts_tens)[1], 3)) # shape=(A*B,3)
marcomangano marked this conversation as resolved.
Show resolved Hide resolved

if rotType == 0 and rot0ang:
# rotating the FFD to be aligned with main axes
for ct_ in range(numpy.shape(pts_vec)[0]):
# here we loop over the pts_vec, rotate them and insert them inplace in pts_vec again
p_ = numpy.copy(pts_vec[ct_ , :])
joanibal marked this conversation as resolved.
Show resolved Hide resolved
p_rot = geo_utils.rotVbyW(p_, rot0axis, numpy.pi / 180 * (rot0ang))
pts_vec[ct_ , :] = p_rot

# getting the bounds of the FFD section
x_min = numpy.min(pts_vec[:, 0])
x_max = numpy.max(pts_vec[:, 0])
y_min = numpy.min(pts_vec[:, 1])
y_max = numpy.max(pts_vec[:, 1])
z_min = numpy.min(pts_vec[:, 2])
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there is a possibility that this can lead to a reference axis that lies outside the volume

for example: if you have a typical ffd with the k direction as the ref axis and you rotate it in the yaw axis (so the wing is swept but the inboard edge isn't aligned with flow direction, then you can end up specifying curve ends that lie outside the FFD volume

though to be fair if you really really cambered our existing approach you could end up with a ref axis that lies outside the ffd volume

not sure if this matters or not but figured I'd bring it up

a bigger problem is that I don't think this is backward compatible with the existing default behavior. Right now the y and z frac depend on the xfrac automatically, with the new setup they're fixed at 0.5. I am surprised that this isn't failing reg tests (do we not use the xfrac option anyplace? or possibly our reg tests use rectangular volumes without pretwist or precamber)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The example you bring up is interesting... I am actually wondering if our current code could incur in such a problem.
I assume that for more complex FFD blocks you might want to "manually" feed the Ref Axis curve rather than relying on this approach.
However, at least with the current fix, the FFD is rotated back in the "correct" alignment before you define the nodes, similarly to what is done for the scaling fix, so this "yawed" FFD box you bring up should not be an issue. The correct set-up depends on the user, once more.
Anyway, this is not a good excuse for a sloppy fix, and we should take care of any corner case that comes up to mind. Let me sleep on this and get back to you!

About the backwards incompatibility, I am not sure it is actually such a huge problem here.
In the current code the y/z coordinates are given straight up by:

refaxisNodes[place+i,1] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],1])
refaxisNodes[place+i,2] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],2])

So it's a mean over the entire section points rather than the average of min/max. I agree it's a bit brute force, and I am open to improve it.
For a parallelepiped FFD the behaviour is identical, but I see how the behaviour for a cambered airfoil might be a bit different.
Maybe I am missing some major red flag here. Again, let me sleep on this!

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm also worried about the changing the default behavior here.
Can we changing it so if only xFraction is given it works the same as before?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I am envisioning something like getting rid of the defaults Fraction=0.5 and instead having an approach similar to before extended to y and z too.
For example:

if xFraction:
       xMin = numpy.min(self.FFD.coef[sectionArr[i+skip,:,:],0])
       xMax = numpy.max(self.FFD.coef[sectionArr[i+skip,:,:],0])
       refaxisNodes[place+i,0] = xFraction*(xMax - xMin) + xMin
else:
       refaxisNodes[place+i,0] = numpy.mean(self.FFD.coef[sectionArr[i+skip,:,:],0])
       
 [same for y and z]

This way the default behavior will be exactly identical

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@marcomangano I think you are right and I was wrong about the way this piece of code was working

z_max = numpy.max(pts_vec[:, 2])

# Temporary ref axis nodes coordinates - aligned with main system of reference
x_nodes = xFraction * (x_max - x_min) + x_min # chordwise
marcomangano marked this conversation as resolved.
Show resolved Hide resolved
y_nodes = y_max - yFraction * (y_max - y_min) # top-bottom
z_nodes = z_max - zFraction * (z_max - z_min) # top-bottom

# These are the FFD ref axis nodes - if the block has not been rotated
nds = [x_nodes, y_nodes, z_nodes]
nds_final = numpy.copy(nds)

if rotType == 0 and rot0ang:
# rotating the non-aligned FFDs back in position
nds_final[:] = geo_utils.rotVbyW(nds, rot0axis, numpy.pi / 180 * (-rot0ang))

# insert the final coordinates in the var to be passed to pySpline:
refaxisNodes[place+i,0] = nds_final[0]
refaxisNodes[place+i,1] = nds_final[1]
refaxisNodes[place+i,2] = nds_final[2]

place += i + 1

# Add additional volumes
Expand All @@ -437,7 +494,7 @@ def addRefAxis(self, name, curve=None, xFraction=None, volumes=None,
curve = pySpline.Curve(X=refaxisNodes, k=2)
nAxis = len(curve.coef)
self.axis[name] = {'curve':curve, 'volumes':volumes,
'rotType':rotType, 'axis':axis,
'rotType':rotType, 'axis':axis, 'rot0ang':rot0ang, 'rot0axis':rot0axis,
'rotAxisVar':rotAxisVar}
else:
raise Error("One of 'curve' or 'xFraction' must be "
Expand Down Expand Up @@ -1188,6 +1245,16 @@ def updateCalculations(self, new_pts, isComplex, config):

for ipt in range(self.nPtAttach):
base_pt = self.refAxis.curves[self.curveIDs[ipt]](self.links_s[ipt])
# Variables for rotType = 0 rotation + scaling
ang = self.axis[self.curveIDNames[ipt]]['rot0ang']
ax_dir = self.axis[self.curveIDNames[ipt]]['rot0axis']
bp_ = numpy.copy(base_pt) # copy of original pointset - will not be rotated
if ang:
ang *= numpy.pi/180 # conv to [rad]
# Rotating the FFD according to inputs
# The FFD points should now be aligned with the main system of reference
base_pt = geo_utils.rotVbyW(bp_, ax_dir, ang)

scale = self.scale[self.curveIDNames[ipt]](self.links_s[ipt])
scale_x = self.scale_x[self.curveIDNames[ipt]](self.links_s[ipt])
scale_y = self.scale_y[self.curveIDNames[ipt]](self.links_s[ipt])
Expand All @@ -1199,12 +1266,32 @@ def updateCalculations(self, new_pts, isComplex, config):
self.curveIDs[ipt]].getDerivative(self.links_s[ipt])
deriv /= geo_utils.euclideanNorm(deriv) # Normalize
new_vec = -numpy.cross(deriv, self.links_n[ipt])
new_vec = geo_utils.rotVbyW(new_vec, deriv, self.rot_theta[
self.curveIDNames[ipt]](self.links_s[ipt])*numpy.pi/180)
if isComplex:
new_pts[ipt] = base_pt + new_vec*scale
new_pts[ipt] = bp_ + new_vec*scale # using "unrotated" bp_ vector
else:
new_pts[ipt] = numpy.real(base_pt + new_vec*scale)
new_pts[ipt] = numpy.real(bp_ + new_vec*scale)

if ang:
# Rotating to be aligned with main sys ref
nv_ = numpy.copy(new_vec)
new_vec = geo_utils.rotVbyW(nv_, ax_dir, ang)

# Apply scaling
new_vec[0] *= scale_x
new_vec[1] *= scale_y
new_vec[2] *= scale_z

if ang:
# Rotating back the scaled pointset to its original position
nv_rot = numpy.copy(new_vec) # nv_rot is scaled and rotated
new_vec = geo_utils.rotVbyW(nv_rot , ax_dir, -ang)

new_vec = geo_utils.rotVbyW(new_vec, deriv, self.rot_theta[self.curveIDNames[ipt]](self.links_s[ipt])*numpy.pi/180)

if isComplex:
new_pts[ipt] = bp_ + new_vec
else:
new_pts[ipt] = numpy.real(bp_ + new_vec)

else:
rotX = geo_utils.rotxM(self.rot_x[
Expand All @@ -1215,6 +1302,11 @@ def updateCalculations(self, new_pts, isComplex, config):
self.curveIDNames[ipt]](self.links_s[ipt]))

D = self.links_x[ipt]
if ang:
raise Warning("if rot0ang != 0, use rotType=0. The derivatives with other rotations are slightly off")
# rotate non-aligned FFDs
D_ = numpy.copy(D)
D = geo_utils.rotVbyW(D_, ax_dir, ang)

rotM = self._getRotMatrix(rotX, rotY, rotZ, rotType)

Expand Down Expand Up @@ -1246,6 +1338,13 @@ def updateCalculations(self, new_pts, isComplex, config):
D[0] *= scale_x
D[1] *= scale_y
D[2] *= scale_z
if ang:
# rotate non-aligned FFDs back to initial position
D_ = numpy.copy(D)
bp_rot = numpy.copy(base_pt) # here base_pt has been rotated
D = geo_utils.rotVbyW(D_ , ax_dir, -ang)
base_pt = geo_utils.rotVbyW(bp_rot, ax_dir, -ang)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am 99% sure it's fine to do this but only if the rotation is truly implemented as a linear operator as it should be

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good point. This came directly from the Mads/Sandy fix so I did not 100% thought it through, I will check the rotVbyW function.

Copy link
Contributor Author

@marcomangano marcomangano Mar 3, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the rotation function ... do you think we should enforce some normalization of the rotation axis?


if isComplex:
new_pts[ipt] = base_pt + D*scale
else:
Expand Down
35 changes: 32 additions & 3 deletions tests/reg_tests/commonUtils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,11 +176,17 @@ def totalSensitivityCS(DVGeo,nPt,ptName):

return dIdxCS

def testSensitivities(DVGeo,refDeriv,handler):
def testSensitivities(DVGeo,refDeriv,handler,pointset=1):
#create test points
points = numpy.zeros([2,3])
points[0,:] = [0.25,0,0]
points[1,:] = [-0.25,0,0]
if pointset == 1:
points[0,:] = [0.25,0,0]
points[1,:] = [-0.25,0,0]
elif pointset == 2:
points[0, :] = [0.25, 0.4, 4]
points[1, :] = [-0.8, 0.2, 7]
else:
raise Warning("Enter a valid pointset")

# add points to the geometry object
ptName = 'testPoints'
Expand Down Expand Up @@ -235,3 +241,26 @@ def testSensitivitiesD8(DVGeo,refDeriv,handler):
dIdx = DVGeo.totalSensitivity(dIdPt,ptName)

handler.root_add_dict("dIdx",dIdx,rtol=1e-7,atol=1e-7)

# --- Adding standard twist and single axis scaling functions ---
# These functions are added for Test 24 but could be extended to other tests

fix_root_sect=1
nRefAxPts = 4

def twist(val, geo):
axis_key = list(geo.axis.keys())[0]
for i in range(fix_root_sect, nRefAxPts):
geo.rot_theta[axis_key].coef[i] = val[i - fix_root_sect]

def thickness(val, geo):
axis_key = list(geo.axis.keys())[0]

for i in range(1, nRefAxPts):
geo.scale_z[axis_key].coef[i] = val[i - fix_root_sect]

def chord(val, geo):
axis_key = list(geo.axis.keys())[0]

for i in range(1, nRefAxPts):
geo.scale_x[axis_key].coef[i] = val[i - fix_root_sect]
31 changes: 31 additions & 0 deletions tests/reg_tests/ref/test_DVGeometry_23.ref
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
{
"RefAxis_nodes_coord": {
"__ndarray__": [
[
-0.4,
-0.09999999999999998,
0.0
],
[
-0.4,
-0.09999999999999998,
3.32194917
],
[
-0.4,
-0.09999999999999998,
5.78384945
],
[
-0.4,
-0.09999999999999998,
8.0
]
],
"dtype": "float64",
"shape": [
4,
3
]
}
}
Loading