# Source code for pygeo.parameterization.DVGeoMulti

# ======================================================================
#         Imports
# ======================================================================
from collections import OrderedDict
import numpy as np
from mpi4py import MPI
from baseclasses.utils import Error
from pysurf import intersectionAPI, curveSearchAPI, utilitiesAPI, adtAPI, tsurf_tools, tecplot_interface
from pysurf import intersectionAPI_cs, curveSearchAPI_cs, utilitiesAPI_cs, adtAPI_cs

[docs]class DVGeometryMulti:
"""
A class for manipulating multiple components using multiple FFDs
and handling design changes near component intersections.

Parameters
----------
comm : MPI.IntraComm, optional
The communicator associated with this geometry object.

checkDVs : bool, optional
Flag to check whether there are duplicate DV names in or across components.

debug : bool, optional
Flag to generate output useful for debugging the intersection setup.

isComplex : bool, optional
Flag to use complex variables for complex step verification.

"""

def __init__(self, comm=MPI.COMM_WORLD, checkDVs=True, debug=False, isComplex=False):

self.compNames = []
self.comps = OrderedDict()
self.DVGeoDict = OrderedDict()
self.points = OrderedDict()
self.comm = comm
self.updated = {}
self.intersectComps = []
self.checkDVs = checkDVs
self.debug = debug
self.complex = isComplex

# Set real or complex Fortran API
if isComplex:
self.dtype = complex
else:
self.dtype = float

[docs]    def addComponent(self, comp, DVGeo, triMesh=None, scale=1.0, bbox=None, pointSetKwargs=None):
"""
Method to add components to the DVGeometryMulti object.

Parameters
----------
comp : str
The name of the component.

DVGeo : DVGeometry
The DVGeometry object defining the component FFD.

triMesh : str, optional
Path to the triangulated mesh file for this component.

scale : float, optional
A multiplicative scaling factor applied to the triangulated mesh coordinates.
Useful for when the scales of the triangulated and CFD meshes do not match.

bbox : dict, optional
Specify a bounding box that is different from the bounds of the FFD.
The keys can include xmin, xmax, ymin, ymax, zmin, zmax.
If any of these are not provided, the FFD bound is used.

pointSetKwargs : dict, optional
Keyword arguments to be passed to the component addPointSet call for the triangulated mesh.

"""

# Assign mutable defaults
if bbox is None:
bbox = {}
if pointSetKwargs is None:
pointSetKwargs = {}

if triMesh is not None:
# We also need to read the triMesh and save the points
nodes, triConn, triConnStack, barsConn = self._readCGNSFile(triMesh)

# scale the nodes
nodes *= scale

# add these points to the corresponding dvgeo
else:
# the user has not provided a triangulated surface mesh for this file
nodes = None
triConn = None
triConnStack = None
barsConn = None

# we will need the bounding box information later on, so save this here
xMin, xMax = DVGeo.FFD.getBounds()

# also we might want to modify the bounding box if the user specified any coordinates
if "xmin" in bbox:
xMin[0] = bbox["xmin"]
if "ymin" in bbox:
xMin[1] = bbox["ymin"]
if "zmin" in bbox:
xMin[2] = bbox["zmin"]
if "xmax" in bbox:
xMax[0] = bbox["xmax"]
if "ymax" in bbox:
xMax[1] = bbox["ymax"]
if "zmax" in bbox:
xMax[2] = bbox["zmax"]

# initialize the component object
self.comps[comp] = component(comp, DVGeo, nodes, triConn, triConnStack, barsConn, xMin, xMax)

# add the name to the list
self.compNames.append(comp)

# also save the DVGeometry pointer in the dictionary we pass back
self.DVGeoDict[comp] = DVGeo

self,
compA,
compB,
dStarA=0.2,
dStarB=0.2,
featureCurves=None,
distTol=1e-14,
project=False,
marchDir=1,
includeCurves=False,
intDir=None,
curveEpsDict=None,
trackSurfaces=None,
excludeSurfaces=None,
remeshBwd=True,
):
"""
Method that defines intersections between components.

Parameters
----------
compA : str
The name of the first component.

compB : str
The name of the second component.

dStarA : float, optional
Distance from the intersection over which the inverse-distance deformation is applied on compA.

dStarB : float, optional
Distance from the intersection over which the inverse-distance deformation is applied on compB.

featureCurves : list or dict, optional
Points on feature curves will remain on the same curve after deformations and projections.
Feature curves can be specified as a list of curve names.
In this case, the march direction for all curves is marchDir.
Alternatively, a dictionary can be provided.
In this case, the keys are the curve names and the values are the march directions for each curve.
See marchDir for the definition of march direction.

distTol : float, optional
Distance tolerance to merge nearby nodes in the intersection curve.

project : bool, optional
Flag to specify whether to project points to curves and surfaces after the deformation step.

marchDir : int, optional
The side of the intersection where the feature curves are remeshed.
The sign determines the direction and the value (1, 2, 3) specifies the axis (x, y, z).
If remeshBwd is True, the other side is also remeshed.
In this case, the march direction only serves to define the 'free end' of the feature curve.
If None, the entire curve is remeshed.
This argument is only used if a list is provided for featureCurves.

includeCurves : bool, optional
Flag to specify whether to include features curves in the inverse-distance deformation.

intDir : int, optional
If there are multiple intersection curves, this specifies which curve to choose.
The sign determines the direction and the value (1, 2, 3) specifies the axis (x, y, z).
For example, -1 specifies the intersection curve as the one that is further in the negative x-direction.

curveEpsDict : dict, optional
Required if using feature curves.
The keys of the dictionary are the curve names and the values are distances.
All points within the specified distance from the curve are considered to be on the curve.

trackSurfaces : dict, optional
Points on tracked surfaces will remain on the same surfaces after deformations and projections.
The keys of the dictionary are the surface names and the values are distances.
All points within the specified distance from the surface are considered to be on the surface.

excludeSurfaces : dict, optional
Points on excluded surfaces are removed from the intersection computations.
The keys of the dictionary are the surface names and the values are distances.
All points within the specified distance from the surface are considered to be on the surface.

remeshBwd : bool, optional
Flag to specify whether to remesh feature curves on the side opposite that
which is specified by the march direction.

"""

# Assign mutable defaults
if featureCurves is None:
featureCurves = []
if curveEpsDict is None:
curveEpsDict = {}
if trackSurfaces is None:
trackSurfaces = {}
if excludeSurfaces is None:
excludeSurfaces = {}

# just initialize the intersection object
self.intersectComps.append(
CompIntersection(
compA,
compB,
dStarA,
dStarB,
featureCurves,
distTol,
self,
project,
marchDir,
includeCurves,
intDir,
curveEpsDict,
trackSurfaces,
excludeSurfaces,
remeshBwd,
self.debug,
self.dtype,
)
)

[docs]    def getDVGeoDict(self):
"""Return a dictionary of component DVGeo objects."""
return self.DVGeoDict

[docs]    def addPointSet(self, points, ptName, compNames=None, comm=None, applyIC=False, **kwargs):
"""
Add a set of coordinates to DVGeometryMulti.
The is the main way that geometry, in the form of a coordinate list, is manipulated.

Parameters
----------
points : array, size (N,3)
The coordinates to embed.
These coordinates should all be inside at least one FFD volume.
ptName : str
A user supplied name to associate with the set of coordinates.
This name will need to be provided when updating the coordinates
or when getting the derivatives of the coordinates.
compNames : list, optional
A list of component names that this point set should be added to.
To ease bookkeepping, an empty point set with ptName will be added to components not in this list.
If a list is not provided, this point set is added to all components.
comm : MPI.IntraComm, optional
Comm that is associated with the added point set. Does not
work now, just added to be consistent with the API of
other DVGeo types.
applyIC : bool, optional
Flag to specify whether this point set will follow the updated intersection curve(s).
This is typically only needed for the CFD surface mesh.

"""

# if compList is not provided, we use all components
if compNames is None:
compNames = self.compNames

# before we do anything, we need to create surface ADTs
# for which the user provided triangulated meshes
for comp in compNames:

# check if we have a trimesh for this component
if self.comps[comp].triMesh:

# Now we build the ADT using pySurf
# Set bounding box for new tree
BBox = np.zeros((2, 3))
useBBox = False

# dummy connectivity data for quad elements since we have all tris

# Compute set of nodal normals by taking the average normal of all
# elements surrounding the node. This allows the meshing algorithms,
# for instance, to march in an average direction near kinks.
)
self.comps[comp].nodal_normals = nodal_normals.T

# Create new tree (the tree itself is stored in Fortran level)
self.comps[comp].nodes.T,
self.comps[comp].triConnStack.T,
BBox.T,
useBBox,
MPI.COMM_SELF.py2f(),
comp,
)

# create the pointset class
self.points[ptName] = PointSet(points, comm=comm)

for comp in self.compNames:
# initialize the list for this component
self.points[ptName].compMap[comp] = []
self.points[ptName].compMapFlat[comp] = []

# we now need to create the component mapping information
for i in range(self.points[ptName].nPts):

# initial flags
inFFD = False
proj = False
projList = []

# loop over components and check if this point is in a single BBox
for comp in compNames:

# apply a small tolerance for the bounding box in case points are coincident with the FFD
boundTol = 1e-16
xMin = self.comps[comp].xMin
xMax = self.comps[comp].xMax
xMin -= np.abs(xMin * boundTol) + boundTol
xMax += np.abs(xMax * boundTol) + boundTol

# check if inside
if (
xMin[0] < points[i, 0] < xMax[0]
and xMin[1] < points[i, 1] < xMax[1]
and xMin[2] < points[i, 2] < xMax[2]
):

# add this component to the projection list
projList.append(comp)

# this point was not inside any other FFD before
if not inFFD:
inFFD = True
inComp = comp
# this point was inside another FFD, so we need to project it...
else:
# set the projection flag
proj = True

# project this point to components, we need to set inComp string
if proj:

# set a high initial distance
dMin2 = 1e10

# loop over the components
for comp in compNames:
# check if this component is in the projList
if comp in projList:

# check if we have an ADT:
if self.comps[comp].triMesh:
# Initialize reference values (see explanation above)
numPts = 1
dist2 = np.ones(numPts, dtype=self.dtype) * 1e10
xyzProj = np.zeros((numPts, 3), dtype=self.dtype)
normProjNotNorm = np.zeros((numPts, 3), dtype=self.dtype)

# Call projection function
points[i].T, comp, dist2, xyzProj.T, self.comps[comp].nodal_normals.T, normProjNotNorm.T
)

# if this is closer than the previous min, take this comp
if dist2 < dMin2:
dMin2 = dist2[0]
inComp = comp

else:
raise Error(
f"The point at (x, y, z) = ({points[i, 0]:.3f}, {points[i, 1]:.3f} {points[i, 2]:.3f})"
+ f"in point set {ptName} is inside multiple FFDs but a triangulated mesh "
+ f"for component {comp} is not provided to determine which component owns this point."
)

# this point was inside at least one FFD. If it was inside multiple,
# we projected it before to figure out which component it should belong to
if inFFD:
# we can add the point index to the list of points inComp owns
self.points[ptName].compMap[inComp].append(i)

# also create a flattened version of the compMap
for j in range(3):
self.points[ptName].compMapFlat[inComp].append(3 * i + j)

# this point is outside any FFD...
else:
raise Error(
f"The point at (x, y, z) = ({points[i, 0]:.3f}, {points[i, 1]:.3f} {points[i, 2]:.3f}) "
+ f"in point set {ptName} is not inside any FFDs."
)

# using the mapping array, add the pointsets to respective DVGeo objects
for comp in self.compNames:
compMap = self.points[ptName].compMap[comp]

# check if this pointset will get the IC treatment
if applyIC:
# loop over the intersections and add pointsets
for IC in self.intersectComps:

# finally, we can deallocate the ADTs
for comp in compNames:
if self.comps[comp].triMesh:

# mark this pointset as up to date
self.updated[ptName] = False

[docs]    def setDesignVars(self, dvDict):
"""
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 dictionary are simply ignored.

"""

# Check if we have duplicate DV names
if self.checkDVs:
dvNames = self.getVarNames()
duplicates = len(dvNames) != len(set(dvNames))
if duplicates:
raise Error(
"There are duplicate DV names in a component or across components. "
"If this is intended, initialize the DVGeometryMulti class with checkDVs=False."
)

