DVGeometryVSP.py

# ======================================================================
#         Imports
# ======================================================================
from __future__ import print_function
import copy
import tempfile
import shutil
import os
import sys
import time
import numpy
from collections import OrderedDict
from scipy import sparse
from scipy.spatial import cKDTree
from mpi4py import MPI
from pyspline import pySpline
import vsp
from cgtlib import lsect

class Error(Exception):
    """
    Format the error message in a box to make it clear this
    was a explicitly raised exception.
    """
    def __init__(self, message):
        msg = '\n+'+'-'*78+'+'+'\n' + '| DVGeometryVSP Error: '
        i = 22
        for word in message.split():
            if len(word) + i + 1 > 78: # Finish line and start new one
                msg += ' '*(78-i)+'|\n| ' + word + ' '
                i = 1 + len(word)+1
            else:
                msg += word + ' '
                i += len(word)+1
        msg += ' '*(78-i) + '|\n' + '+'+'-'*78+'+'+'\n'
        print(msg)
        Exception.__init__(self)

class DVGeometryVSP(object):
    """A class for manipulating VSP geometry

    The purpose of the DVGeometryVSP class is to provide translation
    of the VSP geometry engine to externally supplied surfaces. This
    allows the use of VSP design variables to control the MACH
    framework.

    There are several import limitations:

    1. Since VSP is surface based only, it cannot be used to
    parameterize a geometry that doesn't lie on the surface. This
    means it cannot be used for structural analysis. It is generally
    possible use most of the constraints DVConstraints since most of
    those points lie on the surface.

    2. It cannot handle *moving* intersection. A geometry with static
    intersections is fine as long as the intersection doesn't move

    3. It does not support complex numbers for the complex-step
    method.

    4. It does not surpport separate configurations.

    Parameters
    ----------
    vspFile : str
       filename of .vsp3 file.

    comm : MPI Intra Comm
       Comm on which to build operate the object. This is used to
       perform embarasisngly parallel finite differencing. Defaults to
       MPI.COMM_WORLD.

    scale : float
       A global scale factor from the VSP geometry to incoming (CFD) mesh
       geometry. For example, if the VSP model is in inches, and the CFD
       in meters, scale=0.0254.

    comps : list of strings
       A list of string defining the subset of the VSP components to use when
       exporting the P3D surface files

    Examples
    --------
    The general sequence of operations for using DVGeometry is as follows:
      >>> from pygeo import *
      >>> DVGeo = DVGeometryVSP("wing.vsp3", MPI_COMM_WORLD)
      >>> # Add a set of coordinates Xpt into the object
      >>> DVGeo.addPointSet(Xpt, 'myPoints')
      >>>

    """
    def __init__(self, vspFile, comm=MPI.COMM_WORLD, scale=1.0, comps=[],
                 intersectedComps=None, debug=False):
        self.points = OrderedDict()
        self.pointSets = OrderedDict()
        self.updated = {}
        self.vspScale = scale
        self.comm = comm
        self.vspFile = vspFile
        self.debug = debug
        self.jac = None
        # Load in the VSP model
        vsp.ClearVSPModel()
        vsp.ReadVSPFile(vspFile)

        # Setup the export group set (0) with just the sets we want.
        self.exportSet = 9

        # List of all componets returned from VSP. Note that this
        # order is important...it is the order the comps will be
        # written out in plot3d format.
        allComps = vsp.FindGeoms()

        # If we were not given comps, use all of them
        if comps == []:
            for c in allComps:
                comps.append(vsp.GetContainerName(c))

        # First set the export flag for exportSet to False for everyone
        for comp in allComps:
            vsp.SetSetFlag(comp, self.exportSet, False)

        self.exportComps = []
        for comp in allComps:
            # Check if this one is in our list:
            compName = vsp.GetContainerName(comp)
            if compName in comps:
                vsp.SetSetFlag(comp, self.exportSet, True)
                self.exportComps.append(compName)

        # Create a directory in which we will put the temporary files
        # we need. We *should* use something like tmepfile.mkdtemp()
        # but that behaves badly on pleiades.
        tmpDir = None
        if self.comm.rank == 0:
            tmpDir = './tmpDir_%d_%s'%(MPI.COMM_WORLD.rank, time.time())
            print ('Temp dir is: %s'%tmpDir)
            if not os.path.exists(tmpDir):
                os.makedirs(tmpDir)
        self.tmpDir = self.comm.bcast(tmpDir)

        # Initial list of DVs
        self.DVs = OrderedDict()

        # Run the update. This will also set the conn on the first pass
        self.conn = None
        self.pts0 = None
        self.cumSizes = None
        self.sizes = None
        if self.comm.rank == 0:
            self.pts0, self.conn, self.cumSizes, self.sizes = self._getUpdatedSurface()

        self.pts0 = self.comm.bcast(self.pts0)
        self.conn = self.comm.bcast(self.conn)
        self.cumSizes = self.comm.bcast(self.cumSizes)
        self.sizes = self.comm.bcast(self.sizes)
        self.pts = self.pts0.copy()

        # Finally process theintersection information. We had to wait
        # until all processors have the surface information.
        self.intersectComps = []
        if intersectedComps is None:
            intersectedComps = []
        for i in range(len(intersectedComps)):
            c = intersectedComps[i]
            # Get the index of each of the two comps:
            compIndexA = self.exportComps.index(c[0])
            compIndexB = self.exportComps.index(c[1])
            direction = c[2]
            dStar = c[3]
            extraComps = []
            if len(c) == 5:
                for j in range(len(c[4])):
                    extraComps.append(self.exportComps.index(c[4][j]))

            self.intersectComps.append(CompIntersection(
                compIndexA, compIndexB, extraComps, direction, dStar,
                self.pts, self.cumSizes, self.sizes, self.tmpDir))

    def addPointSet(self, points, ptName, **kwargs):
        """
        Add a set of coordinates to DVGeometry

        The is the main way that geometry, in the form of a coordinate
        list is given to DVGeoemtry to be manipulated.

        Parameters
        ----------
        points : array, size (N,3)
            The coordinates to embed. These cordinates *should* all
            project into the interior of the FFD volume.
        ptName : str
            A user supplied name to associate with the set of
            coordinates. Thisname will need to be provided when
            updating the coordinates or when getting the derivatives
            of the coordinates.
        """

        # save this name so that we can zero out the jacobians properly
        self.points[ptName] = True # ADFlow checks self.points to see
                                   # if something is added or not.

        points = numpy.array(points).real.astype('d')
        # Attach the points to the original surface using the fast ADT
        # projection code from pySpline. Note that we convert to
        # 1-based indexing for conn here.

        if len(points) > 0:
            faceID, uv = pySpline.libspline.adtprojections(
                self.pts0.T, (self.conn+1).T, points.T)
            uv = uv.T
            faceID -= 1 # Convert back to zero-based indexing.
        else:
            faceID = numpy.zeros((0), 'intc')
            uv = numpy.zeros((0, 2), 'intc')

        # From the faceID we can back out what component each one is
        # connected to. This way if we have intersecting components we
        # only change the ones that are apart of the two surfaces.
        cumFaceSizes = numpy.zeros(len(self.sizes) + 1, 'intc')
        for i in range(len(self.sizes)):
            nCellI = self.sizes[i][0]-1
            nCellJ = self.sizes[i][1]-1
            cumFaceSizes[i+1] = cumFaceSizes[i] + nCellI*nCellJ
        compIDs = numpy.searchsorted(cumFaceSizes, faceID, side='right')-1

        # Compute the offsets by evaluating the new points.
        # eval which should be fast enough
        newPts = numpy.zeros_like(points)
        for idim in range(3):
            newPts[:, idim] = \
                              (1-uv[:, 0])*(1 - uv[:, 1]) * self.pts0[self.conn[faceID, 0], idim] + \
                              (  uv[:, 0])*(1 - uv[:, 1]) * self.pts0[self.conn[faceID, 1], idim] + \
                              (  uv[:, 0])*(    uv[:, 1]) * self.pts0[self.conn[faceID, 2], idim] + \
                              (1-uv[:, 0])*(    uv[:, 1]) * self.pts0[self.conn[faceID, 3], idim]
        offset = newPts - points

        # Create the little class with the data
        self.pointSets[ptName] = PointSet(points, faceID, uv, offset,
                                          self.conn, len(self.pts0))

        # Add the points to each of the intersection curve objects if we have any
        for IC in self.intersectComps:
            IC.addPointSet(points, compIDs, ptName)

        self.updated[ptName] = False

    def setDesignVars(self, dvDict, updateJacobian=True):
        """
        Standard routine for setting design variables from a design
        variable dictionary.

        Parameters
        ----------
        dvDict : dict
            Dictionary of design variables. The keys of the dictionary
            must correspond to the design variable names. Any
            additional keys in the dfvdictionary are simply ignored.
            """

        # Just dump in the values
        for key in dvDict:
            if key in self.DVs:
                self.DVs[key].value = dvDict[key]

        # When we set the design variables we will also compute the
        # new surface self.pts. Only one proc needs to do this. This
        # should be the only location where self.pts is set.
        if self.comm.rank == 0:
            self.pts, conn, cumSizes, sizes = self._getUpdatedSurface()
        self.pts = self.comm.bcast(self.pts)

        # We need to give the updated coordinates to each of the
        # intersectComps (if we have any) so they can update the new
        # intersection curve
        for IC in self.intersectComps:
            IC.setSurface(self.pts)

        # We will also compute the jacobian so it is also up to date,
        # provided we are asked for ti
        if updateJacobian:
            self._computeSurfJacobian()

        # Flag all the pointSets as not being up to date:
        for pointSet in self.updated:
            self.updated[pointSet] = False

    def getValues(self):
        """
        Generic routine to return the current set of design
        variables. Values are returned in a dictionary format
        that would be suitable for a subsequent call to setValues()

        Returns
        -------
        dvDict : dict
            Dictionary of design variables
        """
        dvDict = OrderedDict()
        for dvName in self.DVs:
            dvDict[dvName] = self.DVs[dvName].value

        return dvDict

    def update(self, ptSetName, config=None):
        """This is the main routine for returning coordinates that have been
        updated by design variables. Multiple configs are not
        supported.

        Parameters
        ----------
        ptSetName : str
            Name of point-set to return. This must match ones of the
            given in an :func:`addPointSet()` call.
        """

        # Since we are sure that self.pts is always up to date, we can
        # just do the eval + offset.
        n      = len(self.pointSets[ptSetName].points)
        uv     = self.pointSets[ptSetName].uv
        faceID = self.pointSets[ptSetName].faceID
        offset = self.pointSets[ptSetName].offset
        newPts = numpy.zeros((n, 3))

        for idim in range(3):
            newPts[:, idim] = \
                              (1-uv[:, 0])*(1 - uv[:, 1]) * self.pts[self.conn[faceID, 0], idim] + \
                              (  uv[:, 0])*(1 - uv[:, 1]) * self.pts[self.conn[faceID, 1], idim] + \
                              (  uv[:, 0])*(    uv[:, 1]) * self.pts[self.conn[faceID, 2], idim] + \
                              (1-uv[:, 0])*(    uv[:, 1]) * self.pts[self.conn[faceID, 3], idim]
        newPts -= offset

        # Now compute the delta between the nominal new poitns and the
        # original points:
        delta = newPts - self.pointSets[ptSetName].points

        # Potentially modify the delta according to our intersection curves
        for IC in self.intersectComps:
            delta = IC.update(ptSetName, delta)

        # Now get the final newPts with a possibly modified delta.
        newPts = self.pointSets[ptSetName].points + delta

        # Finally flag this pointSet as being up to date:
        self.updated[ptSetName] = True

        return newPts

    def writeVSPFile(self, fileName, exportSet=0):
        """Take the current design and write a new FSP file"""
        vsp.WriteVSPFile(fileName, exportSet)

    def pointSetUpToDate(self, ptSetName):
        """
        This is used externally to query if the object needs to update
        its pointset or not. Essentially what happens, is when
        update() is called with a point set, it the self.updated dict
        entry for pointSet is flagged as true. Here we just return
        that flag. When design variables are set, we then reset all
        the flags to False since, when DVs are set, nothing (in
        general) will up to date anymore.

        Parameters
        ----------
        ptSetName : str
            The name of the pointset to check.
        """
        if ptSetName in self.updated:
            return self.updated[ptSetName]
        else:
            return True

    def getNDV(self):
        """ Return the number of DVs"""
        return len(self.DVs)

    def getVarNames(self):
        """
        Return a list of the design variable names. This is typically
        used when specifying a wrt= argument for pyOptSparse.

        Examples
        --------
        optProb.addCon(.....wrt=DVGeo.getVarNames())
        """
        return list(self.DVs.keys())


    def totalSensitivity(self, dIdpt, ptSetName, comm=None, config=None):
        """
        This function computes sensitivty information.

        Specificly, it computes the following:
        :math:`\\frac{dX_{pt}}{dX_{DV}}^T \\frac{dI}{d_{pt}}

        Parameters
        ----------
        dIdpt : array of size (Npt, 3) or (N, Npt, 3)

            This is the total derivative of the objective or function
            of interest with respect to the coordinates in
            'ptSetName'. This can be a single array of size (Npt, 3)
            **or** a group of N vectors of size (Npt, 3, N). If you
            have many to do, it is faster to do many at once.

        ptSetName : str
            The name of set of points we are dealing with

        comm : MPI.IntraComm
            The communicator to use to reduce the final derivative. If
            comm is None, no reduction takes place.

        Returns
        -------
        dIdxDict : dic
            The dictionary containing the derivatives, suitable for
            pyOptSparse
        """

        # We may not have set the variables so the surf jac might not be computed.
        if self.jac is None:
            self._computeSurfJacobian()

        # Make dIdpt at least 3D
        if len(dIdpt.shape) == 2:
            dIdpt = numpy.array([dIdpt])
        N = dIdpt.shape[0]

        nDV = self.getNDV()
        dIdx_local = numpy.zeros((N, nDV), 'd')

        # The following code computes the final sensitivity product:
        #
        #         T        T      T
        #  pXsurf     pXpt     pI
        #  ------   -------   ----
        #  pXdv      pXsurf    pXpt
        #
        # Where I is the objective, Xpt are the externally coordinates
        # supplied in addPointSet, Xsurf are the coordinates of the
        # VSP plot3d File and Xdv are the design variables.
        nPts = len(self.pts)

        # Extract just the single dIdpt we are working with. Make
        # a copy because we may need to modify it.

        n = 0
        for IC in self.intersectComps:
            n += len(IC.seam)

        dIdSeam = numpy.zeros((N, n*3))
        for i in range(N):
            n = 0
            for IC in self.intersectComps:
                dIdpt_i = dIdpt[i, :, :].copy()
                seamBar = IC.sens(dIdpt_i, ptSetName)
                dIdpt[i, :, :] = dIdpt_i
                dIdSeam[i, 3*n:3*(n+len(seamBar))] = seamBar.flatten()
                n += len(IC.seam)
        dIdpt = dIdpt.reshape((dIdpt.shape[0], dIdpt.shape[1]*3))

        # Now take dIdpt_i back to the plot3D surface:
        tmp = self.pointSets[ptSetName].dPtdXsurf.T.dot(dIdpt.T)

        # Now vstack the result with seamBar as that is far as the
        # forward FD jacobian went.
        tmp = numpy.vstack([tmp, dIdSeam.T])

        # Do final local jacobian transpose product back to the
        # DVs. This is a little trickier since we have a distributed
        # total surface jacobian. First we all reduce the "tmp"
        # variable which is pI/pXsurf. The "Reduce" is important
        # here...that is performance critical. We need the fast numpy
        # version. 
        
        if comm:
            tmp2 = numpy.zeros_like(tmp)
            comm.Reduce(tmp, tmp2, op=MPI.SUM)
            comm.Bcast(tmp2)
        else:
            tmp2 = tmp
        # ---------------------------------------------------------------

        # Remember the jacobian contains the surface poitns *and* the
        # the seam nodes. dIdx_compact is the final derivative for the
        # variables this proc owns. 
        dIdx_compact = self.jac.T.dot(tmp2)

        # Now scatter the dIdx_compact into dIdx_local before we all
        # reduce over the actual DVs. 
        dvKeys = list(self.DVs.keys())
        i = 0 # Counter on local Jac
        for iDV in range(len(dvKeys)):
            # I have to do this one.
            if iDV % self.comm.size == self.comm.rank:
                dIdx_local[:, iDV] = dIdx_compact[i, :]
                i += 1

        dIdx = self.comm.allreduce(dIdx_local, op=MPI.SUM)

        # Now convert to dict:
        dIdxDict = {}
        i = 0
        for dvName in self.DVs:
            dIdxDict[dvName] = numpy.array([dIdx[:, i]]).T
            i += 1

        return dIdxDict


    def addVariable(self, component, group, parm, value=None,
                    lower=None, upper=None, scale=1.0, scaledStep=True,
                    dh=1e-6):
        """
        Add a design variable definition.

        Parameters
        ----------
        component : str
            Name of the VSP component
        group : str
            Name of the VSP group
        parm : str
            Name of the VSP parameter
        value : float or None
            The design variable. If this value is not supplied (None), then
            the current value in the VSP model will be queried and used
        lower : float or None
            Lower bound for the design variable. Use None for no lower bound
        upper : float or None
            Upper bound for the design variable. Use None for no upper bound
        scale : float
            Scale factor
        scaledStep : bool
            Flag to use a scaled step sized based on the initial value of the
            variable. It will remain constant thereafter.
        dh : float
            Step size. When scaledStep is True, the actual step is dh*value. Otherwise
            this actual step is used.
        """

        container_id = vsp.FindContainer(component, 0)
        if container_id == "":
            raise Error('Bad component for DV: %s'%component)

        parm_id = vsp.FindParm(container_id, parm, group)
        if parm_id == "":
            raise Error('Bad group or parm: %s %s'%(group, parm))

        # Now we know the parmID is ok. So we just get the value
        val = vsp.GetParmVal(parm_id)

        dvName = '%s:%s:%s'%(component, group, parm)

        if value is None:
            value = val

        if scaledStep:
            dh = dh * value

            if value == 0:
                raise Error('Initial value is exactly 0. scaledStep option cannot be used'
                            'Specify an explicit dh with scaledStep=False')

        self.DVs[dvName] = vspDV(parm_id, component, group, parm, value, lower,
                                 upper, scale, dh)

    def addVariablesPyOpt(self, optProb):
        """
        Add the current set of variables to the optProb object.

        Parameters
        ----------
        optProb : pyOpt_optimization class
            Optimization problem definition to which variables are added
        """

        for dvName in self.DVs:
            dv = self.DVs[dvName]
            optProb.addVar(dvName, 'c', value=dv.value, lower=dv.lower,
                           upper=dv.upper, scale=dv.scale)

    def printDesignVariables(self):
        """
        Print a formatted list of design variables to the screen
        """
        print ('-'*85)
        print ('%30s%20s%20s%15s'%('Component','Group','Parm','Value'))
        print ('-'*85)
        for dvName in self.DVs:
            DV = self.DVs[dvName]
            print('%30s%20s%20s%15g'%(DV.component, DV.group, DV.parm, DV.value))

    def createDesignFile(self, fileName):
        """Take the current set of design variables and create a .des file"""
        f = open(fileName, 'w')
        f.write('%d\n'%len(self.DVs))
        for dvName in self.DVs:
            DV = self.DVs[dvName]
            f.write('%s:%s:%s:%s:%20.15g\n'%(DV.parmID, DV.component, DV.group, DV.parm, DV.value))
        f.close()

    def writePlot3D(self, fileName):
        """Write the current design to Plot3D file"""

        for dvName in self.DVs:
            DV = self.DVs[dvName]
            # We use float here since sometimes pyoptsparse will give
            # stupid numpy zero-dimensional arrays, which swig does
            # not like.
            vsp.SetParmVal(DV.parmID, float(DV.value))
        vsp.Update()

        # Write the export file.
        vsp.ExportFile(fileName, self.exportSet, vsp.EXPORT_PLOT3D)

