Get the derivative of the dXpt/dDV and the DVConstraint.

[KapilKhanal]

Hi, I am new to pygeo. The provided jacobian from pygeo is in CSR format for TACS. How can I just get the jacobian matrices for each of the pointset and the DVs and DVconstraint. The tutorial for this is lacking compared to others.

Context: To pass the DXpt/dDV to openMDAO in the compute_partial method.

[marcomangano]

Hello Kapil,
It looks like you have already taken a look at the computeTotalJacobian function. That method returns exactly the derivative you need. You can pass the desired pointset as input. CSR is just the sparse matrix format used internally to not store and perform operations on a dense matrix.

If you need to convert that matrix to dense format, you could try with this Scipy function.
Since you need to use this Jacobian in OM, you might rather want to convert the Jacobian in a dictionary format. We have a function for that that you could use as a reference here. That however works specifically on the output of TotalSensitivity (which multiplies (DXpt/dDV)^T by another Jacobian matrix) so it will NOT return the Jacobian you need. You could write a similar function that converts the output of computeTotalJacobian into a dictionary by sorting the columns into respective DV entries.

I hope this helps and does not make you more confused. Tagging @hajdik since she has worked on the MPhys wrapper - which as far as I understand does not explicitly return DXpt/dDV for forward-mode calculations.

[KapilKhanal]

Thank you for the prompt reply.

I did look into the computeTotalJacobian() function but the shape of it even after changing is not what I thought it would be. I am most likely not understanding some fundamentals of pygeo deformation and the relation between pointsets and vertices. To add to the cylinder tutorial in the documentation..this is what I am trying to do..get dX/dDV separately. Would you be able to help me for the cylinder in the pygeo tutorial?

`def deform_cylinder(local_perturb,const_offset_base=0):
FFDfile = "ffdbox.xyz"
DVGeo = DVGeometry(FFDfile)

stlmesh = mesh.Mesh.from_file("cylinder.stl") #read same base cylinder 
# Add the STL mesh vertices as pointsets to the FFD geometry
DVGeo.addPointSet(stlmesh.v0, "mesh_v0")
DVGeo.addPointSet(stlmesh.v1, "mesh_v1")
DVGeo.addPointSet(stlmesh.v2, "mesh_v2")


# Add a local shape design variable to modify the FFD volume
# Perturbations will be applied only to the y (up-down) direction

DVGeo.addLocalDV("local_perturb", lower=0, upper=4, axis="y", scale=1.0)
#DVGeo.addLocalDV("const_offset_base", lower=0.5, upper=10, axis="y", scale=1.0)

coords = DVGeo.getLocalIndex(0)
inds = coords[:, 0, :]
inds2 = coords[:, 1, :]


# Set the local_perturb as the design variable
currentDV = DVGeo.getValues()["local_perturb"]
newDV = currentDV.copy()

#front and read
lower_front_idx = DVGeo.getLocalIndex(0)[:, 0, 0]
lower_rear_idx = DVGeo.getLocalIndex(0)[:, 0, 1]
upper_front_idx = DVGeo.getLocalIndex(0)[:, 1, 0]
upper_rear_idx = DVGeo.getLocalIndex(0)[:, 1, 1]

# Define perturbations for the lower points (front and rear) of the cylinder
const_offset = const_offset_base * np.ones(10)
for idx in [lower_front_idx, lower_rear_idx]:
    newDV[idx] = const_offset + np.sin(np.linspace(0, local_perturb * np.pi, 10)) / 20  + np.linspace(0.05, -0.30, 10)

# Define perturbations for the upper points (front and rear) of the cylinder
const_offset = -const_offset_base * np.ones(10)
for idx in [upper_rear_idx, upper_front_idx]:
    newDV[idx] = const_offset + np.cos(np.linspace(0, local_perturb * np.pi, 10)) / 20 + np.linspace(0.05, -0.10, 10)

# Set the new design variable values and update the FFD control points
DVGeo.setDesignVars({"shape": newDV.copy()})
   # Save the deformed STL mesh
stlmesh.vectors[:, 0, :] = DVGeo.update("mesh_v0")
stlmesh.vectors[:, 1, :] = DVGeo.update("mesh_v1")
stlmesh.vectors[:, 2, :] = DVGeo.update("mesh_v2")

DVGeo.computeTotalJacobian('mesh_v0')
DVGeo.computeTotalJacobian('mesh_v1')
DVGeo.computeTotalJacobian('mesh_v2')
jac0 = DVGeo.JT["mesh_v0"].todense().T #shoulde be shape (N,#dv)..dX/dDV
print(f"shape of jac0 {jac0.shape}")
assert jac0.shape == (len(stlmesh.v0))
jac1 = DVGeo.JT["mesh_v1"].todense().T #dY/dDV
jac2 = DVGeo.JT["mesh_v2"].todense().T #dZ/dDV

return {'mesh':stlmesh,'jac':[jac0,jac1,jac2]} # deformed mesh and the associated jacobians..shape(N,3)`

[KapilKhanal]

@marcomangano any advice for the tutorial?