# loop over the components and set the values
for comp in self.compNames:
self.comps[comp].DVGeo.setDesignVars(dvDict)

# 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.comm)

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

[docs]    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 setDesignVars().

Returns
-------
dvDict : dict
Dictionary of design variables.

"""

dvDict = {}
# we need to loop over each DVGeo object and get the DVs
for comp in self.compNames:
dvDictComp = self.comps[comp].DVGeo.getValues()
# we need to loop over these DVs
for k, v in dvDictComp.items():
dvDict[k] = v

return dvDict

[docs]    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 one of those added in an :func:addPointSet() call.

"""

# get the new points
newPts = np.zeros((self.points[ptSetName].nPts, 3), dtype=self.dtype)

# we first need to update all points with their respective DVGeo objects
for comp in self.compNames:
ptsComp = self.comps[comp].DVGeo.update(ptSetName)

# now save this info with the pointset mapping
ptMap = self.points[ptSetName].compMap[comp]
newPts[ptMap] = ptsComp

# get the delta
delta = newPts - self.points[ptSetName].points

# then apply the intersection treatment
for IC in self.intersectComps:
# check if this IC is active for this ptSet
if ptSetName in IC.points:
delta = IC.update(ptSetName, delta)

# now we are ready to take the delta which may be modified by the intersections
newPts = self.points[ptSetName].points + delta

# now, project the points that were warped back onto the trimesh
for IC in self.intersectComps:
if IC.projectFlag and ptSetName in IC.points:
# new points will be modified in place using the newPts array
IC.project(ptSetName, newPts)

# set the pointset up to date
self.updated[ptSetName] = True

return newPts

[docs]    def pointSetUpToDate(self, ptSetName):
"""
This is used externally to query if the object needs to update its point set or not.
When update() is called with a point set, the self.updated value for pointSet is flagged as True.
We reset all flags to False when design variables are set because nothing (in general) will up to date anymore.
Here we just return that flag.

Parameters
----------
ptSetName : str
The name of the pointset to check.

"""
if ptSetName in self.updated:
return self.updated[ptSetName]
else:
return True

[docs]    def getNDV(self):
"""Return the number of DVs."""
# Loop over components and sum the number of DVs
nDV = 0
for comp in self.compNames:
nDV += self.comps[comp].DVGeo.getNDV()
return nDV

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

Examples
--------

"""
dvNames = []
# create a list of DVs from each comp
for comp in self.compNames:
# first get the list of DVs from this component
varNames = self.comps[comp].DVGeo.getVarNames()

# add the component DVs to the full list
dvNames.extend(varNames)

return dvNames

[docs]    def totalSensitivity(self, dIdpt, ptSetName, comm=None, config=None):
"""
This function computes sensitivity 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, optional
The communicator to use to reduce the final derivative.
If comm is None, no reduction takes place.

config : str or list, optional
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 appies to *ALL* configurations.

Returns
-------
dIdxDict : dict
The dictionary containing the derivatives, suitable for pyOptSparse.

Notes
-----
The child and nDVStore options are only used
internally and should not be changed by the user.

"""

# Compute the total Jacobian for this point set
# TODO: this can be done without converting sparse matrices to dense
# or with a function for mat-vec products with all DVGeos
self._computeTotalJacobian(ptSetName)

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

# create a dictionary to save total sensitivity info that might come out of the ICs
compSensList = []

# if we projected points for any intersection treatment,
# we need to propagate the derivative seed of the projected points
# back to the seeds for the initial points we get after ID-warping
for IC in self.intersectComps:
if IC.projectFlag and ptSetName in IC.points:
# initialize the seed contribution to the intersection seam and feature curves from project_b
IC.seamBarProj[ptSetName] = np.zeros((N, IC.seam0.shape[0], IC.seam0.shape[1]))

# we pass in dIdpt and the intersection object, along with pointset information
# the intersection object adjusts the entries corresponding to projected points
# and passes back dIdpt in place.
compSens = IC.project_b(ptSetName, dIdpt, comm)

# append this to the dictionary list...
compSensList.append(compSens)

# do the transpose multiplication

if self.debug:
print(f"[{self.comm.rank}] finished project_b")

# we need to go through all ICs bec even though some procs might not have points on the intersection,
# communication is easier and we can reduce compSens as we compute them
for IC in self.intersectComps:
if ptSetName in IC.points:
compSens = IC.sens(dIdpt, ptSetName, comm)
# save the sensitivities from the intersection stuff
compSensList.append(compSens)

if self.debug:
print(f"[{self.comm.rank}] finished IC.sens")

# reshape the dIdpt array from [N] * [nPt] * [3] to  [N] * [nPt*3]
dIdpt = dIdpt.reshape((dIdpt.shape[0], dIdpt.shape[1] * 3))

# jacobian for the pointset
jac = self.points[ptSetName].jac

# this is the mat-vec product for the remaining seeds.
# this only contains the effects of the FFD motion,
# projections and intersections are handled separately in compSens
dIdxT_local = jac.T.dot(dIdpt.T)
dIdx_local = dIdxT_local.T

# If we have a comm, globaly reduce with sum
if comm:
dIdx = comm.allreduce(dIdx_local, op=MPI.SUM)
else:
dIdx = dIdx_local

# use respective DVGeo's convert to dict functionality
dIdxDict = OrderedDict()
dvOffset = 0
for comp in self.compNames:
DVGeo = self.comps[comp].DVGeo
nDVComp = DVGeo.getNDV()

# we only do this if this component has at least one DV
if nDVComp > 0:
# this part of the sensitivity matrix is owned by this dvgeo
dIdxComp = DVGeo.convertSensitivityToDict(dIdx[:, dvOffset : dvOffset + nDVComp])

for k, v in dIdxComp.items():
dIdxDict[k] = v

# also increment the offset
dvOffset += nDVComp

# finally, we can add the contributions from triangulated component meshes
for compSens in compSensList:
# loop over the items of compSens, which are guaranteed to be in dIdxDict
for k, v in compSens.items():
# these will bring in effects from projections and intersection computations
dIdxDict[k] += v

if self.debug:
print(f"[{self.comm.rank}] finished DVGeo.totalSensitivity")

return dIdxDict

self,
optProb,
globalVars=True,
localVars=True,
sectionlocalVars=True,
ignoreVars=None,
freezeVars=None,
comps=None,
):
"""
Add the current set of variables to the optProb object.

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

globalVars : bool
Flag specifying whether global variables are to be added

localVars : bool
Flag specifying whether local variables are to be added

ignoreVars : list of strings
List of design variables the user DOESN'T want to use
as optimization variables.

freezeVars : list of string
List of design variables the user WANTS to add as optimization
variables, but to have the lower and upper bounds set at the current
variable. This effectively eliminates the variable, but it the variable
is still part of the optimization.

comps : list
List of components we want to add the DVs of.
If no list is provided, we will add DVs from all components.

"""

# If no list was provided, we use all components
if comps is None:
comps = self.compNames

# We can simply loop over all DV objects and call their respective addVariablesPyOpt function
for comp in comps:
optProb,
globalVars=globalVars,
localVars=localVars,
sectionlocalVars=sectionlocalVars,
ignoreVars=ignoreVars,
freezeVars=freezeVars,
)

[docs]    def getLocalIndex(self, iVol, comp):
"""Return the local index mapping that points to the global coefficient list for a given volume.

Parameters
----------

iVol : int
Index specifying the FFD volume.

comp : str
Name of the component.

"""

# Call this on the component DVGeo
DVGeo = self.comps[comp].DVGeo
return DVGeo.FFD.topo.lIndex[iVol].copy()

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

# this function reads the unstructured CGNS grid in filename and returns
# node coordinates and element connectivities.
# Here, only the root proc reads the cgns file, broadcasts node and connectivity info.

# only root proc reads the file
if self.comm.rank == 0:
# use the default routine in tsurftools
nodes, sectionDict = tsurf_tools.getCGNSsections(filename, comm=MPI.COMM_SELF)

# Convert the nodes to complex if necessary
nodes = nodes.astype(self.dtype)

triConn = {}
triConnStack = np.zeros((0, 3), dtype=np.int8)
barsConn = {}

for part in sectionDict:
if "triaConnF" in sectionDict[part].keys():
# this is a surface, read the tri connectivities
triConn[part.lower()] = sectionDict[part]["triaConnF"]
triConnStack = np.vstack((triConnStack, sectionDict[part]["triaConnF"]))

if "barsConn" in sectionDict[part].keys():
# this is a curve, save the curve connectivity
barsConn[part.lower()] = sectionDict[part]["barsConn"]

print(f"The {filename} mesh has {len(nodes)} nodes and {len(triConnStack)} elements.")
else:
# create these to recieve the data
nodes = None
triConn = None
triConnStack = None
barsConn = None

# each proc gets the nodes and connectivities
nodes = self.comm.bcast(nodes, root=0)
triConn = self.comm.bcast(triConn, root=0)
triConnStack = self.comm.bcast(triConnStack, root=0)
barsConn = self.comm.bcast(barsConn, root=0)

return nodes, triConn, triConnStack, barsConn

def _computeTotalJacobian(self, ptSetName):
"""
This routine computes the total jacobian. It takes the jacobians
from respective DVGeo objects and also computes the jacobians for
the intersection seams. We then use this information in the
totalSensitivity function.

"""

# number of design variables
nDV = self.getNDV()

# allocate space for the jacobian.
jac = np.zeros((self.points[ptSetName].nPts * 3, nDV))

# ptset
ptSet = self.points[ptSetName]

dvOffset = 0
# we need to call computeTotalJacobian from all comps and get the jacobians for this pointset
for comp in self.compNames:

# number of design variables
nDVComp = self.comps[comp].DVGeo.getNDV()

# call the function to compute the total jacobian
self.comps[comp].DVGeo.computeTotalJacobian(ptSetName)

if self.comps[comp].DVGeo.JT[ptSetName] is not None:
# We convert to dense storage
# This is not the best way to do this
# TODO: use sparse storage for these
compJ = self.comps[comp].DVGeo.JT[ptSetName].todense().T

# do it (kinda) vectorized
jac[ptSet.compMapFlat[comp], dvOffset : dvOffset + nDVComp] = compJ[:, :]

# increment the offset
dvOffset += nDVComp

# now we can save this jacobian in the pointset
ptSet.jac = jac

class component:
def __init__(self, name, DVGeo, nodes, triConn, triConnStack, barsConn, xMin, xMax):

# save the info
self.name = name
self.DVGeo = DVGeo
self.nodes = nodes
self.triConn = triConn
self.triConnStack = triConnStack
self.barsConn = barsConn
self.xMin = xMin
self.xMax = xMax

# also a dictionary for DV names
self.dvDict = {}

# set a flag for triangulated meshes
if nodes is None:
self.triMesh = False
else:
self.triMesh = True

def updateTriMesh(self):
# update the triangulated surface mesh
self.nodes = self.DVGeo.update("triMesh")

class PointSet:
def __init__(self, points, comm):
self.points = points
self.nPts = len(self.points)
self.compMap = OrderedDict()
self.compMapFlat = OrderedDict()
self.comm = comm

class CompIntersection:
def __init__(
self,
compA,
compB,
dStarA,
dStarB,
featureCurves,
distTol,
DVGeo,
project,
marchDir,
includeCurves,
intDir,
curveEpsDict,
trackSurfaces,
excludeSurfaces,
remeshBwd,
debug,
dtype,
):
"""
Class to store information required for an intersection.
Here, we use some Fortran code from pySurf.
Internally, we store the indices and weights of the points that this intersection will modify.
This code is not super efficient because it is in Python.

See the documentation for addIntersection in DVGeometryMulti for the API.