# ----------------------------------------------------------------------
#        THE REMAINDER OF THE FUNCTIONS NEED NOT BE CALLED BY THE USER
# ----------------------------------------------------------------------

    def _getUpdatedSurface(self, fName='export.x'):
        """Return the updated surface for the currently set design variables,
        converted to an unstructured format. Note that this routine is
        safe to call in an embarrassing parallel fashion. That is it
        can be called with different DVs set on different
        processors and they won't interfere with each other.
        """
        import time

        # Set each of the DVs. We have the parmID stored so its easy.
        for dvName in self.DVs:
            DV = self.DVs[dvName]
            # We use float here since sometimes pyoptsparse will give
            # stupid numpy zero-dimensional arrays, which swig does
            # not like.
            vsp.SetParmVal(DV.parmID, float(DV.value))
        vsp.Update()

        # Write the export file.
        fName = os.path.join(self.tmpDir, fName)
        vsp.ExportFile(fName, self.exportSet, vsp.EXPORT_PLOT3D)

        # Now we load in this updated plot3d file and compute the
        # connectivity if not already done.
        computeConn = False
        if self.conn is None:
            computeConn = True
        pts, conn, cumSizes, sizes = self._readPlot3DSurfFile(fName, computeConn)

        # Apply the VSP scale
        newPts = pts * self.vspScale

        return newPts, conn, cumSizes, sizes

    def _readPlot3DSurfFile(self, fileName, computeConn=False):
        """Read a scii plot3d file and return the points and connectivity"""

        f = open(fileName, 'r')
        nSurf = numpy.fromfile(f, 'int', count=1, sep=' ')[0]
        sizes = numpy.fromfile(f, 'int', count=3*nSurf, sep=' ').reshape((nSurf, 3))
        nPts = 0; nElem = 0
        cumSizes = numpy.zeros(nSurf + 1, 'intc')
        for i in range(nSurf):
            curSize = sizes[i, 0]*sizes[i, 1]
            nPts += curSize
            cumSizes[i+1] = cumSizes[i] + curSize
            nElem += (sizes[i, 0]-1)*(sizes[i, 1]-1)

        # Generate the uncompacted point and connectivity list:
        pts = numpy.zeros((nPts, 3))
        if computeConn:
            conn = numpy.zeros((nElem, 4), dtype='intc')
        else:
            conn = self.conn

        nodeCount = 0
        elemCount = 0
        for iSurf in range(nSurf):
            curSize = sizes[iSurf, 0]*sizes[iSurf, 1]
            for idim in range(3):
                pts[nodeCount:nodeCount+curSize, idim] = (
                    numpy.fromfile(f, 'float', curSize, sep=' '))
            if computeConn:
                # Add in the connectivity.
                iSize = sizes[iSurf, 0]
                for j in range(sizes[iSurf, 1]-1):
                    for i in range(sizes[iSurf, 0]-1):
                        conn[elemCount, 0] = nodeCount + j*iSize + i
                        conn[elemCount, 1] = nodeCount + j*iSize + i+1
                        conn[elemCount, 2] = nodeCount + (j+1)*iSize + i+1
                        conn[elemCount, 3] = nodeCount + (j+1)*iSize + i
                        elemCount += 1
            nodeCount += curSize

        return pts, conn, cumSizes, sizes

    def __del__(self):
        # if self.comm.rank == 0 and not self.debug:
        #     if os.path.exists(self.tmpDir):
        #         shuitl.rmtree(self.tmpDir)
        pass

    def _computeSurfJacobian(self):
        """This routine comptues the jacobian of the VSP surface with respect
        to the design variables. Since this can be somewhat costly, we
        can compute it in an embarassingly parallel fashion across the
        COMM this object was created on.
        """

        # Determine the size of the 'surface' jacobian. It actually
        # contains the surface from VSP as well as any requried
        # intersection curves.

        ptsRef = self.pts.flatten()
        nPts = len(self.pts)

        seamsRef = numpy.zeros((0, 3))
        for IC in self.intersectComps:
            seamsRef = numpy.vstack([seamsRef, IC.seam])
        nSeam = len(seamsRef)
        seamsRef = seamsRef.flatten()

        dvKeys = list(self.DVs.keys())

        # Determine how many variables we'll perturb. 
        n = 0
        for iDV in range(len(dvKeys)):
            # I have to do this one.
            if iDV % self.comm.size == self.comm.rank:
                n += 1

        # The section of the jacobian that I own:
        self.jac = numpy.zeros((3*(nPts + nSeam), n))

        i = 0 # Counter on local Jac
        reqs = []
        for iDV in range(len(dvKeys)):
            # I have to do this one.
            if iDV % self.comm.size == self.comm.rank:

                # Step size for this particular DV
                dh =  self.DVs[dvKeys[iDV]].dh

                # Perturb the DV
                dvSave = self.DVs[dvKeys[iDV]].value
                self.DVs[dvKeys[iDV]].value += dh

                # Get the updated points. Note that this is the only
                # time we need to call _getUpdatedSurface in parallel,
                # so we have to supply a unique file name.
                pts, conn, cumSizes, sizes = self._getUpdatedSurface('dv_perturb_%d.x'%iDV)
                self.jac[0:nPts*3, i] = (pts.flatten() - ptsRef)/dh

                # Do any required intersections:
                seams = numpy.zeros((0, 3))
                for IC in self.intersectComps:
                    IC.setSurface(pts)
                    seams = numpy.vstack([seams, IC.seam])

                self.jac[nPts*3:, i] = (seams.flatten() - seamsRef)/dh

                # Reset the DV
                self.DVs[dvKeys[iDV]].value = dvSave

                i += 1

        # Restore the seams
        for IC in self.intersectComps:
            IC.setSurface(self.pts)


