diff --git a/pygeo/DVConstraints.py b/pygeo/DVConstraints.py index e35b866f..c483aeb8 100644 --- a/pygeo/DVConstraints.py +++ b/pygeo/DVConstraints.py @@ -4250,7 +4250,10 @@ def evalFunctions(self, funcs): cons.extend(self.jac[key].dot(self.DVGeo.DV_listLocal[key].value)) elif key in self.DVGeo.DV_listSectionLocal: cons.extend(self.jac[key].dot(self.DVGeo.DV_listSectionLocal[key].value)) - + elif key in self.DVGeo.DV_listSpanwiseLocal: + cons.extend(self.jac[key].dot(self.DVGeo.DV_listSpanwiseLocal[key].value)) + else: + raise Error(f"con {self.name} diffined wrt {key}, but {key} not found in DVGeo") funcs[self.name] = numpy.array(cons).real.astype('d') def evalFunctionsSens(self, funcsSens): @@ -4326,6 +4329,27 @@ def _finalize(self): self.ncon += len(cons) self.vizConIndices[key] = cons + # Section local shape variables + for key in self.DVGeo.DV_listSpanwiseLocal: + if self.config is None or self.config in self.DVGeo.DV_listSpanwiseLocal[key].config: + + # end for (indSet loop) + cons = self.DVGeo.DV_listSpanwiseLocal[key].mapIndexSets(self.indSetA,self.indSetB) + ncon = len(cons) + if ncon > 0: + # Now form the jacobian: + ndv = self.DVGeo.DV_listSpanwiseLocal[key].nVal + jacobian = numpy.zeros((ncon, ndv)) + for i in range(ncon): + jacobian[i, cons[i][0]] = self.factorA[i] + jacobian[i, cons[i][1]] = self.factorB[i] + self.jac[key] = jacobian + + # Add to the number of constraints and store indices which + # we need for tecplot visualization + self.ncon += len(cons) + self.vizConIndices[key] = cons + # with-respect-to are just the keys of the jacobian self.wrt = list(self.jac.keys()) diff --git a/pygeo/DVGeometry.py b/pygeo/DVGeometry.py index 5d78dae8..d5804d77 100644 --- a/pygeo/DVGeometry.py +++ b/pygeo/DVGeometry.py @@ -101,6 +101,7 @@ def __init__(self, fileName, complex=False, child=False, faceFreeze=None, name=N self.DV_listGlobal = OrderedDict() # Global Design Variable List self.DV_listLocal = OrderedDict() # Local Design Variable List self.DV_listSectionLocal = OrderedDict() # Local Normal Design Variable List + self.DV_listSpanwiseLocal = OrderedDict() # Local Normal Design Variable List # Coefficient rotation matrix dict for Section Local variables self.coefRotM = {} @@ -143,13 +144,15 @@ def __init__(self, fileName, complex=False, child=False, faceFreeze=None, name=N self.dCoefdXdvl = None # derivative counters for offsets - self.nDV_T = None + self.nDV_T = None # total number of design variables self.nDVG_T = None self.nDVL_T = None self.nDVSL_T = None - self.nDVG_count = 0 - self.nDVL_count = 0 - self.nDVSL_count = 0 + self.nDVSW_T = None + self.nDVG_count = 0 # number of global (G) variables + self.nDVL_count = 0 # number of local (L) variables + self.nDVSL_count = 0 # number of section (SL) local variables + self.nDVSW_count = 0 # number of spanwise (SW) local variables # The set of user supplied axis. self.axis = OrderedDict() @@ -658,7 +661,7 @@ def addGeoDVGlobal(self, dvName, value, func, lower=None, upper=None, # if the parent DVGeometry object has a name attribute, prepend it if self.name is not None: dvName = self.name + '_' + dvName - + if type(config) == str: config = [config] self.DV_listGlobal[dvName] = geoDVGlobal( @@ -762,6 +765,172 @@ def addGeoDVLocal(self, dvName, lower=None, upper=None, scale=1.0, return self.DV_listLocal[dvName].nVal + def addGeoDVSpanwiseLocal(self, dvName, spanIndex, axis='y', lower=None, upper=None, + scale=1.0, pointSelect=None, volList=None, config=None): + """ + Add one or more spanwise local design variables to the DVGeometry + object. Spanwise local variables are alternative form of local shape + variables used to apply equal DV changes in a chosen direction. + Some scenarios were this could be useful are: + + 1. 2D airfoil shape optimization. Because adflow works with 3D meshes, + 2D problems are represented my a mesh a single cell wide. Therefor, + to change the 2D representation of the airfoil both sides of the + mesh must be moved equally. This can be done with the addition of + linear constraints on a set of local shape variables, however this + approach requires more DVs than necessary (which complicates DV + sweeps) and the constaints are only enforced to a tolerance. Using + spanwise local design variables insures the airfoil is always + correctly represented in the 3D mesh using the correct amount of + design variables. + + 2. 3D wing optimization with constant airfoil shape. If the initial + wing geometry has a constant airfoil shape and constant chord, then + spanwise local dvs can be used to change the airfoil shape of the + wing while still keeping it constant along the span of the wing. + + Parameters + ---------- + dvName : str + A unique name to be given to this design variable group + + spanIndex : str, ('i', 'j', 'k') + the axis of the FFD along which the DVs are constant + all shape variables + + axis : str. Default is `y` + The coordinate directions to move. Permissible values are `x`, + `y` and `z`. If more than one direction is required, use multiple + calls to addGeoDVLocal with different axis values. + + lower : float + The lower bound for the variable(s). This will be applied to + all shape variables + + upper : float + The upper bound for the variable(s). This will be applied to + all shape variables + + scale : flot + The scaling of the variables. A good approximate scale to + start with is approximately 1.0/(upper-lower). This gives + variables that are of order ~1.0. + + pointSelect : pointSelect object. Default is None Use a + pointSelect object to select a subset of the total number + of control points. See the documentation for the + pointSelect class in geo_utils. Using pointSelect discards everything in + volList. + + volList : list + Use the control points on the volume indicies given in volList. + You should use pointSelect = None, otherwise this will not work. + + config : str or list + Define what configurations this design variable will be applied to + Use a string for a single configuration or a list for multiple + configurations. The default value of None implies that the design + variable applies to *ALL* configurations. + + Returns + ------- + N : int + The number of design variables added. + + Examples + -------- + >>> # Add all spanwise local variables + >>> # moving in the y direction, within +/- 0.5 units + >>> DVGeo.addGeoDVSpanwiseLocal("shape", 'k', lower=-0.5, upper=0.5, axis="z", scale=1.0) + """ + if type(config) == str: + config = [config] + + if pointSelect is not None: + if pointSelect.type != 'ijkBounds': + pts, ind = pointSelect.getPoints(self.FFD.coef) + else: + pts, ind = pointSelect.getPoints_ijk(self) + elif volList is not None: + if self.FFD.symmPlane is not None: + volListTmp = [] + for vol in volList: + volListTmp.append(vol) + for vol in volList: + volListTmp.append(vol+self.FFD.nVol/2) + volList = volListTmp + + volList = numpy.atleast_1d(volList).astype('int') + ind = [] + for iVol in volList: + ind.extend(self.FFD.topo.lIndex[iVol].flatten()) + ind = geo_utils.unique(ind) + else: + # Just take'em all + volList = numpy.arange(self.FFD.nVol) + ind = numpy.arange(len(self.FFD.coef)) + + + secLink = numpy.zeros(self.FFD.coef.shape[0], dtype=int) + secTransform = [numpy.eye(3)] + + if type(spanIndex) is str: + spanIndex = [spanIndex]*len(volList) + elif type(spanIndex) is list: + if len(spanIndex) != len(volList): + raise Error('If a list is given for spanIndex, the length must be' + ' equal to the length of volList.') + + ijk_2_idx = {'i':0, 'j':1, 'k':2} + + + volDVMap = [] + for ivol in volList: + + spanIdx = ijk_2_idx[spanIndex[ivol]] + lIndex = self.FFD.topo.lIndex[ivol] + + topo_shape = lIndex.shape + + # remove the span axis since all dv in that axis are linked + n_linked_coef = topo_shape[spanIdx] + dvs_shape = numpy.delete(topo_shape, spanIdx) + + # get total number of dvs + n_dvs = numpy.product(dvs_shape) + + # make a map from dvs to the ind that are controlled by that dv. + # (phrased another way) map from dv to all ind in the same span size position + dv_to_coef_ind = numpy.zeros((n_dvs,n_linked_coef), dtype='intc') + + + + # slice lIndex to get the indices of the coeffs that are in the same + # spanwise position + dv_idx = 0 + for i in range(dvs_shape[0]): + for j in range(dvs_shape[1]): + # no need to use fancy axis manipulation, since it doesn't need + # to be fast and if statements are expressive + if spanIndex[ivol] == 'i': + coef_ind = lIndex[:, i, j] + elif spanIndex[ivol] == 'j': + coef_ind = lIndex[i, :, j] + elif spanIndex[ivol] == 'k' : + coef_ind = lIndex[i, j, :] + + dv_to_coef_ind[dv_idx] = coef_ind + dv_idx += 1 + + # the for this volume is complete and can be added to the list of maps + volDVMap.append(dv_to_coef_ind) + + self.DV_listSpanwiseLocal[dvName] = geoDVSpanwiseLocal(dvName, lower, upper, + scale, axis, volDVMap, self.masks, + config) + + return self.DV_listSpanwiseLocal[dvName].nVal + def addGeoDVSectionLocal(self, dvName, secIndex, lower=None, upper=None, scale=1.0, axis=1, pointSelect=None, volList=None, orient0=None, orient2='svd', config=None): @@ -1118,6 +1287,16 @@ def setDesignVars(self, dvDict): len(vals_to_set))) self.DV_listSectionLocal[key].value = vals_to_set + if key in self.DV_listSpanwiseLocal: + vals_to_set = numpy.atleast_1d(dvDict[key]).astype('D') + if len(vals_to_set) != self.DV_listSpanwiseLocal[key].nVal: + raise Error('Incorrect number of design variables \ + for DV: %s.\nExpecting %d variables and received \ + %d variabes'%(key, self.DV_listSpanwiseLocal[key].nVal, + len(vals_to_set))) + self.DV_listSpanwiseLocal[key].value = vals_to_set + + # Jacobians are, in general, no longer up to date self.zeroJacobians(self.ptSetNames) @@ -1166,6 +1345,9 @@ def getValues(self): for key in self.DV_listSectionLocal: dvDict[key] = self.DV_listSectionLocal[key].value + # and now the Spanwise local DVs + for key in self.DV_listSpanwiseLocal: + dvDict[key] = self.DV_listSpanwiseLocal[key].value # Now call getValues on the children. This way the # returned dictionary will include the variables from @@ -1412,6 +1594,10 @@ def update(self, ptSetName, childDelta=True, config=None): numpy.put(self.FFD.coef[:, 1], self.ptAttachInd, temp[:, 1]) numpy.put(self.FFD.coef[:, 2], self.ptAttachInd, temp[:, 2]) + # Now add in the spanwise local DVs + for key in self.DV_listSpanwiseLocal: + self.DV_listSpanwiseLocal[key](self.FFD.coef, config) + # Now add in the section local DVs for key in self.DV_listSectionLocal: self.DV_listSectionLocal[key](self.FFD.coef, self.coefRotM, config) @@ -1438,6 +1624,8 @@ def update(self, ptSetName, childDelta=True, config=None): numpy.put(tempCoef[:, 2], self.ptAttachInd, new_pts[:, 2]) # Apply just the complex part of the local varibales + for key in self.DV_listSpanwiseLocal: + self.DV_listSpanwiseLocal[key].updateComplex(tempCoef, config) for key in self.DV_listSectionLocal: self.DV_listSectionLocal[key].updateComplex(tempCoef, self.coefRotM, config) for key in self.DV_listLocal: @@ -1548,7 +1736,7 @@ def convertSensitivityToDict(self, dIdx, out1D=False): """ # compute the various DV offsets - DVCountGlobal, DVCountLocal, DVCountSecLoc = self._getDVOffsets() + DVCountGlobal, DVCountLocal, DVCountSecLoc, DVCountSpanLoc = self._getDVOffsets() i = DVCountGlobal dIdxDict = {} @@ -1559,6 +1747,18 @@ def convertSensitivityToDict(self, dIdx, out1D=False): else: dIdxDict[dv.name] = dIdx[:, i:i+dv.nVal] i += dv.nVal + + + i = DVCountSpanLoc + for key in self.DV_listSpanwiseLocal: + dv = self.DV_listSpanwiseLocal[key] + if out1D: + dIdxDict[dv.name] = numpy.ravel(dIdx[:, i:i+dv.nVal]) + else: + dIdxDict[dv.name] = dIdx[:, i:i+dv.nVal] + i += dv.nVal + + i = DVCountSecLoc for key in self.DV_listSectionLocal: dv = self.DV_listSectionLocal[key] @@ -1567,6 +1767,7 @@ def convertSensitivityToDict(self, dIdx, out1D=False): else: dIdxDict[dv.name] = dIdx[:, i:i+dv.nVal] i += dv.nVal + i = DVCountLocal for key in self.DV_listLocal: dv = self.DV_listLocal[key] @@ -1607,24 +1808,32 @@ def convertDictToSensitivity(self, dIdxDict): dIdx : array Flattened array of length getNDV(). """ - DVCountGlobal, DVCountLocal, DVCountSecLoc = self._getDVOffsets() + DVCountGlobal, DVCountLocal, DVCountSecLoc, DVCountSpanLoc = self._getDVOffsets() dIdx = numpy.zeros(self.nDV_T, self.dtype) i = DVCountGlobal for key in self.DV_listGlobal: dv = self.DV_listGlobal[key] dIdx[i:i+dv.nVal] = dIdxDict[dv.name] i += dv.nVal + i = DVCountLocal for key in self.DV_listLocal: dv = self.DV_listLocal[key] dIdx[i:i+dv.nVal] = dIdxDict[dv.name] i += dv.nVal + i = DVCountSecLoc for key in self.DV_listSectionLocal: dv = self.DV_listSectionLocal[key] dIdx[i:i+dv.nVal] = dIdxDict[dv.name] i += dv.nVal + i = DVCountSpanLoc + for key in self.DV_listSpanLocal: + dv = self.DV_listSpanLocal[key] + dIdx[i:i+dv.nVal] = dIdxDict[dv.name] + i += dv.nVal + #Note: not sure if this works with (multiple) sibling child FFDs for iChild in range(len(self.children)): childdIdx = self.children[iChild].convertDictToSensitivity(dIdxDict) @@ -1644,6 +1853,7 @@ def getVarNames(self): names = list(self.DV_listGlobal.keys()) names.extend(list(self.DV_listLocal.keys())) names.extend(list(self.DV_listSectionLocal.keys())) + names.extend(list(self.DV_listSpanwiseLocal.keys())) # Call the children recursively for iChild in range(len(self.children)): @@ -1765,6 +1975,8 @@ def totalSensitivityProd(self, vec, ptSetName, comm=None, child=False, for key in names: if key in self.DV_listGlobal: dv = self.DV_listGlobal[key] + elif key in self.DV_listSpanwiseLocal: + dv = self.DV_listSpanwiseLocal[key] elif key in self.DV_listSectionLocal: dv = self.DV_listSectionLocal[key] else: @@ -1846,6 +2058,8 @@ def totalSensitivityTransProd(self, vec, ptSetName, comm=None, child=False, for key in names: if key in self.DV_listGlobal: dv = self.DV_listGlobal[key] + elif key in self.DV_listSpanwiseLocal: + dv = self.DV_listSpanwiseLocal[key] elif key in self.DV_listSectionLocal: dv = self.DV_listSectionLocal[key] else: @@ -1865,6 +2079,9 @@ def computeDVJacobian(self,config=None): # This is going to be DENSE in general J_attach = self._attachedPtJacobian(config=config) + # Compute local normal jacobian + J_spanwiselocal = self._spanwiselocalDVJacobian(config=config) + # Compute local normal jacobian J_sectionlocal = self._sectionlocalDVJacobian(config=config) @@ -1881,6 +2098,12 @@ def computeDVJacobian(self,config=None): if J_attach is not None: J_temp = sparse.lil_matrix(J_attach) + if J_spanwiselocal is not None: + if J_temp is None: + J_temp = sparse.lil_matrix(J_spanwiselocal) + else: + J_temp += J_spanwiselocal + if J_sectionlocal is not None: if J_temp is None: J_temp = sparse.lil_matrix(J_sectionlocal) @@ -1986,7 +2209,7 @@ def computeTotalJacobianCS(self, ptSetName, config=None): for child in self.children: child.nPts[ptSetName] = self.nPts[ptSetName] - DVGlobalCount, DVLocalCount, DVSecLocCount = self._getDVOffsets() + DVGlobalCount, DVLocalCount, DVSecLocCount, DVSpanLocCount = self._getDVOffsets() h = 1e-40j @@ -2012,6 +2235,24 @@ def computeTotalJacobianCS(self, ptSetName, config=None): self.DV_listGlobal[key].value[j] = refVal self._unComplexifyCoef() + for key in self.DV_listSpanwiseLocal: + for j in range(self.DV_listSpanwiseLocal[key].nVal): + if self.isChild: + self.FFD.coef = refFFDCoef.copy() + self.coef = refCoef.copy() + self.refAxis.coef = refCoef.copy() + self.refAxis._updateCurveCoef() + + refVal = self.DV_listSpanwiseLocal[key].value[j] + + self.DV_listSpanwiseLocal[key].value[j] += h + deriv = numpy.imag(self._update_deriv_cs(ptSetName,config=config).flatten())/numpy.imag(h) + + self.JT[ptSetName][DVSpanLocCount,:] = deriv + + DVSpanLocCount += 1 + self.DV_listSpanwiseLocal[key].value[j] = refVal + for key in self.DV_listSectionLocal: for j in range(self.DV_listSectionLocal[key].nVal): if self.isChild: @@ -2066,7 +2307,7 @@ def computeTotalJacobianCS(self, ptSetName, config=None): return def addVariablesPyOpt(self, optProb, globalVars=True, localVars=True, - sectionlocalVars=True, ignoreVars=None, freezeVars=None): + sectionlocalVars=True, spanwiselocalVars=True, ignoreVars=None, freezeVars=None): """ Add the current set of variables to the optProb object. @@ -2080,6 +2321,12 @@ def addVariablesPyOpt(self, optProb, globalVars=True, localVars=True, localVars : bool Flag specifying whether local variables are to be added + + sectionlocalVars : bool + Flag specifying whether section local variables are to be added + + spanwiselocalVars : bool + Flag specifying whether spanwiselocal variables are to be added ignoreVars : list of strings List of design variables the user DOESN'T want to use @@ -2099,10 +2346,13 @@ def addVariablesPyOpt(self, optProb, globalVars=True, localVars=True, # Add design variables from the master: varLists = OrderedDict([('globalVars',self.DV_listGlobal), ('localVars',self.DV_listLocal), - ('sectionlocalVars',self.DV_listSectionLocal)]) + ('sectionlocalVars',self.DV_listSectionLocal), + ('spanwiselocalVars',self.DV_listSpanwiseLocal)]) for lst in varLists: - if lst == 'globalVars' and globalVars or lst=='localVars' and localVars \ - or lst=='sectionlocalVars' and sectionlocalVars: + if lst == 'globalVars' and globalVars\ + or lst=='localVars' and localVars \ + or lst=='sectionlocalVars' and sectionlocalVars\ + or lst=='spanwiselocalVars' and spanwiselocalVars: for key in varLists[lst]: if key not in ignoreVars: dv = varLists[lst][key] @@ -2117,7 +2367,7 @@ def addVariablesPyOpt(self, optProb, globalVars=True, localVars=True, # Add variables from the children for child in self.children: - child.addVariablesPyOpt(optProb, globalVars, localVars, sectionlocalVars, + child.addVariablesPyOpt(optProb, globalVars, localVars, sectionlocalVars, spanwiselocalVars, ignoreVars, freezeVars) def writeTecplot(self, fileName): @@ -2369,7 +2619,7 @@ def demoDesignVars(self, directory, includeLocal=True, includeGlobal=True, if self.ptSetNames: pointSet = self.ptSetNames[0] else: - raise Error('DVGeo must have a point set to update for' + raise Error('DVGeo must have a point set to update for ' 'demoDesignVars to work.') else: writePointSet = True @@ -2385,11 +2635,19 @@ def demoDesignVars(self, directory, includeLocal=True, includeGlobal=True, continue lower = geo.DV_listLocal[key].lower upper = geo.DV_listLocal[key].upper + + elif key in geo.DV_listSpanwiseLocal: + if not includeLocal: + continue + lower = geo.DV_listSpanwiseLocal[key].lower + upper = geo.DV_listSpanwiseLocal[key].upper + elif key in geo.DV_listSectionLocal: if not includeLocal: continue lower = geo.DV_listSectionLocal[key].lower upper = geo.DV_listSectionLocal[key].upper + elif key in geo.DV_listGlobal: if not includeGlobal: continue @@ -2610,9 +2868,9 @@ def _getRotMatrix(self, rotX, rotY, rotZ, rotType): def _getNDV(self): """Return the actual number of design variables, global + local - + section local + + section local + spanwise local """ - return self._getNDVGlobal() + self._getNDVLocal() + self._getNDVSectionLocal() + return self._getNDVGlobal() + self._getNDVLocal() + self._getNDVSectionLocal() + self._getNDVSpanwiseLocal() def getNDV(self): """ @@ -2664,6 +2922,19 @@ def _getNDVSectionLocal(self): return nDV + def _getNDVSpanwiseLocal(self): + """ + Get total number of local variables, inclding any children + """ + nDV = 0 + for key in self.DV_listSpanwiseLocal: + nDV += self.DV_listSpanwiseLocal[key].nVal + + for child in self.children: + nDV += child._getNDVSpanwiseLocal() + + return nDV + def _getNDVSelf(self): """ Get total number of local and global variables, not including @@ -2704,6 +2975,17 @@ def _getNDVSectionLocalSelf(self): return nDV + def _getNDVSpanwiseLocalSelf(self): + """ + Get total number of local variables, not including + children + """ + nDV = 0 + for key in self.DV_listSpanwiseLocal: + nDV += self.DV_listSpanwiseLocal[key].nVal + + return nDV + def _getDVOffsets(self): ''' return the global and local DV offsets for this FFD @@ -2711,16 +2993,17 @@ def _getDVOffsets(self): # figure out the split between local and global Variables # All global vars at all levels come first - # then section local vars and then local vars. + # then spanwise, then section local vars and then local vars. # Parent Vars come before child Vars # get the global and local DV numbers on the parents if we don't have them if self.nDV_T==None or self.nDVG_T == None or self.nDVL_T==None \ - or self.nDVSL_T==None: + or self.nDVSL_T==None or self.nDVSW_T==None: self.nDV_T = self._getNDV() self.nDVG_T = self._getNDVGlobal() self.nDVL_T = self._getNDVLocal() self.nDVSL_T = self._getNDVSectionLocal() + self.nDVSW_T = self._getNDVSpanwiseLocal() self.nDVG_count = 0 self.nDVSL_count = self.nDVG_T self.nDVL_count = self.nDVG_T + self.nDVSL_T @@ -2728,6 +3011,7 @@ def _getDVOffsets(self): nDVG = self._getNDVGlobalSelf() nDVL = self._getNDVLocalSelf() nDVSL = self._getNDVSectionLocalSelf() + nDVSW = self._getNDVSpanwiseLocalSelf() # Set the total number of global and local DVs into any children of this parent for child in self.children: @@ -2737,16 +3021,19 @@ def _getDVOffsets(self): child.nDVG_T = self.nDVG_T child.nDVL_T = self.nDVL_T child.nDVSL_T = self.nDVSL_T + child.nDVSW_T = self.nDVSW_T child.nDVG_count = self.nDVG_count + nDVG child.nDVL_count = self.nDVL_count + nDVL child.nDVSL_count = self.nDVSL_count + nDVSL + child.nDVSW_count = self.nDVSW_count + nDVSL # Increment the counters for the children nDVG += child._getNDVGlobalSelf() nDVL += child._getNDVLocalSelf() nDVSL += child._getNDVSectionLocalSelf() + nDVSW += child._getNDVSpanwiseLocalSelf() - return self.nDVG_count, self.nDVL_count, self.nDVSL_count + return self.nDVG_count, self.nDVL_count, self.nDVSL_count, self.nDVSW_count def _update_deriv(self, iDV=0, h=1.0e-40j, oneoverh=1.0/1e-40, config=None, localDV=False): @@ -2860,6 +3147,10 @@ def _update_deriv_cs(self, ptSetName, config=None): numpy.put(self.FFD.coef[:, 2], self.ptAttachInd, new_pts[:, 2]) # Apply the real and complex parts separately + for key in self.DV_listSpanwiseLocal: + self.DV_listSpanwiseLocal[key](self.FFD.coef, self.coefRotM, config) + self.DV_listSpanwiseLocal[key].updateComplex(self.FFD.coef, self.coefRotM, config) + for key in self.DV_listSectionLocal: self.DV_listSectionLocal[key](self.FFD.coef, self.coefRotM, config) self.DV_listSectionLocal[key].updateComplex(self.FFD.coef, self.coefRotM, config) @@ -2988,7 +3279,7 @@ def computeTotalJacobianFD(self, ptSetName, config=None): for child in self.children: child.nPts[ptSetName] = self.nPts[ptSetName] - DVGlobalCount, DVLocalCount, DVSecLocCount = self._getDVOffsets() + DVGlobalCount, DVLocalCount, DVSecLocCount, DVSpanLocCount = self._getDVOffsets() h = 1e-6 @@ -3014,6 +3305,26 @@ def computeTotalJacobianFD(self, ptSetName, config=None): DVGlobalCount += 1 self.DV_listGlobal[key].value[j] = refVal + for key in self.DV_listSpanwiseLocal: + for j in range(self.DV_listSpanwiseLocal[key].nVal): + if self.isChild: + self.FFD.coef = refFFDCoef.copy() + self.coef = refCoef.copy() + self.refAxis.coef = refCoef.copy() + self.refAxis._updateCurveCoef() + + refVal = self.DV_listSpanwiseLocal[key].value[j] + + self.DV_listSpanwiseLocal[key].value[j] += h + coordsph = self.update(ptSetName, childDelta=False, config=config).flatten() + + deriv = (coordsph-coords0)/h + self.JT[ptSetName][DVSpanLocCount,:]=deriv + + DVSpanLocCount += 1 + self.DV_listSpanwiseLocal[key].value[j] = refVal + + for key in self.DV_listSectionLocal: for j in range(self.DV_listSectionLocal[key].nVal): if self.isChild: @@ -3142,6 +3453,102 @@ def _attachedPtJacobian(self, config): return Jacobian + + def _spanwiselocalDVJacobian(self, config=None): + """ + Return the derivative of the coefficients wrt the local normal design + variables + """ + # This is relatively straight forward, since the matrix is + # entirely one's or zeros + nDV = self._getNDVSpanwiseLocalSelf() + self._getDVOffsets() + + if nDV != 0: + Jacobian = sparse.lil_matrix((self.nPtAttachFull*3, self.nDV_T)) + + # Create the storage arrays for the information that must be + # passed to the children + + for iChild in range(len(self.children)): + N = self.FFD.embededVolumes['child%d_axis'%(iChild)].N + self.children[iChild].dXrefdXdvl = numpy.zeros((N*3, self.nDV_T)) + + N = self.FFD.embededVolumes['child%d_coef'%(iChild)].N + self.children[iChild].dCcdXdvl = numpy.zeros((N*3, self.nDV_T)) + + iDVSpanwiseLocal = self.nDVSW_count + for key in self.DV_listSpanwiseLocal: + dv = self.DV_listSpanwiseLocal[key] + + # check that the dv is active for this config + if dv.config is None or config is None or any(c0 == config for c0 in dv.config): + nVal = dv.nVal + + #apply this dv to FFD + self.DV_listSpanwiseLocal[key](self.FFD.coef, config) + + # loop over value of the dv + # (for example a single shape dv may have 20 values that + # control the shape of the FFD at 20 points) + for j in range(nVal): + coefs = dv.dv_to_coefs[j] # affected control points of FFD + + # this is map from dvs to coef + for coef in coefs: + irow = coef*3 + dv.axis + # *3 because the jacobian has a row for each x,y,z of the FFD + # It is basically + # row number = coef index * n dimensions + dimension index + + # value of FFD node location = x0 + dv_SWLocal[j] + # so partial(FFD node location)/partial(dv_SWLocal) = 1 + # for each node effected by the dv_SWLocal[j] + Jacobian[irow, iDVSpanwiseLocal] = 1.0 + + for iChild in range(len(self.children)): + # Get derivatives of child ref axis and FFD control + # points w.r.t. parent's FFD control points + dXrefdCoef = self.FFD.embededVolumes['child%d_axis'%(iChild)].dPtdCoef + dCcdCoef = self.FFD.embededVolumes['child%d_coef'%(iChild)].dPtdCoef + + # derivative of Change in the FFD coef due to DVs + # same as Jacobian above, but differnt ordering + dCoefdXdvl = numpy.zeros(self.FFD.coef.shape,dtype='d') + + for coef in coefs: + dCoefdXdvl[coef, dv.axis] = 1.0 + + + dXrefdXdvl = numpy.zeros((dXrefdCoef.shape[0]*3),'d') + dCcdXdvl = numpy.zeros((dCcdCoef.shape[0]*3),'d') + + dXrefdXdvl[0::3] = dXrefdCoef.dot(dCoefdXdvl[:, 0]) + dXrefdXdvl[1::3] = dXrefdCoef.dot(dCoefdXdvl[:, 1]) + dXrefdXdvl[2::3] = dXrefdCoef.dot(dCoefdXdvl[:, 2]) + + dCcdXdvl[0::3] = dCcdCoef.dot(dCoefdXdvl[:, 0]) + dCcdXdvl[1::3] = dCcdCoef.dot(dCoefdXdvl[:, 1]) + dCcdXdvl[2::3] = dCcdCoef.dot(dCoefdXdvl[:, 2]) + + # TODO: the += here is to allow recursion check this with multiple nesting + # levels + self.children[iChild].dXrefdXdvl[:, iDVLocal] += dXrefdXdvl + self.children[iChild].dCcdXdvl[:, iDVLocal] += dCcdXdvl + + iDVSpanwiseLocal += 1 + else: + iDVSpanwiseLocal += self.DV_listSectionLocal[key].nVal + + # end if config check + # end for + else: + Jacobian = None + + + return Jacobian + + def _sectionlocalDVJacobian(self, config=None): """ Return the derivative of the coefficients wrt the local normal design @@ -3483,7 +3890,7 @@ def checkDerivatives(self, ptSetName): h = 1e-6 # figure out the split between local and global Variables - DVCountGlob, DVCountLoc, DVCountSecLoc = self._getDVOffsets() + DVCountGlob, DVCountLoc, DVCountSecLoc, DVCountSpanLoc = self._getDVOffsets() for key in self.DV_listGlobal: for j in range(self.DV_listGlobal[key].nVal): @@ -3580,6 +3987,38 @@ def checkDerivatives(self, ptSetName): DVCountSecLoc += 1 self.DV_listSectionLocal[key].value[j] = refVal + for key in self.DV_listSpanwiseLocal: + for j in range(self.DV_listSpanwiseLocal[key].nVal): + + print('========================================') + print(' SpanwiseLocalVar(%s), Value(%d) '%(key, j)) + print('========================================') + + if self.isChild: + self.FFD.coef = refFFDCoef.copy() + self.coef = refCoef.copy() + self.refAxis.coef = self.coef.copy() + self.refAxis._updateCurveCoef() + + refVal = self.DV_listSpanwiseLocal[key].value[j] + + self.DV_listSpanwiseLocal[key].value[j] += h + coordsph = self.update(ptSetName).flatten() + + deriv = (coordsph-coords0)/h + + for ii in range(len(deriv)): + relErr = (deriv[ii] - Jac[DVCountSpanLoc, ii])/( + 1e-16 + Jac[DVCountSpanLoc, ii]) + absErr = deriv[ii] - Jac[DVCountSpanLoc,ii] + + if abs(relErr) > h and abs(absErr) > h: + print(ii, deriv[ii], Jac[DVCountSpanLoc, ii], relErr, absErr) + # print(ii, deriv[ii], Jac[DVCountSpanLoc, ii], relErr, absErr) + + DVCountSpanLoc += 1 + self.DV_listSpanwiseLocal[key].value[j] = refVal + for child in self.children: child.checkDerivatives(ptSetName) @@ -3725,7 +4164,7 @@ def sectionFrame(self, sectionIndex, sectionTransform, sectionLink, ivol=0, ax2 /= numpy.linalg.norm(ax2) else: raise Error('orient2 must be \'svd\' or \'ffd\'') - + # Options for choosing in-plane axes # 1. Align axis '0' with projection of the given vector on section # plane. @@ -3899,6 +4338,104 @@ def _convertTo1D(value, dim1): raise Error('The size of the 1D array was the incorret shape') +class geoDVSpanwiseLocal(geoDVLocal): + + def __init__(self, dvName, lower, upper, scale, axis, vol_dv_to_coefs, mask, config): + + """Create a set of geometric design variables which change the shape + of a surface surface_id. Local design variables change the surface + in all three axis. + See addGeoDVLocal for more information + """ + + self.dv_to_coefs = [] + + + # add all the coefs to a flat array, but check that it isn't masked first + for ivol in range(len(vol_dv_to_coefs)): + for loc_dv in range(len(vol_dv_to_coefs[ivol])): + coefs = vol_dv_to_coefs[ivol][loc_dv] + + loc_dv_to_coefs = [] + + #loop through each of coefs to see if it is masked + for coef in coefs: + if mask[coef]==False: + loc_dv_to_coefs.append(coef) + + + self.dv_to_coefs.append(loc_dv_to_coefs) + + + if 'x' == axis.lower(): + self.axis = 0 + elif 'y' == axis.lower(): + self.axis = 1 + elif 'z' == axis.lower(): + self.axis = 2 + else: + raise NotImplementedError + + self.nVal = len(self.dv_to_coefs) + self.value = numpy.zeros(self.nVal, 'D') + + self.name = dvName + self.lower = None + self.upper = None + self.config = config + + if lower is not None: + self.lower = _convertTo1D(lower, self.nVal) + if upper is not None: + self.upper = _convertTo1D(upper, self.nVal) + if scale is not None: + self.scale = _convertTo1D(scale, self.nVal) + + + def __call__(self, coef, config): + """When the object is called, apply the design variable values to + coefficients""" + if self.config is None or config is None or any(c0 == config for c0 in self.config): + for i in range(self.nVal): + coef[self.dv_to_coefs[i], self.axis] += self.value[i].real + + return coef + + def updateComplex(self, coef, config): + if self.config is None or config is None or any(c0 == config for c0 in self.config): + for i in range(self.nVal): + coef[self.dv_to_coefs[i], self.axis] += self.value[i].imag*1j + + return coef + + def mapIndexSets(self,indSetA,indSetB): + ''' + Map the index sets from the full coefficient indices to the local set. + ''' + + cons = [] + for j in range(len(indSetA)): + # Try to find this index # in the coefList (temp) + up = None + down = None + + # Note: We are doing inefficient double looping here + for idx_dv, coefs in enumerate(self.dv_to_coefs): + + for coef in coefs: + + if coef == indSetA[j]: + up = idx_dv + if coef == indSetB[j]: + down = idx_dv + + # If we haven't found up AND down do nothing + if up is not None and down is not None: + cons.append([up, down]) + + return cons + + class geoDVSectionLocal(object): def __init__(self, dvName, lower, upper, scale, axis, coefListIn, mask, diff --git a/tests/reg_tests/ref/test_Cylinder_spanwise_dvs.ref b/tests/reg_tests/ref/test_Cylinder_spanwise_dvs.ref new file mode 100644 index 00000000..99019bcf --- /dev/null +++ b/tests/reg_tests/ref/test_Cylinder_spanwise_dvs.ref @@ -0,0 +1,75 @@ +{ + "funcs1": { + "DVCon1_leradius_constraints_0": { + "__ndarray__": [ + 0.9999962914850281, + 0.999995760958035, + 0.9999998221258696, + 0.9999993319711181, + 1.0000049447559047 + ], + "dtype": "float64", + "shape": [ + 5 + ] + } + }, + "funcs2": { + "DVCon1_leradius_constraints_0": { + "__ndarray__": [ + 1.0274767606779724, + 1.0389909669351234, + 1.0523915707811478, + 1.0676870537795249, + 1.0849072587266027 + ], + "dtype": "float64", + "shape": [ + 5 + ] + } + }, + "funcsSens": { + "DVCon1_leradius_constraints_0": { + "shape": { + "__ndarray__": [ + [ + -0.25686901113785543, + 0.2564858233025938, + -0.2564858233025938, + 0.25686901113785543 + ], + [ + -0.3281090443426573, + 0.32610410944473645, + -0.25341375902391494, + 0.25541869392183564 + ], + [ + -0.4081754150387783, + 0.40504333499679473, + -0.24239958498968017, + 0.24553166503166368 + ], + [ + -0.49722701453849344, + 0.4934592622297184, + -0.22327209583838797, + 0.2270398481471629 + ], + [ + -0.5954377625541436, + 0.5915177149773407, + -0.19586588913351255, + 0.1997859367103156 + ] + ], + "dtype": "float64", + "shape": [ + 5, + 4 + ] + } + } + } +} \ No newline at end of file diff --git a/tests/reg_tests/test_Cylinder.py b/tests/reg_tests/test_Cylinder.py index 9a2c23ba..0197d5ea 100644 --- a/tests/reg_tests/test_Cylinder.py +++ b/tests/reg_tests/test_Cylinder.py @@ -142,3 +142,60 @@ def scale_circle(val, geo): print(funcsSens) + def train_spanwise_dvs(self, train=True, refDeriv=True): + self.test_spanwise_dvs(train=train, refDeriv=refDeriv) + + + def test_spanwise_dvs(self, train=False, refDeriv=False): + refFile = os.path.join(self.base_path,'ref/test_Cylinder_spanwise_dvs.ref') + + with BaseRegTest(refFile, train=train) as handler: + handler.root_print("Test 1: Basic FFD, global DVs") + radius = 1.0 + height = 10.0 + + DVCon = DVConstraints() + surf = self.make_cylinder_mesh(radius, height) + DVCon.setSurface(surf) + # DVCon.writeSurfaceTecplot('cylinder_surface.dat') + + ffd_name = os.path.join(self.base_path,'../inputFiles/cylinder_ffd.xyz') + self.make_ffd(ffd_name, radius, height) + DVGeo = DVGeometry(ffd_name) + + DVGeo.addGeoDVSpanwiseLocal("shape", 'i', lower=-0.5, upper=0.5, axis="y", scale=1.0) + + size = DVGeo._getNDVSpanwiseLocal() + DVCon.setDVGeo(DVGeo) + + + leList = [[0, 0, 0 ], [-radius/2, 0, height]] + xAxis = [-1, 0, 0] + yAxis = [0, 1, 0] + DVCon.addLERadiusConstraints(leList, nSpan=5, axis=yAxis, + chordDir=xAxis, scaled=False) + # DVCon.writeTecplot('cylinder_constraints.dat') + + funcs = {} + DVCon.evalFunctions(funcs) + print(funcs) + handler.root_add_dict('funcs1', funcs, rtol=1e-6, atol=1e-6) + + numpy.random.seed(0) + DVGeo.setDesignVars({'shape':(numpy.random.rand(size) - 0.5)}) + + funcs = {} + DVCon.evalFunctions(funcs) + handler.root_add_dict('funcs2', funcs, rtol=1e-6, atol=1e-6) + print(funcs) + + funcsSens = {} + DVCon.evalFunctionsSens(funcsSens) + print(funcsSens) + handler.root_add_dict('funcsSens', funcsSens, rtol=1e-6, atol=1e-6) + print(funcsSens) + + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/reg_tests/test_DVGeometry.py b/tests/reg_tests/test_DVGeometry.py index 0a34a675..41f30656 100644 --- a/tests/reg_tests/test_DVGeometry.py +++ b/tests/reg_tests/test_DVGeometry.py @@ -4,8 +4,9 @@ import numpy from baseclasses import BaseRegTest import commonUtils -from pygeo import geo_utils, DVGeometry +from pygeo import geo_utils, DVGeometry, DVConstraints from parameterized import parameterized +from stl import mesh class RegTestPyGeo(unittest.TestCase): @@ -805,6 +806,59 @@ def test_22(self, train=False, refDeriv=False): handler.par_add_norm('norm', norm_diff, rtol=1e-7, atol=1e-7) os.remove(copyName) + + def test_spanwise_dvs(self, train=False, refDeriv=False): + """ + Test spanwise_dvs + """ + # refFile = os.path.join(self.base_path,'ref/test_DVGeometry_spanwise_dvs.ref') + # with BaseRegTest(refFile, train=train) as handler: + # handler.root_print("Test spanwise local variables writing function") + + meshfile = os.path.join(self.base_path, '../inputFiles/c172.stl') + ffdfile = os.path.join(self.base_path, '../inputFiles/c172.xyz') + testmesh = mesh.Mesh.from_file(meshfile) + # test mesh dim 0 is triangle index + # dim 1 is each vertex of the triangle + # dim 2 is x, y, z dimension + + # create a DVGeo object with a few local thickness variables + DVGeo = DVGeometry(ffdfile) + DVGeo.addGeoDVSpanwiseLocal("shape", 'i', lower=-0.5, upper=0.5, axis="y", scale=1.0) + + + # create a DVConstraints object for the wing + DVCon =DVConstraints() + DVCon.setDVGeo(DVGeo) + p0 = testmesh.vectors[:,0,:] / 1000 + v1 = testmesh.vectors[:,1,:] / 1000 - p0 + v2 = testmesh.vectors[:,2,:] / 1000 - p0 + DVCon.setSurface([p0, v1, v2]) + + + leList = [[0.7, 0.0, 0.1],[0.7, 0.0, 2.4]] + teList = [[0.9, 0.0, 0.1],[0.9, 0.0, 2.4]] + + nSpan = 10 + nChord = 10 + name = "thickness_con" + DVCon.addThicknessConstraints2D(leList, teList, nSpan, nChord, name=name) + + size = DVGeo._getNDVSpanwiseLocal() + + numpy.random.seed(0) + DVGeo.setDesignVars({'shape':(numpy.random.rand(size) - 0.5)}) + + funcs = {} + DVCon.evalFunctions(funcs) + # print(funcs) + + for i in range(nChord): + for j in range(nSpan-1): + numpy.testing.assert_allclose(funcs[name][i*nChord + j+1], funcs[name][i*nChord + j], rtol=2e-15) + + + def train_23_xyzFraction(self, train=True): self.test_23_xyzFraction(train=train)