"""

# same communicator with DVGeo
self.comm = DVGeo.comm

# counter for outputting curves etc at each update
self.counter = 0

# flag that determines if we will remesh the other side of the feature curves on compB
self.remeshBwd = remeshBwd

# Flag for debug ouput
self.debug = debug

# Set real or complex Fortran APIs
self.dtype = dtype
if dtype == float:
self.curveSearchAPI = curveSearchAPI.curvesearchapi
self.intersectionAPI = intersectionAPI.intersectionapi
self.utilitiesAPI = utilitiesAPI.utilitiesapi
self.mpi_type = MPI.DOUBLE
elif dtype == complex:
self.curveSearchAPI = curveSearchAPI_cs.curvesearchapi
self.intersectionAPI = intersectionAPI_cs.intersectionapi
self.utilitiesAPI = utilitiesAPI_cs.utilitiesapi
self.mpi_type = MPI.C_DOUBLE_COMPLEX

# tolerance used for each curve when mapping nodes to curves
self.curveEpsDict = {}
for k, v in curveEpsDict.items():
self.curveEpsDict[k.lower()] = v

# beginning and end element indices for each curve
self.seamBeg = {}
self.seamEnd = {}

# indices of nodes to be projected to curves.
self.curveProjIdx = {}

# lists to track which feature curve is on which comp
self.curvesOnA = []
self.curvesOnB = []

# a dict to to keep track which curves get points mapped to them
# keys will be pointSetNames, then the values will be another dict,
# where the keys are curve names and values will be a bool value
self.curveProjFlag = {}

# dicts to keep track of the coordinates and counts of points projected to curves
self.curvePtCounts = {}
self.curvePtCoords = {}

# dict to keep track of the total number of points on each curve
self.nCurvePts = {}

# dictionary to keep the seam seeds that come from curve projections
self.seamBarProj = {}

# dictionaries to save the indices of points mapped to surfaces for each comp
self.surfIdxA = {}
self.surfIdxB = {}

# names of compA and compB must be provided
self.compA = DVGeo.comps[compA]
self.compB = DVGeo.comps[compB]

self.dStarA = dStarA
self.dStarB = dStarB
# self.halfdStar = self.dStar/2.0
self.points = OrderedDict()

# Make surface names lowercase
self.trackSurfaces = {}
for k, v in trackSurfaces.items():
self.trackSurfaces[k.lower()] = v

self.excludeSurfaces = {}
for k, v in excludeSurfaces.items():
if k.lower() in self.trackSurfaces:
raise Error(f"Surface {k} cannot be in both trackSurfaces and excludeSurfaces.")
self.excludeSurfaces[k.lower()] = v

# process the feature curves

# list to save march directions
marchDirs = []

# list to save curve names where we remesh all the curve
self.remeshAll = []

# if a list is provided, we use this and the marchdir information
if type(featureCurves) is list:
self.featureCurveNames = featureCurves
for i in range(len(self.featureCurveNames)):
self.featureCurveNames[i] = self.featureCurveNames[i].lower()
# get one march dir per curve
marchDirs.append(marchDir)

else:
# if a dict is provided, the marchdirs are the dict values
# we save this info in lists
self.featureCurveNames = []
# save the curve name and march direction information
for k, v in featureCurves.items():
self.featureCurveNames.append(k.lower())
marchDirs.append(v)

# now loop over the feature curves and flip if necessary
for ii, curveName in enumerate(self.featureCurveNames):
# figure out which comp owns this curve...
if curveName in self.compB.barsConn:
curveComp = self.compB
self.curvesOnB.append(curveName)
elif curveName in self.compA.barsConn:
curveComp = self.compA
self.curvesOnA.append(curveName)
else:
raise Error(f"Curve {curveName} does not belong in {self.compA.name} or {self.compB.name}.")

# sort the feature curve
newConn, newMap = tsurf_tools.FEsort(curveComp.barsConn[curveName].tolist())

# we only want to have a single curve
if len(newConn) > 1:
raise Error(f"The curve {curveName} generated more than one curve with FESort.")

# get the connectivity
newConn = newConn[0]

# we may also need to flip the curve
curveNodes = curveComp.nodes

if marchDirs[ii] is None:
# we remesh all of this curve
self.remeshAll.append(curveName)
curveComp.barsConn[curveName] = newConn
else:
# get the direction we want to march
mdir = abs(marchDirs[ii]) - 1
msign = np.sign(marchDirs[ii])

# check if we need to flip
if msign * curveNodes[newConn[0][0]][mdir] > msign * curveNodes[newConn[0][1]][mdir]:
# flip on both axes
newConn = np.flip(newConn, axis=0)
newConn = np.flip(newConn, axis=1)

# save the new connectivity
curveComp.barsConn[curveName] = newConn

self.distTol = distTol

# flag to determine if we want to project nodes after intersection treatment
self.projectFlag = project

# create the dictionary if we are projecting.
if project:
self.projData = {}
if includeCurves:
# dict to save the data related to projection to curves
self.curveProjData = {}

# flag to include feature curves in ID-warping
self.incCurves = includeCurves

# direction to pick if we have multiple intersection curves
self.intDir = intDir

# only the node coordinates will be modified for the intersection calculations because we have calculated and saved all the connectivity information
if self.comm.rank == 0:
print(f"Computing initial intersection between {compA} and {compB}")
self.seam0 = self._getIntersectionSeam(self.comm, firstCall=True)
self.seam = self.seam0.copy()

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

# get the updated surface coordinates
self._getUpdatedCoords()

self.seam = self._getIntersectionSeam(comm)

def addPointSet(self, pts, ptSetName, compMap, comm):

# Figure out which points this intersection object has to deal with

# Use pySurf to project the point on curve
# Get number of points
nPoints = len(pts)

# Initialize references if user provided none
dist2 = np.ones(nPoints, dtype=self.dtype) * 1e10
xyzProj = np.zeros((nPoints, 3), dtype=self.dtype)
tanProj = np.zeros((nPoints, 3), dtype=self.dtype)
elemIDs = np.zeros((nPoints), dtype="int32")

# Only call the Fortran code if we have at least one point
if nPoints > 0:
# This will modify xyzProj, tanProj, dist2, and elemIDs if we find better projections than dist2.
# Remember that we should adjust some indices before calling the Fortran code
# Remember to use [:] to don't lose the pointer (elemIDs is an input/output variable)
elemIDs[:] = (
elemIDs + 1
)  # (we need to do this separetely because Fortran will actively change elemIDs contents.
self.curveSearchAPI.mindistancecurve(
pts.T, self.nodes0.T, self.conn0.T + 1, xyzProj.T, tanProj.T, dist2, elemIDs
)

# Adjust indices back to Python standards
elemIDs[:] = elemIDs - 1

# dist2 has the array of squared distances
d = np.sqrt(dist2)

indices = []
factors = []
for i in range(len(pts)):
# figure out which component this point is mapped to
if i in compMap[self.compA.name]:
# component A owns this
dStar = self.dStarA
else:
# comp B owns this point
dStar = self.dStarB

# then get the halfdStar for that component
halfdStar = dStar / 2.0

if d[i] < dStar:

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

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

# Get all points included in the intersection computation
intersectPts = pts[indices]
nPoints = len(intersectPts)

if self.projectFlag:
# Create the dictionaries to save projection data
self.projData[ptSetName] = {
# We need one dictionary for each component
"compA": {"surfaceInd": {}},
"compB": {"surfaceInd": {}},
}

if nPoints > 0 and self.excludeSurfaces:

# Associate points with the excluded surfaces
for surface in self.excludeSurfaces:
surfaceEps = self.excludeSurfaces[surface]
self.associatePointsToSurface(intersectPts, ptSetName, surface, surfaceEps)

# Combine the excluded indices using a set to avoid duplicates
excludeSet = set()
for surface in self.excludeSurfaces:
if surface in self.compA.triConn:
# Pop this surface from the saved data
surfaceInd = self.projData[ptSetName]["compA"]["surfaceInd"].pop(surface)
elif surface in self.compB.triConn:
surfaceInd = self.projData[ptSetName]["compB"]["surfaceInd"].pop(surface)

excludeSet.update(surfaceInd)

# Invert excludeSet to get the points we want to keep
oneToN = set(range(nPoints))
includeSet = oneToN.difference(excludeSet)

# Keep only the points not associated with the excluded surfaces
indices = [indices[i] for i in includeSet]
factors = [factors[i] for i in includeSet]

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

# now we need to figure out which components we are projecting to if projection is enabled
# this can be done faster above but whatever
if self.projectFlag:
flagA = False
flagB = False

indices = self.points[ptSetName][1]

# create the list we use to map the points to projection components
indA = []
indB = []

# maybe we can do this vectorized
for ind in indices:

# check compA
if ind in compMap[self.compA.name]:
flagA = True
indA.append(ind)

# check compB
if ind in compMap[self.compB.name]:
flagB = True
indB.append(ind)

# Save the flags and indices
self.projData[ptSetName]["compA"]["flag"] = flagA
self.projData[ptSetName]["compA"]["ind"] = indA
self.projData[ptSetName]["compB"]["flag"] = flagB
self.projData[ptSetName]["compB"]["ind"] = indB

# Associate points with the tracked surfaces
for surface in self.trackSurfaces:
surfaceEps = self.trackSurfaces[surface]
if surface in self.compA.triConn:
compPoints = pts[indA]
elif surface in self.compB.triConn:
compPoints = pts[indB]
else:

# This proc has some points to project
if len(compPoints) > 0:
self.associatePointsToSurface(compPoints, ptSetName, surface, surfaceEps)

# if we include the feature curves in the warping, we also need to project the added points to the intersection and feature curves and determine how the points map to the curves
if self.incCurves:

# convert the list to an array
# we specify the dtype because numpy cannot know the type when 'indices' is empty
indices = np.array(indices, dtype="intc")

# get the coordinates of all points affected by this intersection
ptsToCurves = pts[indices]

# project these to the combined curves
# Use pySurf to project the point on curve
# Get number of points
nPoints = len(ptsToCurves)

# Initialize references if user provided none
dist2 = np.ones(nPoints, dtype=self.dtype) * 1e10
xyzProj = np.zeros((nPoints, 3), dtype=self.dtype)
tanProj = np.zeros((nPoints, 3), dtype=self.dtype)
elemIDs = np.zeros((nPoints), dtype="int32")

# Only call the Fortran code if we have at least one point
if nPoints > 0:
# This will modify xyzProj, tanProj, dist2, and elemIDs if we find better projections than dist2.
# Remember that we should adjust some indices before calling the Fortran code
# Remember to use [:] to don't lose the pointer (elemIDs is an input/output variable)
elemIDs[:] = elemIDs + 1
# (we need to do this separetely because Fortran will actively change elemIDs contents.
self.curveSearchAPI.mindistancecurve(
ptsToCurves.T, self.seam0.T, self.seamConn.T + 1, xyzProj.T, tanProj.T, dist2, elemIDs
)

# Adjust indices back to Python standards
elemIDs[:] = elemIDs - 1

# dist2 has the array of squared distances
d = np.sqrt(dist2)

# get the names of all curves including the intersection
allCurves = ["intersection"]
for curveName in self.featureCurveNames:
allCurves.append(curveName)

# track the points that dont get associated with any curve
# get a full masking array with zeros
allNodesBool = np.zeros(len(elemIDs))

# dict to save the pt indices
self.curveProjIdx[ptSetName] = {}
# dict to save other data
self.curveProjData[ptSetName] = {}

# now loop over feature curves and use the epsilon that each curve has
# to determine which points maps to which curve
for curveName in allCurves:

# get the epsilon for this curve
# we map the points closer than eps to this curve
eps = self.curveEpsDict[curveName]

# also get the range of element IDs this curve owns
seamBeg = self.seamBeg[curveName]
seamEnd = self.seamEnd[curveName]

# this returns a bool array of indices that satisfy the conditions
# we check for elemIDs because we projected to all curves at once
curveBool = np.all([d < eps, elemIDs >= seamBeg, elemIDs < seamEnd], axis=0)

# get the indices of the points mapped to this element
idxs = np.nonzero(curveBool)

# save the indices. idx has the indices of the "indices" array
# that had the indices of points that get any intersection treatment
self.curveProjIdx[ptSetName][curveName] = np.array(indices[idxs])

if self.debug:
ptCoords = ptsToCurves[idxs]
tecplot_interface.write_tecplot_scatter(
f"{curveName}.plt", curveName, ["X", "Y", "Z"], ptCoords
)

# also update the masking array
# we will use this to figure out the indices that did not get attached to any curves
allNodesBool = np.any([curveBool, allNodesBool], axis=0)

# create an empty dict to save aux data later on
self.curveProjData[ptSetName][curveName] = {}

# negate the surface mask and get indices
surfPtIdx = np.nonzero(np.logical_not(allNodesBool))

# figure out which of these surfNodes live only on components A and B
allSurfIdx = np.array(indices[surfPtIdx[0]])

# component A
# this is the local indices of the points affected

# component B

# initialize the bool dict for this pointset
self.curveProjFlag[ptSetName] = {}
self.curvePtCounts[ptSetName] = {}
self.curvePtCoords[ptSetName] = {}
self.nCurvePts[ptSetName] = {}

# we need to figure out if we have any points mapped to curves on comp A
for curveName in allCurves:

# get the indices mapped to this curve, on this proc
idxs = self.curveProjIdx[ptSetName][curveName]

# call the utility function
nPtsTotal, nPtsProcs, curvePtCoords = self._commCurveProj(pts, idxs, comm)

# save the displacements and points
self.curvePtCounts[ptSetName][curveName] = nPtsProcs
self.curvePtCoords[ptSetName][curveName] = curvePtCoords

# also save the total number for convenience
self.nCurvePts[ptSetName][curveName] = nPtsTotal

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

# original coordinates of the added pointset
pts = self.points[ptSetName][0]
# indices of the points that get affected by this intersection
indices = self.points[ptSetName][1]
# factors for each node in pointSet
factors = self.points[ptSetName][2]

# coordinates for the remeshed curves
# we use the initial seam coordinates here
coor = self.seam0
# bar connectivity for the remeshed elements
conn = self.seamConn
# deltas for each point (nNode, 3) in size
if self.seam.shape == self.seam0.shape:
dr = self.seam - self.seam0
else:
# The topology has changed so we do not update the intersection
# This will most likely break the mesh but allows
# 1) the mesh to be output for visualization
# 2) the optimization to continue after raising a fail flag
if self.comm.rank == 0:
print("The intersection topology has changed. The intersection will not be updated.")
return delta

# define an epsilon to avoid dividing by zero later on
eps = 1e-50

# loop over the points that get affected
for i in range(len(factors)):
# j is the index of the point in the full set we are working with.
j = indices[i]

# coordinates of the original point
rp = pts[j]

# Run vectorized weighted interpolation

# get the two end points for the line elements
r0 = coor[conn[:, 0]]
r1 = coor[conn[:, 1]]

# get the deltas for two end points
dr0 = dr[conn[:, 0]]
dr1 = dr[conn[:, 1]]

# compute a, b, and c coefficients
a = (r1[:, 0] - r0[:, 0]) ** 2 + (r1[:, 1] - r0[:, 1]) ** 2 + (r1[:, 2] - r0[:, 2]) ** 2
b = 2 * (
(r1[:, 0] - r0[:, 0]) * (r0[:, 0] - rp[0])
+ (r1[:, 1] - r0[:, 1]) * (r0[:, 1] - rp[1])
+ (r1[:, 2] - r0[:, 2]) * (r0[:, 2] - rp[2])
)
c = (r0[:, 0] - rp[0]) ** 2 + (r0[:, 1] - rp[1]) ** 2 + (r0[:, 2] - rp[2]) ** 2

# distances for each element
dists = np.sqrt(np.maximum(a, 0.0))

# compute some re-occurring terms
# the determinant can be zero or negative, but it CANNOT be positive
# this is because the quadratic that defines the distance from the line cannot have two roots.
# if the point is on the line, the quadratic will have a single root...
det = b * b - 4 * a * c
# these might be negative 1e-20sth so clip them...
# these will be strictly zero or greater than zero.
# Numerically, they cannot be negative bec. we are working with real numbers
sabc = np.sqrt(np.maximum(a + b + c, 0.0))
sc = np.sqrt(np.maximum(c, 0.0))

# denominators on the integral evaluations
# add an epsilon so that these terms never become zero
# det <= 0, sabc and sc >= 0, therefore the den1 and den2 should be <=0
den1 = det * sabc - eps
den2 = det * sc - eps

# integral evaluations
eval1 = (-2 * (2 * a + b) / den1 + 2 * b / den2) * dists
eval2 = ((2 * b + 4 * c) / den1 - 4 * c / den2) * dists

# denominator only gets one integral
den = np.sum(eval1)

# do each direction separately
interp = np.zeros(3, dtype=self.dtype)
for iDim in range(3):
# numerator gets two integrals with the delta components
num = np.sum((dr1[:, iDim] - dr0[:, iDim]) * eval2 + dr0[:, iDim] * eval1)
# final result
interp[iDim] = num / 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 update_d(self, ptSetName, dPt, dSeam):

"""forward mode differentiated version of the update routine.
Note that dPt and dSeam are both one dimensional arrays
"""
pts = self.points[ptSetName][0]
indices = self.points[ptSetName][1]
factors = self.points[ptSetName][2]

# we need to reshape the arrays for simpler code
dSeam = dSeam.reshape((len(self.seam0), 3))

for i in range(len(factors)):
# j is the index of the point in the full set we are
# working with.
j = indices[i]

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

# Now the delta is replaced by 1-factor times the weighted
# interp of the seam * factor of the original:
dPt[j * 3 + iDim] = factors[i] * dPt[j * 3 + iDim] + (1 - factors[i]) * interp_d[iDim]

return

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

# original coordinates of the added pointset
pts = self.points[ptSetName][0]
# indices of the points that get affected by this intersection
indices = self.points[ptSetName][1]
# factors for each node in pointSet
factors = self.points[ptSetName][2]

# coordinates for the remeshed curves
# we use the initial seam coordinates here
coor = self.seam0
# bar connectivity for the remeshed elements
conn = self.seamConn

# define an epsilon to avoid dividing by zero later on
eps = 1e-50

# if we are handling more than one function,
# seamBar will contain the seeds for each function separately
seamBar = np.zeros((dIdPt.shape[0], self.seam0.shape[0], self.seam0.shape[1]))

# if we have the projection flag, then we need to add the contribution to seamBar from that
if self.projectFlag:
seamBar += self.seamBarProj[ptSetName]

for k in range(dIdPt.shape[0]):
for i in range(len(factors)):

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

# coordinates of the original point
rp = pts[j]

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

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

# get the two end points for the line elements
r0 = coor[conn[:, 0]]
r1 = coor[conn[:, 1]]
# compute a, b, and c coefficients
a = (r1[:, 0] - r0[:, 0]) ** 2 + (r1[:, 1] - r0[:, 1]) ** 2 + (r1[:, 2] - r0[:, 2]) ** 2
b = 2 * (
(r1[:, 0] - r0[:, 0]) * (r0[:, 0] - rp[0])
+ (r1[:, 1] - r0[:, 1]) * (r0[:, 1] - rp[1])
+ (r1[:, 2] - r0[:, 2]) * (r0[:, 2] - rp[2])
)
c = (r0[:, 0] - rp[0]) ** 2 + (r0[:, 1] - rp[1]) ** 2 + (r0[:, 2] - rp[2]) ** 2
# distances for each element
dists = np.sqrt(np.maximum(a, 0.0))

# compute some re-occurring terms
det = b * b - 4 * a * c
sabc = np.sqrt(np.maximum(a + b + c, 0.0))
sc = np.sqrt(np.maximum(c, 0.0))
# denominators on the integral evaluations
den1 = det * sabc - eps
den2 = det * sc - eps
# integral evaluations
eval1 = (-2 * (2 * a + b) / den1 + 2 * b / den2) * dists
eval2 = ((2 * b + 4 * c) / den1 - 4 * c / den2) * dists

# denominator only gets one integral
den = np.sum(eval1)

evalDiff = eval1 - eval2

# do each direction separately
for iDim in range(3):
# seeds for the r0 point
seamBar[k, conn[:, 0], iDim] += localVal[iDim] * evalDiff / den

# seeds for the r1 point
seamBar[k, conn[:, 1], iDim] += localVal[iDim] * eval2 / den

# seamBar is the bwd seeds for the intersection curve...
# it is N,nseampt,3 in size
# now call the reverse differentiated seam computation
compSens = self._getIntersectionSeam_b(seamBar, comm)

return compSens

def project(self, ptSetName, newPts):
# we need to build ADTs for both components if we have any components that lie on either
# we also need to save ALL intermediate variables for gradient computations in reverse mode

# get the comm for this point set
comm = self.points[ptSetName][3]

self.comm.Barrier()

# check if we need to worry about either surface
# we will use these flags to figure out if we need to do warping.
# we need to do the comm for the updated curves regardless
flagA = False
flagB = False
if len(self.curvesOnA) > 0:
flagA = True
if len(self.curvesOnB) > 0:
flagB = True

# do the pts on the intersection outside the loop
nptsg = self.nCurvePts[ptSetName]["intersection"]

# the deltas for these points are zero. they should already be on the intersection
# also get the initial coordinates of these points. We use this during warping
# intersection curve will be on both components
if flagA:
deltaA = np.zeros((nptsg, 3), dtype=self.dtype)
curvePtCoordsA = self.curvePtCoords[ptSetName]["intersection"].copy()
else:
deltaA = np.zeros((0, 3), dtype=self.dtype)
curvePtCoordsA = np.zeros((0, 3), dtype=self.dtype)

if flagB:
deltaB = np.zeros((nptsg, 3), dtype=self.dtype)
curvePtCoordsB = self.curvePtCoords[ptSetName]["intersection"].copy()
else:
deltaB = np.zeros((0, 3), dtype=self.dtype)
curvePtCoordsB = np.zeros((0, 3), dtype=self.dtype)

# loop over the feature curves that we need to project
for curveName in self.featureCurveNames:

# get the indices of points we need to project
idx = self.curveProjIdx[ptSetName][curveName]

# these are the updated coordinates that will be projected to the curve
ptsOnCurve = newPts[idx, :].copy()

if self.debug:
tecplot_interface.write_tecplot_scatter(
f"{curveName}_warped_pts.plt", "intersection", ["X", "Y", "Z"], ptsOnCurve
)

# conn of the current curve
seamBeg = self.seamBeg[curveName]
seamEnd = self.seamEnd[curveName]
curveConn = self.seamConn[seamBeg:seamEnd]

# Project these to the combined curves using pySurf
# Get number of points
nPoints = ptsOnCurve.shape[0]

if self.debug:
print(f"[{self.comm.rank}] curveName: {curveName}, nPoints on the fwd pass: {nPoints}")

# Initialize references if user provided none
dist2 = np.ones(nPoints, dtype=self.dtype) * 1e10
xyzProj = np.zeros((nPoints, 3), dtype=self.dtype)
tanProj = np.zeros((nPoints, 3), dtype=self.dtype)
elemIDs = np.zeros((nPoints), dtype="int32")

# only call the Fortran code if we have at least one point
if nPoints > 0:
# This will modify xyzProj, tanProj, dist2, and elemIDs if we find better projections than dist2.
# Remember that we should adjust some indices before calling the Fortran code
# Remember to use [:] to don't lose the pointer (elemIDs is an input/output variable)
elemIDs[:] = (
elemIDs + 1
)  # (we need to do this separetely because Fortran will actively change elemIDs contents.
ptsOnCurve.T, self.seam.T, curveConn.T + 1, xyzProj.T, tanProj.T, dist2, elemIDs
)

# Adjust indices back to Python standards
elemIDs[:] = elemIDs - 1

# we only have the curvemask if we do the projection on this proc

# save some information for gradient comp
self.curveProjData[ptSetName][curveName]["xyz"] = ptsOnCurve.copy()
self.curveProjData[ptSetName][curveName]["coor"] = self.seam.copy()
self.curveProjData[ptSetName][curveName]["barsConn"] = curveConn.copy()
self.curveProjData[ptSetName][curveName]["xyzProj"] = xyzProj.copy()
self.curveProjData[ptSetName][curveName]["tanProj"] = tanProj.copy()
self.curveProjData[ptSetName][curveName]["elemIDs"] = elemIDs.copy()

# get the delta for the points on this proc
deltaLocal = xyzProj - ptsOnCurve

if self.debug:
tecplot_interface.write_tecplot_scatter(
f"{curveName}_projected_pts.plt", curveName, ["X", "Y", "Z"], xyzProj
)

# update the point coordinates on this processor.
# we do not need to do any communication for this
# since newPts is the final coordinates of the points we just projected
newPts[idx] = xyzProj

# communicate the deltas
if comm:
sizes = self.curvePtCounts[ptSetName][curveName]
disp = np.array([np.sum(sizes[:i]) for i in range(comm.size)], dtype="intc")

# save these for grad comp
self.curveProjData[ptSetName][curveName]["sizes"] = sizes
self.curveProjData[ptSetName][curveName]["disp"] = disp

# sendbuf
deltaLocal = deltaLocal.flatten()
sendbuf = [deltaLocal, sizes[comm.rank] * 3]

# recvbuf
nptsg = self.nCurvePts[ptSetName][curveName]
deltaGlobal = np.zeros(nptsg * 3, dtype=self.dtype)

recvbuf = [deltaGlobal, sizes * 3, disp * 3, self.mpi_type]

# do an allgatherv
comm.Allgatherv(sendbuf, recvbuf)

# reshape into a nptsg,3 array
deltaGlobal = deltaGlobal.reshape((nptsg, 3))

else:
# we dont have a comm, so this is a "serial" pointset
deltaGlobal = deltaLocal

# also save the sizes and disp stuff as if we have one proc
self.curveProjData[ptSetName][curveName]["sizes"] = self.curvePtCounts[ptSetName][curveName]
self.curveProjData[ptSetName][curveName]["disp"] = [0]

# we only add the deltaLocal to deltaA if this curve is on compA,
# and we have points on compA surface
if curveName in self.curvesOnA and flagA:
# stack the deltas
deltaA = np.vstack((deltaA, deltaGlobal))

# also stack the original coordinates for warping
curvePtCoordsNew = self.curvePtCoords[ptSetName][curveName]
curvePtCoordsA = np.vstack((curvePtCoordsA, curvePtCoordsNew))

# do the same for compB
# we can use elif bec. one curve cannot be on both comps
elif curveName in self.curvesOnB and flagB:
# stack the deltas
deltaB = np.vstack((deltaB, deltaGlobal))

# also stack the original coordinates for warping
curvePtCoordsNew = self.curvePtCoords[ptSetName][curveName]
curvePtCoordsB = np.vstack((curvePtCoordsB, curvePtCoordsNew))

self.comm.Barrier()

# then, we warp all of the nodes that were affected by the intersection treatment
# using the deltas from the previous project to curve step

if flagA:
self._warpSurfPts(self.points[ptSetName][0], newPts, self.surfIdxA[ptSetName], curvePtCoordsA, deltaA)

if flagB:
self._warpSurfPts(self.points[ptSetName][0], newPts, self.surfIdxB[ptSetName], curvePtCoordsB, deltaB)

# save some info for the sens. computations
self.curveProjData[ptSetName]["curvePtCoordsA"] = curvePtCoordsA
self.curveProjData[ptSetName]["curvePtCoordsB"] = curvePtCoordsB

# get the flags for components
flagA = self.projData[ptSetName]["compA"]["flag"]
flagB = self.projData[ptSetName]["compB"]["flag"]

# Initialize component-wide projection indices
indAComp = self.projData[ptSetName]["compA"]["ind"].copy()
indBComp = self.projData[ptSetName]["compB"]["ind"].copy()

# call the actual driver with the info to prevent code multiplication
if flagA:
# First project points on the tracked surfaces
surfaceIndA = self.projData[ptSetName]["compA"]["surfaceInd"]
for surface in surfaceIndA:
surfaceInd = surfaceIndA[surface]

# Get the subset of indices that is associated with this surface
indA = [self.projData[ptSetName]["compA"]["ind"][i] for i in surfaceInd]

# Remove these points from the component-wide projection indices
for ind in indA:
indAComp.remove(ind)

# get the points using the mapping
ptsA = newPts[indA]
# call the projection routine with the info
# this returns the projected points and we use the same mapping to put them back in place
newPts[indA] = self._projectToComponent(
ptsA, self.compA, self.projData[ptSetName][surface], surface=surface
)

# Project remaining points to the component as a whole
if indAComp:
ptsA = newPts[indAComp]
newPts[indAComp] = self._projectToComponent(ptsA, self.compA, self.projData[ptSetName]["compA"])

# do the same for B
if flagB:
surfaceIndB = self.projData[ptSetName]["compB"]["surfaceInd"]
for surface in surfaceIndB:
surfaceInd = surfaceIndB[surface]
indB = [self.projData[ptSetName]["compB"]["ind"][i] for i in surfaceInd]
for ind in indB:
indBComp.remove(ind)
ptsB = newPts[indB]
newPts[indB] = self._projectToComponent(
ptsB, self.compB, self.projData[ptSetName][surface], surface=surface
)

if indBComp:
ptsB = newPts[indBComp]
newPts[indBComp] = self._projectToComponent(ptsB, self.compB, self.projData[ptSetName]["compB"])

# Store component-wide indices for derivative computation
self.projData[ptSetName]["compA"]["indAComp"] = indAComp
self.projData[ptSetName]["compB"]["indBComp"] = indBComp

def project_b(self, ptSetName, dIdpt, comm):
# call the functions to propagate ad seeds bwd
# we need to build ADTs for both components if we have any components that lie on either
# we also need to save ALL intermediate variables for gradient computations in reverse mode

# number of functions we have
N = dIdpt.shape[0]

# get the flags for components
flagA = self.projData[ptSetName]["compA"]["flag"]
flagB = self.projData[ptSetName]["compB"]["flag"]

# Initialize dictionaries to accumulate triangulated mesh sensitivities
compSens_local = {}
compSensA = {}
compSensB = {}

# call the actual driver with the info to prevent code multiplication
if flagA:

# Project remaining points to the component as a whole
indAComp = self.projData[ptSetName]["compA"]["indAComp"]
if indAComp:
dIdptA = dIdpt[:, indAComp]
dIdpt[:, indAComp], compSensA = self._projectToComponent_b(
dIdptA, self.compA, self.projData[ptSetName]["compA"]
)

# First project points on the tracked surfaces
surfaceIndA = self.projData[ptSetName]["compA"]["surfaceInd"]
for surface in surfaceIndA:
surfaceInd = surfaceIndA[surface]

# Get the subset of indices that is associated with this surface
indA = [self.projData[ptSetName]["compA"]["ind"][i] for i in surfaceInd]

# get the points using the mapping
dIdptA = dIdpt[:, indA]
# call the projection routine with the info
# this returns the projected points and we use the same mapping to put them back in place
dIdpt[:, indA], compSensA_temp = self._projectToComponent_b(
dIdptA, self.compA, self.projData[ptSetName][surface], surface=surface
)

# Accumulate triangulated mesh sensitivities
for k, v in compSensA_temp.items():
try:
compSensA[k] += v
except KeyError:
compSensA[k] = v

for k, v in compSensA.items():
compSens_local[k] = v

# set the compSens entries to all zeros on these procs
else:
# get the values from each DVGeo
xA = self.compA.DVGeo.getValues()

# loop over each entry in xA and xB and create a dummy zero gradient array for all
for k, v in xA.items():
# create the zero array:
zeroSens = np.zeros((N, v.shape[0]))
compSens_local[k] = zeroSens

# do the same for B
if flagB:

indBComp = self.projData[ptSetName]["compB"]["indBComp"]
if indBComp:
dIdptB = dIdpt[:, indBComp]
dIdpt[:, indBComp], compSensB = self._projectToComponent_b(
dIdptB, self.compB, self.projData[ptSetName]["compB"]
)

surfaceIndB = self.projData[ptSetName]["compB"]["surfaceInd"]
for surface in surfaceIndB:
surfaceInd = surfaceIndB[surface]
indB = [self.projData[ptSetName]["compB"]["ind"][i] for i in surfaceInd]
dIdptB = dIdpt[:, indB]
dIdpt[:, indB], compSensB_temp = self._projectToComponent_b(
dIdptB, self.compB, self.projData[ptSetName][surface], surface=surface
)

for k, v in compSensB_temp.items():
try:
compSensB[k] += v
except KeyError:
compSensB[k] = v

for k, v in compSensB.items():
compSens_local[k] = v

# set the compSens entries to all zeros on these procs
else:
# get the values from each DVGeo
xB = self.compB.DVGeo.getValues()

# loop over each entry in xA and xB and create a dummy zero gradient array for all
for k, v in xB.items():
# create the zero array:
zeroSens = np.zeros((N, v.shape[0]))
compSens_local[k] = zeroSens

# finally sum the results across procs if we are provided with a comm
if comm:
compSens = {}
# because the results are in a dictionary, we need to loop over the items and sum
for k in compSens_local:
compSens[k] = comm.allreduce(compSens_local[k], op=MPI.SUM)
else:
# we can just pass the dictionary
compSens = compSens_local

# now we do the warping

# get the comm for this point set
ptSetComm = self.points[ptSetName][3]

# check if any processor did any warping on compA
curvePtCoordsA = self.curveProjData[ptSetName]["curvePtCoordsA"]
curvePtCoordsB = self.curveProjData[ptSetName]["curvePtCoordsB"]
if ptSetComm:
nCurvePtCoordsAG = ptSetComm.allreduce(len(curvePtCoordsA), op=MPI.MAX)
nCurvePtCoordsBG = ptSetComm.allreduce(len(curvePtCoordsB), op=MPI.MAX)
else:
nCurvePtCoordsAG = len(curvePtCoordsA)
nCurvePtCoordsBG = len(curvePtCoordsB)

# check if we need to worry about either surface
# we will use these flags to figure out if we need to do warping.
# we need to do the comm for the updated curves regardless
flagA = False
flagB = False
if len(self.curvesOnA) > 0:
flagA = True
if len(self.curvesOnB) > 0:
flagB = True

if ptSetComm:
rank = ptSetComm.rank
else:
rank = 0

# call the bwd warping routine
# deltaA_b is the seed for the points projected to curves
if flagA:
deltaA_b_local = self._warpSurfPts_b(
dIdpt, self.points[ptSetName][0], self.surfIdxA[ptSetName], curvePtCoordsA
)
else:
deltaA_b_local = np.zeros((N, nCurvePtCoordsAG, 3))

# do the same for comp B
if flagB:
deltaB_b_local = self._warpSurfPts_b(
dIdpt, self.points[ptSetName][0], self.surfIdxB[ptSetName], curvePtCoordsB
)
else:
deltaB_b_local = np.zeros((N, nCurvePtCoordsBG, 3))

# reduce seeds for both
if ptSetComm:
deltaA_b = ptSetComm.allreduce(deltaA_b_local, op=MPI.SUM)
deltaB_b = ptSetComm.allreduce(deltaB_b_local, op=MPI.SUM)
# no comm, local is global
else:
deltaA_b = deltaA_b_local
deltaB_b = deltaB_b_local

# remove the seeds for the intersection, the disps are zero (constant) so no need to diff. those
if flagA:
# this just returns the seeds w/o the intersection seeds
deltaA_b = deltaA_b[:, self.nCurvePts[ptSetName]["intersection"] :]
if flagB:
# this just returns the seeds w/o the intersection seeds
deltaB_b = deltaB_b[:, self.nCurvePts[ptSetName]["intersection"] :]

# loop over the curves
for curveName in self.featureCurveNames:

# sizes and displacements for this curve
sizes = self.curveProjData[ptSetName][curveName]["sizes"]
disp = self.curveProjData[ptSetName][curveName]["disp"]

# we get the seeds from compA seeds
if curveName in self.curvesOnA:

if flagA:
# contribution on this proc
deltaBar = deltaA_b[:, disp[rank] : disp[rank] + sizes[rank]].copy()

# this proc does not have any pts projected, so set the seed to zero
else:
deltaBar = np.zeros((N, 0, 3))

# remove the seeds for this curve from deltaA_b seeds
deltaA_b = deltaA_b[:, disp[-1] + sizes[-1] :]

# seeds from compB
elif curveName in self.curvesOnB:

if flagB:
# contribution on this proc
deltaBar = deltaB_b[:, disp[rank] : disp[rank] + sizes[rank]].copy()

# this proc does not have any pts projected, so set the seed to zero
else:
deltaBar = np.zeros((N, 0, 3))

# remove the seeds for this curve from deltaA_b seeds
deltaB_b = deltaB_b[:, disp[-1] + sizes[-1] :]

else:
print("This should not happen")

# now we have the local seeds of projection points for all functions in xyzProjb

# get the indices of points we need to project
idx = self.curveProjIdx[ptSetName][curveName]
nPoints = len(idx)

# get some data from the fwd run
xyz = self.curveProjData[ptSetName][curveName]["xyz"]
coor = self.curveProjData[ptSetName][curveName]["coor"]
barsConn = self.curveProjData[ptSetName][curveName]["barsConn"]
xyzProj = self.curveProjData[ptSetName][curveName]["xyzProj"]
tanProj = self.curveProjData[ptSetName][curveName]["tanProj"]
elemIDs = self.curveProjData[ptSetName][curveName]["elemIDs"]

# we dont use tangents so seeds are zero
tanProjb = np.zeros_like(tanProj)

if nPoints > 0:

# run the bwd projection for everyfunction
for k in range(N):

# contribution from delta
xyzProjb = deltaBar[k].copy()

# add the contribution from dIdpt for the idx points themselves
xyzProjb += dIdpt[k, idx]

# Call Fortran code (This will accumulate seeds in xyzb and self.coorb)
xyzb_new, coorb_new = self.curveSearchAPI.mindistancecurve_b(
xyz.T,
coor.T,
barsConn.T + 1,
xyzProj.T,
xyzProjb.T,
tanProj.T,
tanProjb.T,
elemIDs + 1,
)

# Accumulate derivatives with the correct k

# replace the seed in dIdpt. we subtract deltaBar here
# because delta was equal to the projected points minus the original points.
# So the delta seeds contribute with a negative sign
# to the seeds of the points before projection
dIdpt[k, idx, :] = xyzb_new.T - deltaBar[k]

# add the seed to the seam seed
self.seamBarProj[ptSetName][k, :, :] += coorb_new.T

return compSens

def _commCurveProj(self, pts, indices, comm):
"""
This function will get the points, indices, and comm.
This function is called once for each feature curve.
The indices are the indices of points that was mapped to this curve.
We compute how many points we have mapped to this curve globally.
Furthermore, we compute the displacements.
Finally, we communicate the initial coordinates of these points.
These will later be used in the point-based warping.