class vspDV(object):

    def __init__(self, parmID, component, group, parm, value, lower, upper, scale, dh):
        """Inernal class for storing VSP design variable information"""
        self.parmID = parmID
        self.component = component
        self.group = group
        self.parm = parm
        self.value = value
        self.lower = lower
        self.upper = upper
        self.dh = dh
        self.scale = scale

class PointSet(object):
    def __init__(self, points, faceID, uv, offset, conn, nSurf):
        self.points = points
        self.faceID = faceID
        self.uv = uv
        self.offset = offset

        # We need to compute the (constant) dPtdXsurf. This is the
        # sensivitiy of the the points of this proc set with respect
        # to the positions on the VSP surface (Xsurf). This is
        # computed directly from the (linear) shape functions, so we
        # have all the information we need. Since this is (very)
        # sparse, we'll use a scipy sparse matrix.

        jac = sparse.lil_matrix((len(self.points)*3, nSurf*3))

        # Loop over each point:
        for i in range(len(self.points)):

            # This point depends on the the 4 nodes of its quad.
            for idim in range(3):
                jac[i*3+idim, 3*conn[self.faceID[i], 0]+idim] = (1 - uv[i, 0])*(1 - uv[i, 1])
                jac[i*3+idim, 3*conn[self.faceID[i], 1]+idim] = (    uv[i, 0])*(1 - uv[i, 1])
                jac[i*3+idim, 3*conn[self.faceID[i], 2]+idim] = (    uv[i, 0])*(    uv[i, 1])
                jac[i*3+idim, 3*conn[self.faceID[i], 3]+idim] = (1 - uv[i, 0])*(    uv[i, 1])

        # Convert to csc becuase we'll be doing a transpose product on it.
        self.dPtdXsurf = jac.tocsc()


