Added coordinate transformations and a few fixes

pygeo/parameterization/DVGeo.py

@@ -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
@@ -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,
@@ -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
 
@@ -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
 
         """
 
@@ -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
@@ -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")
@@ -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")
@@ -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")
@@ -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.
 
@@ -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
         -------
@@ -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")
@@ -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
 
@@ -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):
@@ -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)
 
@@ -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)
 
@@ -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

pygeo/parameterization/designVars.py

@@ -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
 

pygeo/pyBlock.py

@@ -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.