"""

# only do this fancy stuff if this is a "parallel" pointset
if comm:
nproc = comm.size

# communicate the counts
sizes = np.array(comm.allgather(len(indices)), dtype="intc")

# total number of points
nptsg = np.sum(sizes)

# get the displacements
disp = np.array([np.sum(sizes[:i]) for i in range(nproc)], dtype="intc")

# sendbuf
ptsLocal = pts[indices].flatten()
sendbuf = [ptsLocal, len(indices) * 3]

# recvbuf
ptsGlobal = np.zeros(3 * nptsg, dtype=self.dtype)

recvbuf = [ptsGlobal, sizes * 3, disp * 3, self.mpi_type]

# do an allgatherv
comm.Allgatherv(sendbuf, recvbuf)

# reshape into a nptsg,3 array
curvePtCoords = ptsGlobal.reshape((nptsg, 3))

# this is a "serial" pointset, so the results are just local
else:
nptsg = len(indices)
sizes = [nptsg]
curvePtCoords = pts[indices]

return nptsg, sizes, curvePtCoords

def _warpSurfPts(self, pts0, ptsNew, indices, curvePtCoords, delta):
"""
This function warps points using the displacements from curve projections.

pts0: The original surface point coordinates.
ptsNew: Updated surface pt coordinates. We will add the warped delta to these inplace.
indices: Indices of the points that we will use for this operation.
curvePtCoords: Original coordinates of points on curves.
delta: Displacements of the points on curves after projecting them.