class CompIntersection(object):
    def __init__(self, compA, compB, extraComps, direction, dStar, pts, cumSizes, sizes, tmpDir):
        '''Class to store information required for an intersection.  The order
        of compA and compB are important: LSect will give a different
        result for compA intersecting compB than compB intersecting
        compA.

        Input
        -----
        compA , int : Index of the surface in the plot3D file
        compB , int : Index of the surface in the plot3D file
        extraComps, list : Indexes of other comps to move as well 
        direction, str: Coordiante index direction to use for
            lsect. Must be 'j' or 'k'
        dStar, real : Radius over which to attenuate the deformation
        pts , real (N, 3) : All coordinates that are being used for this model

        Internally we will store the indices and the weights of the
        points that this intersection will have to modify. In general,
        all this code is not super efficient since it's all python,
        but it should not be more of a bottleneck than VSP is itself
        in doing the export.
        '''
        self.compA = compA
        self.compB = compB
        self.extraComps = extraComps
        self.dir = direction.lower()
        self.dStar = dStar
        self.halfdStar = dStar/2.0
        self.pts = pts
        self.cumSizes = cumSizes
        self.sizes = sizes
        if direction not in ['j', 'k']:
            raise Error("Direction must be given as 'j' or 'k' for intersection")
        if direction == 'j':
            self.dir = 1
        else:
            self.dir = 2

        # First generate the initial intersection curve.
        s = self.sizes[self.compA]
        blkA = self.pts[self.cumSizes[self.compA]:self.cumSizes[self.compA+1], :].reshape(
            (s[0], s[1], s[2], 3), order='f')

        s = self.sizes[self.compB]
        blkB = self.pts[self.cumSizes[self.compB]:self.cumSizes[self.compB+1], :].reshape(
            (s[0], s[1], s[2], 3), order='f')

        nmax = max(blkB.shape[0], blkB.shape[1], blkA.shape[0], blkA.shape[1])
        seam, nOut = lsect.lsect_wrap(blkA, blkB, self.dir, nmax)
        self.seam0 = seam[0:nOut, :]
        self.seam = self.seam0.copy()
        self.ptSets = {}

    def setSurface(self, pts):
        """ This set the new udpated surface on which we need to comptue the new intersection curve"""

        s = self.sizes[self.compA]
        blkA = pts[self.cumSizes[self.compA]:self.cumSizes[self.compA+1], :].reshape(
            (s[0], s[1], s[2], 3), order='f')

        s = self.sizes[self.compB]
        blkB = pts[self.cumSizes[self.compB]:self.cumSizes[self.compB+1], :].reshape(
            (s[0], s[1], s[2], 3), order='f')

        nmax = max(blkB.shape[0], blkB.shape[1], blkA.shape[0], blkA.shape[1])
        seam, nOut = lsect.lsect_wrap(blkA, blkB, self.dir, nmax)
        self.seam = seam[0:nOut, :]

    def addPointSet(self, pts, compIDs, ptSetName):

        # Figure out which points this intersection object has to deal with.
        tree = cKDTree(self.seam0)

        # Find all the distances. This is pretty fast.
        d, index = tree.query(pts)
        indices = []
        factors = []
        for i in range(len(pts)):
            if compIDs[i] == self.compA or compIDs[i] == self.compB or compIDs[i] in self.extraComps:
                if d[i] < self.dStar:

                    # Compute the factor
                    if d[i] < self.halfdStar:
                        factor = .5*(d[i]/self.halfdStar)**3
                    else:
                        factor = .5*(2-((self.dStar - d[i])/self.halfdStar)**3)

                    # Save the index and factor
                    indices.append(i)
                    factors.append(factor)

        # Save the affected indices and the factor in the little dictionary
        self.ptSets[ptSetName] = [pts.copy(), indices, factors]

    def update(self, ptSetName, delta):

        """Update the delta in ptSetName with our correction. The delta need
        to be supplied as we will be changing it and returning them
        """
        pts     = self.ptSets[ptSetName][0]
        indices = self.ptSets[ptSetName][1]
        factors = self.ptSets[ptSetName][2]
        seamDiff = self.seam - self.seam0
        for i in range(len(factors)):

            # j is the index of the point in the full set we are
            # working with.
            j = indices[i]

            #Run the weighted interp:
            # num = numpy.zeros(3)
            # den = 0.0
            # for k in range(len(seamDiff)):
            #     rr = pts[j] - self.seam0[k]
            #     LdefoDist = 1.0/numpy.sqrt(rr[0]**2 + rr[1]**2 + rr[2]**2+1e-16)
            #     LdefoDist3 = LdefoDist**3
            #     Wi = LdefoDist3
            #     Si = seamDiff[k]
            #     num = num + Wi*Si
            #     den = den + Wi

            # interp = num / den

            # Do it vectorized
            rr = pts[j] - self.seam0
            LdefoDist = (1.0/numpy.sqrt(rr[:,0]**2 + rr[:,1]**2 + rr[:,2]**2+1e-16))
            LdefoDist3 = LdefoDist**3
            Wi = LdefoDist3
            den = numpy.sum(Wi)
            interp = numpy.zeros(3)
            for iDim in range(3):
                interp[iDim] = numpy.sum(Wi*seamDiff[:, iDim])/den

            # Now the delta is replaced by 1-factor times the weighted
            # interp of the seam * factor of the original:

            delta[j] = factors[i]*delta[j] + (1-factors[i])*interp

        return delta

    def sens(self, dIdPt, ptSetName):
        # Return the reverse accumulation of dIdpt on the seam
        # nodes. Also modifies the dIdp array accordingly.

        pts     = self.ptSets[ptSetName][0]
        indices = self.ptSets[ptSetName][1]
        factors = self.ptSets[ptSetName][2]

        seamBar = numpy.zeros_like(self.seam0)
        for i in range(len(factors)):

            # j is the index of the point in the full set we are
            # working with.
            j = indices[i]

            # This is the local seed (well the 3 seeds for the point)
            localVal = dIdPt[j]*(1 - factors[i])

            # Scale the dIdpt by the factor..dIdpt is input/output
            dIdPt[j] *= factors[i]

            # Do it vectorized
            rr = pts[j] - self.seam0
            LdefoDist = (1.0/numpy.sqrt(rr[:,0]**2 + rr[:,1]**2 + rr[:,2]**2+1e-16))
            LdefoDist3 = LdefoDist**3
            Wi = LdefoDist3
            den = numpy.sum(Wi)
            interp = numpy.zeros(3)
            for iDim in range(3):
                seamBar[:, iDim] += Wi*localVal[iDim]/den

        return seamBar