"""

# Return if curvePtCoords is empty
if not np.any(curvePtCoords):
return

for j in indices:

# point coordinates with the baseline design
# this is the point we will warp
ptCoords = pts0[j]

# Vectorized point-based warping
rr = ptCoords - curvePtCoords
LdefoDist = 1.0 / np.sqrt(rr[:, 0] ** 2 + rr[:, 1] ** 2 + rr[:, 2] ** 2 + 1e-16)
LdefoDist3 = LdefoDist**3
Wi = LdefoDist3
den = np.sum(Wi)
interp = np.zeros(3, dtype=self.dtype)
for iDim in range(3):
interp[iDim] = np.sum(Wi * delta[:, iDim]) / den

# finally, update the coord in place
ptsNew[j] = ptsNew[j] + interp

def _warpSurfPts_b(self, dIdPt, pts0, indices, curvePtCoords):

# seeds for delta
deltaBar = np.zeros((dIdPt.shape[0], curvePtCoords.shape[0], 3))

# Return zeros if curvePtCoords is empty
if not np.any(curvePtCoords):
return deltaBar

for k in range(dIdPt.shape[0]):

for j in indices:

# point coordinates with the baseline design
# this is the point we will warp
ptCoords = pts0[j]

# local seed for 3 coords
localVal = dIdPt[k, j]

# Vectorized point-based warping
rr = ptCoords - curvePtCoords
LdefoDist = 1.0 / np.sqrt(rr[:, 0] ** 2 + rr[:, 1] ** 2 + rr[:, 2] ** 2 + 1e-16)
LdefoDist3 = LdefoDist**3
Wi = LdefoDist3
den = np.sum(Wi)

for iDim in range(3):
deltaBar[k, :, iDim] += Wi * localVal[iDim] / den

# return the seeds for the delta vector
return deltaBar

def _projectToComponent(self, pts, comp, projDict, surface=None):
# We build an ADT for this component using pySurf
# Set bounding box for new tree
BBox = np.zeros((2, 3))
useBBox = False

# dummy connectivity data for quad elements since we have all tris

if surface is not None:
# Use the triConn for just this surface
triConn = comp.triConn[surface]
# Set the adtID as the surface name
else:
# Use the stacked triConn for the whole component
triConn = comp.triConnStack
# Set the adtID as the component name

# Compute set of nodal normals by taking the average normal of all
# elements surrounding the node. This allows the meshing algorithms,
# for instance, to march in an average direction near kinks.
comp.nodal_normals = nodal_normals.T

# Create new tree (the tree itself is stored in Fortran level)
)

# project
numPts = pts.shape[0]
dist2 = np.ones(numPts, dtype=self.dtype) * 1e10
xyzProj = np.zeros((numPts, 3), dtype=self.dtype)
normProjNotNorm = np.zeros((numPts, 3), dtype=self.dtype)

if self.debug:
print(f"[{self.comm.rank}] Projecting to component {comp.name}, pts.shape = {pts.shape}")

# Call projection function
pts.T, adtID, dist2, xyzProj.T, comp.nodal_normals.T, normProjNotNorm.T
)

elementID = elementID - 1
uvw = uvw.T

# normalize the normals
normProj = tsurf_tools.normalize(normProjNotNorm)

# save the data
projDict["procID"] = procID.copy()
projDict["elementType"] = elementType.copy()
projDict["elementID"] = elementID.copy()
projDict["uvw"] = uvw.copy()
projDict["dist2"] = dist2.copy()
projDict["normProjNotNorm"] = normProjNotNorm.copy()
projDict["normProj"] = normProj.copy()

# also save the original and projected points
projDict["xyz"] = pts.copy()
projDict["xyzProj"] = xyzProj.copy()

# return projected points
return xyzProj

def _projectToComponent_b(self, dIdpt, comp, projDict, surface=None):

# We build an ADT for this component using pySurf
# Set bounding box for new tree
BBox = np.zeros((2, 3))
useBBox = False

# dummy connectivity data for quad elements since we have all tris

if surface is not None:
# Use the triConn for just this surface
triConn = comp.triConn[surface]
# Set the adtID as the surface name
else:
# Use the stacked triConn for the whole component
triConn = comp.triConnStack
# Set the adtID as the component name

# Compute set of nodal normals by taking the average normal of all
# elements surrounding the node. This allows the meshing algorithms,
# for instance, to march in an average direction near kinks.
comp.nodal_normals = nodal_normals.T

# Create new tree (the tree itself is stored in Fortran level)
)

# also extract the projection data we have from the fwd pass
procID = projDict["procID"]
elementType = projDict["elementType"]
elementID = projDict["elementID"]
uvw = projDict["uvw"]
dist2 = projDict["dist2"]
normProjNotNorm = projDict["normProjNotNorm"]

# get the original and projected points too
xyz = projDict["xyz"]
xyzProj = projDict["xyzProj"]

# also, we dont care about the normals, so the AD seeds for them should (?) be zero
normProjb = np.zeros_like(normProjNotNorm)

# also create the dIdtp for the triangulated surface nodes
dIdptComp = np.zeros((dIdpt.shape[0], comp.nodes.shape[0], 3))

# now propagate the ad seeds back for each function
for i in range(dIdpt.shape[0]):
# the derivative seeds for the projected points
xyzProjb = dIdpt[i].copy()

# Compute derivatives of the normalization process
normProjNotNormb = tsurf_tools.normalize_b(normProjNotNorm, normProjb)

# Call projection function
# ATTENTION: The variable "xyz" here in Python corresponds to the variable "coor" in the Fortran code.
# On the other hand, the variable "coor" here in Python corresponds to the variable "adtCoor" in Fortran.
# I could not change this because the original ADT code already used "coor" to denote nodes that should be
# projected.

xyz.T,
procID,
elementType,
elementID + 1,
uvw.T,
dist2,
xyzProj.T,
xyzProjb.T,
comp.nodal_normals.T,
normProjNotNorm.T,
normProjNotNormb.T,
)

# Transpose results to make them consistent
xyzb = xyzb.T
coorb = coorb.T

# Put the reverse ad seed back into dIdpt
dIdpt[i] = xyzb
# Also save the triangulated surface node seeds
dIdptComp[i] = coorb

# Now we are done with the ADT

# Call the total sensitivity of the component's DVGeo
compSens = comp.DVGeo.totalSensitivity(dIdptComp, "triMesh")

# the entries in dIdpt is replaced with AD seeds of initial points that were projected
# we also return the total sensitivity contributions from components' triMeshes
return dIdpt, compSens

def _getUpdatedCoords(self):
# this code returns the updated coordinates

# first comp a
self.compA.updateTriMesh()

# then comp b
self.compB.updateTriMesh()

return

def _getIntersectionSeam(self, comm, firstCall=False):
# we can parallelize here. each proc gets one intersection, but needs re-structuring of some of the code.

# this function computes the intersection curve, cleans up the data and splits the curve based on features or curves specified by the user.

# create the dictionary to save all intermediate variables for reverse differentiation
self.seamDict = {}

# Call pySurf with the quad information.
dummyConn = np.zeros((0, 4))

# compute the intersection curve, in the first step we just get the array sizes to hide allocatable arrays from python
arraySizes = self.intersectionAPI.computeintersection(
self.compA.nodes.T,
self.compA.triConnStack.T,
dummyConn.T,
self.compB.nodes.T,
self.compB.triConnStack.T,
dummyConn.T,
self.distTol,
comm.py2f(),
)

# Retrieve results from Fortran if we have an intersection
if np.max(arraySizes[1:]) > 0:

# Second Fortran call to retrieve data from the CGNS file.
intersectionArrays = self.intersectionAPI.retrievedata(*arraySizes)

# We need to do actual copies, otherwise data will be overwritten if we compute another intersection.
# We subtract one to make indices consistent with the Python 0-based indices.
# We also need to transpose it since Python and Fortran use different orderings to store matrices in memory.

intNodes = np.array(intersectionArrays[0]).T
# The last entry of intNodes is always 0,0,0 for some reason, even for a proper intersection
# This is probably a bug in pySurf
# It has no effect here because the computation relies on connectivity, not individual points
barsConn = np.array(intersectionArrays[1]).T - 1
parentTria = np.array(intersectionArrays[2]).T - 1

# Save these intermediate variables
self.seamDict["barsConn"] = barsConn
self.seamDict["parentTria"] = parentTria

else:
raise Error(f"The components {self.compA.name} and {self.compB.name} do not intersect.")

# Release memory used by Fortran
self.intersectionAPI.releasememory()

# Sort the output
newConn, newMap = tsurf_tools.FEsort(barsConn.tolist())

# newConn might have multiple intersection curves
if len(newConn) == 1:
# we have a single intersection curve, just take this.
seamConn = newConn[0].copy()
# We have multiple intersection curves
else:
if self.intDir is None:
# we have multiple intersection curves but the user did not specify which direction to pick
for i in range(len(newConn)):
curvename = f"{self.compA.name}_{self.compB.name}_{i}"
tecplot_interface.writeTecplotFEdata(intNodes, newConn[i], curvename, curvename)
raise Error(
f"More than one intersection curve between comps {self.compA.name} and {self.compB.name}. "
+ "The curves are written as Tecplot files in the current directory. "
+ "Try rerunning after specifying intDir for the intersection."
)

# the user did specify which direction to pick
else:
int_centers = np.zeros(len(newConn), dtype=self.dtype)
# we will figure out the locations of these points and pick the one closer to the user picked direction
for i in range(len(newConn)):
# get all the points
int_pts = intNodes[newConn[i]][:, 0]

# average the values
# the API uses a 1 based indexing, but here, we convert to a zero based indexing
int_centers[i] = np.average(int_pts[abs(self.intDir) - 1])

# multiply the values with the sign of intDir
int_centers *= np.sign(self.intDir)

# get the argmax
int_index = np.argmax(int_centers)

# this is the intersection seam
seamConn = newConn[int_index].copy()

# Get the number of elements
nElem = seamConn.shape[0]

# now that we have a continuous, ordered seam connectivity in seamConn, we can try to detect features

# we need to track the nodes that are closest to the supplied feature curves
breakList = []
curveBeg = {}
curveBegCoor = {}

# loop over the feature curves
for curveName in self.featureCurveNames:

# we need to initialize the dictionary here
# to get the intermediate output from mindistancecurve call
self.seamDict[curveName] = {}

# if this curve is on compB, we use it to track intersection features
if curveName in self.compB.barsConn and curveName not in self.remeshAll:

# get the curve connectivity
curveConn = self.compB.barsConn[curveName]

# Use pySurf to project the point on curve
# First, we need to get a list of nodes that define the intersection
intNodesOrd = intNodes[seamConn[:, 0]]

# Get number of points
nPoints = len(intNodesOrd)

# Initialize references
dist2 = np.ones(nPoints, dtype=self.dtype) * 1e10
xyzProj = np.zeros((nPoints, 3), dtype=self.dtype)
tanProj = np.zeros((nPoints, 3), dtype=self.dtype)
elemIDs = np.zeros((nPoints), dtype="int32")

# then find the closest point to the curve

# only call the Fortran code if we have at least one point
# this is redundant but it is how its done in pySurf
if nPoints > 0:
# This will modify xyzProj, tanProj, dist2, and elemIDs if we find better projections than dist2.
# Remember that we should adjust some indices before calling the Fortran code
# Remember to use [:] to don't lose the pointer (elemIDs is an input/output variable)
elemIDs[:] = (
elemIDs + 1
)  # (we need to do this separetely because Fortran will actively change elemIDs contents.
intNodesOrd.T, self.compB.nodes.T, curveConn.T + 1, xyzProj.T, tanProj.T, dist2, elemIDs
)

# Adjust indices back to Python standards
elemIDs[:] = elemIDs - 1

self.seamDict[curveName]["elemIDs"] = elemIDs.copy()
self.seamDict[curveName]["intNodesOrd"] = intNodesOrd.copy()
self.seamDict[curveName]["xyzProj"] = xyzProj.copy()
self.seamDict[curveName]["tanProj"] = tanProj.copy()
self.seamDict[curveName]["dist2"] = dist2.copy()
self.seamDict[curveName]["projPtIndx"] = seamConn[:, 0][np.argmin(dist2)].copy()

# now, find the index of the smallest distance
breakList.append(np.argmin(dist2))

# also get which element is the closest to the feature point
curveBeg[curveName] = elemIDs[np.argmin(dist2)]

# get which point on this element we projected to.
curveBegCoor[curveName] = xyzProj[np.argmin(dist2)]

else:
# if it is not on compB, we still need to set up some variables so that we remesh the whole curve
# set the beginning to the first element
curveBeg[curveName] = 0

# number of features we detected. This will be equal to the number of feature curves on compB
nFeature = len(breakList)

# if this is the first call,
if firstCall:
# we also save the initial curve with nodes and connectivities for distance calculations
self.conn0 = seamConn.copy()
self.nodes0 = intNodes.copy()
self.nFeature = nFeature
else:
if nFeature != self.nFeature:
raise Error("Number of features on the intersection curve has changed.")

# flip
# we want breakList to be in increasing order...
ii = 0
for i in range(nFeature):
# we loop over the breaklist elements and check if the element index is going up or down
if breakList[i] < breakList[np.mod(i + 1, nFeature)]:
ii += 1

# now check if we need to flip the curve
if ii == 1:  # we need at least 2 features where the element number increases...
# we need to reverse the order of our feature curves
# and we will flip the elements too so keep track of this change
breakList = np.mod(seamConn.shape[0] - np.array(breakList), seamConn.shape[0])

# and we need to invert the curves themselves
seamConn = np.flip(seamConn, axis=0)
seamConn = np.flip(seamConn, axis=1)

# roll so that the first breakList entry is the first node
seamConn = np.roll(seamConn, -breakList[0], axis=0)

# also adjust the feature indices
breakList = np.mod(breakList - breakList[0], nElem)

# get the number of elements between each feature
curveSizes = []
for i in range(nFeature - 1):
curveSizes.append(np.mod(breakList[i + 1] - breakList[i], nElem))
# check the last curve outside the loop
curveSizes.append(np.mod(breakList[0] - breakList[-1], nElem))

# copy the curveSizes for the first call
if firstCall:
self.nElems = curveSizes[:]

# now loop over the curves between the feature nodes. We will remesh them separately to retain resolution between curve features, and just append the results since the features are already ordered
curInd = 0
seam = np.zeros((0, 3), dtype=self.dtype)
finalConn = np.zeros((0, 2), dtype="int32")
for i in range(nFeature):
# just use the same number of points *2 for now
nNewNodes = self.nElems[i] + 1
coor = intNodes
barsConn = seamConn[curInd : curInd + curveSizes[i]]
curInd += curveSizes[i]
method = "linear"
spacing = "linear"
initialSpacing = 0.1
finalSpacing = 0.1

# re-sample the curve (try linear for now), to get N number of nodes on it spaced linearly
# Call Fortran code. Remember to adjust transposes and indices
newCoor, newBarsConn = self.utilitiesAPI.remesh(
nNewNodes, coor.T, barsConn.T + 1, method, spacing, initialSpacing, finalSpacing
)
newCoor = newCoor.T
newBarsConn = newBarsConn.T - 1

# add these n -resampled nodes back to back in seam and return a copy of the array
# we don't need the connectivity info for now? we just need the coords
# first increment the new connectivity by number of coordinates already in seam
newBarsConn += len(seam)

# now stack the nodes
seam = np.vstack((seam, newCoor))

# and the conn
finalConn = np.vstack((finalConn, newBarsConn))

if firstCall:
# save the beginning and end indices of these elements
self.seamBeg["intersection"] = 0
self.seamEnd["intersection"] = len(finalConn)

# save stuff to the dictionary for sensitivity computations...
self.seamDict["intNodes"] = intNodes.copy()
self.seamDict["seamConn"] = seamConn.copy()
self.seamDict["curveSizes"] = curveSizes.copy()
# size of the intersection seam w/o any feature curves
self.seamDict["seamSize"] = len(seam)
self.seamDict["curveBegCoor"] = curveBegCoor.copy()

# Output the intersection curve
if self.comm.rank == 0 and self.debug:
curvename = f"{self.compA.name}_{self.compB.name}_{self.counter}"
tecplot_interface.writeTecplotFEdata(intNodes, seamConn, curvename, curvename)

# we need to re-mesh feature curves if the user wants...
if self.incCurves:

# we need to set up some variables
if firstCall:
self.nNodeFeature = {}
self.distFeature = {}

remeshedCurves = np.zeros((0, 3), dtype=self.dtype)
remeshedCurveConn = np.zeros((0, 2), dtype="int32")

# loop over each curve, figure out what nodes get re-meshed, re-mesh, and append to seam...
for curveName in self.featureCurveNames:

# figure out which comp owns this curve...
if curveName in self.compB.barsConn:
curveComp = self.compB
dStarComp = self.dStarB
elif curveName in self.compA.barsConn:
curveComp = self.compA
dStarComp = self.dStarA

# connectivity for this curve.
curveConn = curveComp.barsConn[curveName]

# we already have the element that is closest to the intersection
# or if this curve does not start from the intersection,
# this is simply the first element
elemBeg = curveBeg[curveName]

# now lets split this element so that we get a better initial point...
# this has to be on compB
if curveName in curveBegCoor:
# save the original coordinate of the first point
ptBegSave = self.compB.nodes[curveConn[elemBeg, 0]].copy()
# and replace this with the starting point we want
self.compB.nodes[curveConn[elemBeg, 0]] = curveBegCoor[curveName].copy()

# compute the element lengths starting from elemBeg
firstNodes = curveComp.nodes[curveConn[elemBeg:, 0]]
secondNodes = curveComp.nodes[curveConn[elemBeg:, 1]]
diff = secondNodes - firstNodes
dist2 = diff[:, 0] ** 2 + diff[:, 1] ** 2 + diff[:, 2] ** 2
elemDist = np.sqrt(dist2)

# get the cumulative distance
cumDist = np.cumsum(elemDist)

# if no marchdir for this curve, we use all of it
if curveName in self.remeshAll:
# we remesh all of this curve
elemBeg = 0
elemEnd = len(curveConn)
if firstCall:
self.nNodeFeature[curveName] = elemEnd + 1
else:
# do the regular thing
if firstCall:
# compute the distances from curve nodes to intersection seam
curvePts = curveComp.nodes[curveConn[elemBeg:, 0]]

# Get number of points
nPoints = len(curvePts)

# Initialize references if user provided none
dist2 = np.ones(nPoints, dtype=self.dtype) * 1e10
xyzProj = np.zeros((nPoints, 3), dtype=self.dtype)
tanProj = np.zeros((nPoints, 3), dtype=self.dtype)
elemIDs = np.zeros((nPoints), dtype="int32")

# then find the closest point to the curve

# Only call the Fortran code if we have at least one point
if nPoints > 0:
# This will modify xyzProj, tanProj, dist2, and elemIDs if we find better projections than dist2.
# Remember that we should adjust some indices before calling the Fortran code
# Remember to use [:] to don't lose the pointer (elemIDs is an input/output variable)
elemIDs[:] = (
elemIDs + 1
)  # (we need to do this separetely because Fortran will actively change elemIDs contents.
curvePts.T, self.nodes0.T, self.conn0.T + 1, xyzProj.T, tanProj.T, dist2, elemIDs
)

dNodes = np.sqrt(dist2)

# number of elements to use, subtract one to get the correct element count
nElem = (np.abs(dNodes - dStarComp * 1.3)).argmin() - 1

# we want to be one after the actual distance, so correct if needed
if dNodes[nElem] < dStarComp * 1.3:
nElem += 1

elemEnd = elemBeg + nElem

# get the total curve distance from elemBeg to this element.
distCurve = cumDist[nElem]

# save this distance as the remesh distance
self.distFeature[curveName] = distCurve

# also save how many nodes we have, we want 2 times this when re-meshing
self.nNodeFeature[curveName] = nElem + 1

else:
# figure out how many elements we need to go in this direction
elemEnd = (np.abs(cumDist - self.distFeature[curveName])).argmin() + elemBeg

# get the new connectivity data between the initial and final elements
curveConnTrim = curveConn[elemBeg : elemEnd + 1]

# remesh the new connectivity curve, using nNode*2 times nodes
nNewNodes = self.nNodeFeature[curveName]
coor = curveComp.nodes
barsConn = curveConnTrim
method = "linear"
spacing = "linear"
initialSpacing = 0.1
finalSpacing = 0.1

# now re-sample the curve (try linear for now), to get N number of nodes on it spaced linearly
# Call Fortran code. Remember to adjust transposes and indices
newCoor, newBarsConn = self.utilitiesAPI.remesh(
nNewNodes, coor.T, barsConn.T + 1, method, spacing, initialSpacing, finalSpacing
)
newCoor = newCoor.T
newBarsConn = newBarsConn.T - 1

# increment the connectivitiy data
newBarsConn += len(remeshedCurves)

# append this new curve to the featureCurve data
remeshedCurves = np.vstack((remeshedCurves, newCoor))
remeshedCurveConn = np.vstack((remeshedCurveConn, newBarsConn))

# number of new nodes added in the opposite direction
nNewNodesReverse = 0
if elemBeg > 0 and self.remeshBwd:
# also re-mesh the initial part of the curve, to prevent any negative volumes there
curveConnTrim = curveConn[:elemBeg]

nNewNodesReverse = self.nNodeFeature[curveName]
coor = self.compB.nodes
barsConn = curveConnTrim
method = "linear"
spacing = "linear"
initialSpacing = 0.1
finalSpacing = 0.1

# now re-sample the curve (try linear for now), to get N number of nodes on it spaced linearly
# Call Fortran code. Remember to adjust transposes and indices
newCoor, newBarsConn = self.utilitiesAPI.remesh(
nNewNodesReverse, coor.T, barsConn.T + 1, method, spacing, initialSpacing, finalSpacing
)
newCoor = newCoor.T
newBarsConn = newBarsConn.T - 1

newBarsConn = newBarsConn + len(remeshedCurves)

remeshedCurves = np.vstack((remeshedCurves, newCoor))
remeshedCurveConn = np.vstack((remeshedCurveConn, newBarsConn))

if curveName in curveBegCoor:
# finally, put the modified initial and final points back in place.
self.compB.nodes[curveConn[elemBeg, 0]] = ptBegSave.copy()

# save some info for gradient computations later on
self.seamDict[curveName]["nNewNodes"] = nNewNodes
self.seamDict[curveName]["nNewNodesReverse"] = nNewNodesReverse
self.seamDict[curveName]["elemBeg"] = elemBeg
self.seamDict[curveName]["elemEnd"] = elemEnd
# this includes the initial coordinates of the points for each curve
# that has a modified initial curve
self.seamDict["curveBegCoor"] = curveBegCoor

if firstCall:
# save the beginning and end indices of these elements
self.seamBeg[curveName] = (
len(finalConn) + len(remeshedCurveConn) - (nNewNodes + nNewNodesReverse) + 2
)
self.seamEnd[curveName] = len(finalConn) + len(remeshedCurveConn)

# Output the feature curves
if self.comm.rank == 0 and self.debug:
curvename = f"featureCurves_{self.counter}"
tecplot_interface.writeTecplotFEdata(remeshedCurves, remeshedCurveConn, curvename, curvename)

# now we are done going over curves,
# so we can append all the new curves to the "seam",
# which now contains the intersection, and re-meshed feature curves

# increment the conn from curves
remeshedCurveConn += len(seam)
# stack the nodes
seam = np.vstack((seam, remeshedCurves))
# stack the conn
finalConn = np.vstack((finalConn, remeshedCurveConn))

# save the connectivity
self.seamConn = finalConn

self.counter += 1

return seam.copy()

def _getIntersectionSeam_b(self, seamBar, comm):
# seamBar contains all the bwd seeds for all coordinates in self.seam

# seam bar has shape [N, nSeamPt, 3]
# seeds for N functions
N = seamBar.shape[0]
# n points in total in the combined seam
# 3 coordinates for each point

# allocate the space for component coordinate seeds
coorAb = np.zeros((N, self.compA.nodes.shape[0], self.compA.nodes.shape[1]))
coorBb = np.zeros((N, self.compB.nodes.shape[0], self.compB.nodes.shape[1]))

# first, extract the actual intersection coordinates from feature curves.
# we might not have any feature curves but the intersection curve will be there
seamSize = self.seamDict["seamSize"]

# coordinates of the beginning points for feature curves on compB
curveBegCoor = self.seamDict["curveBegCoor"]

intBar = seamBar[:, :seamSize, :]
curveBar = seamBar[:, seamSize:, :]

# dictionary to save the accumulation of curve projection seeds
curveProjb = {}

# check if we included feature curves
if self.incCurves:

# offset for the derivative seeds for this curve
iBeg = 0

# loop over each curve
for curveName in self.featureCurveNames:
# get the fwd data
curveDict = self.seamDict[curveName]
nNewNodes = curveDict["nNewNodes"]
elemBeg = curveDict["elemBeg"]
elemEnd = curveDict["elemEnd"]

# get the derivative seeds
newCoorb = curveBar[:, iBeg : iBeg + nNewNodes, :].copy()
iBeg += nNewNodes

# figure out which comp owns this curve...
if curveName in self.compB.barsConn:
curveComp = self.compB
coorb = coorBb

elif curveName in self.compA.barsConn:
curveComp = self.compA
coorb = coorAb

# connectivity for this curve.
curveConn = curveComp.barsConn[curveName]

# adjust the first coordinate of the curve
if curveName in curveBegCoor:
# save the original coordinate of the first point
ptBegSave = self.compB.nodes[curveConn[elemBeg, 0]].copy()
# and replace this with the starting point we want
self.compB.nodes[curveConn[elemBeg, 0]] = curveBegCoor[curveName].copy()

# get the coordinates of points
coor = curveComp.nodes

# connectivity for this curve.
curveConn = curveComp.barsConn[curveName]

# trim the connectivity data
barsConn = curveComp.barsConn[curveName][elemBeg : elemEnd + 1]

# constant inputs
method = "linear"
spacing = "linear"
initialSpacing = 0.1
finalSpacing = 0.1

cb = np.zeros((N, coor.shape[0], coor.shape[1]))

# loop over functions
for ii in range(N):
# Call Fortran code. Remember to adjust transposes and indices
_, _, cbi = self.utilitiesAPI.remesh_b(
nNewNodes - 1,
coor.T,
newCoorb[ii].T,
barsConn.T + 1,
method,
spacing,
initialSpacing,
finalSpacing,
)
# derivative seeds for the coordinates.
cb[ii] = cbi.T.copy()

# check if we adjusted the initial coordinate of the curve w/ a seam coordinate
if elemBeg > 0:
if self.remeshBwd:
# first, we need to do the re-meshing of the other direction

# get the fwd data
nNewNodes = curveDict["nNewNodesReverse"]

# get the derivative seeds
newCoorb = curveBar[:, iBeg : iBeg + nNewNodes, :].copy()
iBeg += nNewNodes

# bars conn is everything up to elemBeg
barsConn = curveComp.barsConn[curveName][:elemBeg]

# loop over functions
for ii in range(N):
# Call Fortran code. Remember to adjust transposes and indices
_, _, cbi = self.utilitiesAPI.remesh_b(
nNewNodes - 1,
coor.T,
newCoorb[ii].T,
barsConn.T + 1,
method,
spacing,
initialSpacing,
finalSpacing,
)
# derivative seeds for the coordinates.
cb[ii] += cbi.T.copy()

# the first seed is for the projected point...
projb = cb[:, curveConn[elemBeg, 0], :].copy()

# zero out the seed of the replaced node
cb[:, curveConn[elemBeg, 0], :] = np.zeros((N, 3))

# put the modified initial and final points back in place.
self.compB.nodes[curveConn[elemBeg, 0]] = ptBegSave.copy()

# we need to call the curve projection routine to propagate the seed...
intNodesOrd = self.seamDict[curveName]["intNodesOrd"]
elemIDs = self.seamDict[curveName]["elemIDs"]
xyzProj = self.seamDict[curveName]["xyzProj"]
tanProj = self.seamDict[curveName]["tanProj"]
dist2 = self.seamDict[curveName]["dist2"]

# we need the full bars conn for this
barsConn = curveComp.barsConn[curveName]

# allocate zero seeds
xyzProjb = np.zeros_like(xyzProj)
tanProjb = np.zeros_like(tanProj)

curveProjb[curveName] = np.zeros((N, 3))

for ii in range(N):

# the only nonzero seed is indexed by argmin dist2
xyzProjb[np.argmin(dist2)] = projb[ii].copy()

xyzb_new, coorb_new = self.curveSearchAPI.mindistancecurve_b(
intNodesOrd.T,
self.compB.nodes.T,
barsConn.T + 1,
xyzProj.T,
xyzProjb.T,
tanProj.T,
tanProjb.T,
elemIDs + 1,
)

# add the coorb_new to coorBb[ii] since coorb_new has the seeds from mindistancecurve_b
coorBb[ii] += coorb_new.T

# xyzb_new is the seed for the intersection seam node
# instead of saving the array full of zeros, we just save the entry we know is nonzero
curveProjb[curveName][ii] = xyzb_new.T[np.argmin(dist2)]

# all the remaining seeds in coorb live on the component tri-mesh...
coorb += cb

# now we only have the intersection seam...
nFeature = len(self.nElems)
intNodes = self.seamDict["intNodes"]
seamConn = self.seamDict["seamConn"]
curveSizes = self.seamDict["curveSizes"]

# seeds for the original intersection
intNodesb = np.zeros((N, intNodes.shape[0], intNodes.shape[1]))

# loop over each feature and propagate the sensitivities
curInd = 0
curSeed = 0
for i in range(nFeature):
# just use the same number of points *2 for now
nNewElems = self.nElems[i]
nNewNodes = nNewElems + 1
coor = intNodes
barsConn = seamConn[curInd : curInd + curveSizes[i]]

curInd += curveSizes[i]
method = "linear"
spacing = "linear"
initialSpacing = 0.1
finalSpacing = 0.1

for ii in range(N):
newCoorb = intBar[ii, curSeed : curSeed + nNewNodes, :]
# re-sample the curve (try linear for now), to get N number of nodes on it spaced linearly
# Call Fortran code. Remember to adjust transposes and indices
newCoor, newBarsConn, cb = self.utilitiesAPI.remesh_b(
nNewElems, coor.T, newCoorb.T, barsConn.T + 1, method, spacing, initialSpacing, finalSpacing
)
intNodesb[ii] += cb.T

curSeed += nNewNodes

# add the contributions from the curve projection if we have any
for curveName, v in curveProjb.items():
# get the index
idx = self.seamDict[curveName]["projPtIndx"]
if self.debug:
for ii in range(N):
intNodesb[ii, idx] += v[ii]

dummyConn = np.zeros((0, 4))
barsConn = self.seamDict["barsConn"]
parentTria = self.seamDict["parentTria"]

# do the reverse intersection computation to get the seeds of coordinates
for ii in range(N):
cAb, cBb = self.intersectionAPI.computeintersection_b(
self.compA.nodes.T,
self.compA.triConnStack.T,
dummyConn.T,
self.compB.nodes.T,
self.compB.triConnStack.T,
dummyConn.T,
intNodes.T,
intNodesb[ii].T,
barsConn.T + 1,
parentTria.T + 1,
self.distTol,
)

coorAb[ii] += cAb.T
coorBb[ii] += cBb.T

# get the total sensitivities from both components
compSens_local = {}
compSensA = self.compA.DVGeo.totalSensitivity(coorAb, "triMesh")
for k, v in compSensA.items():
compSens_local[k] = v

compSensB = self.compB.DVGeo.totalSensitivity(coorBb, "triMesh")
for k, v in compSensB.items():
compSens_local[k] = v

# finally sum the results across procs if we are provided with a comm
if comm:
compSens = {}
# because the results are in a dictionary, we need to loop over the items and sum
for k, v in compSens_local.items():
compSens[k] = comm.allreduce(v, op=MPI.SUM)
else:
# we can just pass the dictionary
compSens = compSens_local

return compSens

def associatePointsToSurface(self, points, ptSetName, surface, surfaceEps):
projDict = {}
# Compute the distance from each point to this surface and store it in projDict
if surface in self.compA.triConn:
self._projectToComponent(points, self.compA, projDict, surface=surface)
comp = "compA"
elif surface in self.compB.triConn:
self._projectToComponent(points, self.compB, projDict, surface=surface)
comp = "compB"
else:

# Identify the points that are within the given tolerance from this surface
# surfaceInd contains indices of the provided points not the entire point set
surfaceDist = np.sqrt(np.array(projDict["dist2"]))
surfaceInd = [ind for ind, value in enumerate(surfaceDist) if (value < surfaceEps)]

# Output the points associated with this surface
if self.debug:
data = [np.append(points[i], surfaceDist[i]) for i in surfaceInd]
tecplot_interface.write_tecplot_scatter(
f"{surface}_points_{self.comm.rank}.plt", f"{surface}", ["X", "Y", "Z", "dist"], data
)

# Save the indices only if there is at least one point
if surfaceInd:
self.projData[ptSetName][comp]["surfaceInd"][surface] = surfaceInd
# Initialize a data dictionary for this surface
self.projData[ptSetName][surface] = {}