Sensitivities of a modal eigenvalue

dolfin-adjoint
1

Dear,

I’m working with topological optimization and would like to analyze situations involving dynamics problems.

To begin with, analyzing a problem of: Free Vibration and eigenvalue problem.

According to Topology Optimization | SpringerLink, to maximize eigenvalue we have:

\max_{\rho } \ \beta

s.t.:\ \lambda_{i} \geq \beta ,\ i=1,...,N_{dof,}

\left( K-\lambda_{i} M\right) \Phi_{i} =0,\ i=1,...,N_{dof},

\sum^{N}_{e=1} v_{e}\rho_{e} \leq V,\ \ \ \ 0<\rho_{min} \leq \rho_{e} \leq 1,\ \ e=1,...,N

And the sensitivities of a single modal eigenvalue are simple found as:

\frac{\partial \lambda_{i} }{\partial \rho_{e} } =\Phi^{T}_{i} \left[ \frac{\partial K}{\partial \rho_{e} } -\lambda_{i} \frac{\partial M}{\partial \rho_{e} } \right] \Phi_{i}

So, considering the optimization approach made using the dolfin-adjoint, I need to solve the “forward” function using modal analysis, solving the eigenvalues and eigenvectors and “assemble” the J = “eigenvalues” and deliver to the Jhat = ReducedFunctional( J, m).

But my problem is that I tried the possible “conversions” and I can’t “send” r(eingenvalue) to J.

And I would still like to understand how to insert the sensitivity equation into fenics. The fact that considering the optimization method, the sensitivity equation would be the compliance itself to be calculated in the adjoint method.

I recently read a work, however, I believe that the approach mixed modal analysis with forced analysis and I was even more doubtful.

Could someone help me with this “conversion” that is generated by petsc(r - eigenvalue) to send to the dolfin adjoint?

I’m using the code summarized below for reasons I haven’t published it yet.
https://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.html

My optimization code is based on the following example:
http://www.dolfin-adjoint.org/en/latest/documentation/poisson-topology/poisson-topology.html

Obs.:
if I consider the following situation:

    N_eig = 1   # number of eigenvalues
    eigensolver.solve(N_eig)
   
    # Extraction
    r, c, rx, cx = eigensolver.get_eigenpair(0)
    # r = eigenvalue
    # rx = eigenvector
    
    # Initialize function and assign eigenvector
    eigenvector = Function(V, name = "Eigenvector")
    eigenvector.vector()[:] = rx
    eigenvalue = Function(V, name = "Eigenvalue")
    eigenvalue.vector()[:] = r
...

    J = eigenvalue
    m = Control(rho_n)
    Jhat = ReducedFunctional(J, m)
...

generates the error:

Traceback (most recent call last):
  File "OTE_modal_Dynamic_cantilever_3D_versao00_1_noload.py", line 207, in <module>
    rho_opt = solver.solve()
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/optimization/ipopt_solver.py", line 199, in solve
    results = self.ipopt_problem.solve(guess)
  File "src/cyipopt.pyx", line 581, in cyipopt.problem.solve
  File "src/cyipopt.pyx", line 617, in cyipopt.objective_cb
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/reduced_functional_numpy.py", line 36, in __call__
    return self.rf.__call__(self.set_local(m_copies, m_array))
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/reduced_functional.py", line 137, in __call__
    func_value = self.scale * self.functional.block_variable.checkpoint
TypeError: unsupported operand type(s) for *: 'float' and 'NoneType'
2

Obs.:
if I consider the following situation:

    N_eig = 1   # number of eigenvalues
    eigensolver.solve(N_eig)
   
    # Extraction
    r, c, rx, cx = eigensolver.get_eigenpair(0)
    # r = eigenvalue
    # rx = eigenvector
    
    # Initialize function and assign eigenvector
    eigenvector = Function(V, name = "Eigenvector")
    eigenvector.vector()[:] = rx
    eigenvalue = Function(V, name = "Eigenvalue")
    eigenvalue.vector()[:] = r
...

    J = r
    m = Control(rho_n)
    Jhat = ReducedFunctional(J, m)

generates the error:

Traceback (most recent call last):
  File "OTE_modal_Dynamic_cantilever_3D_versao00_1_noload.py", line 229, in <module>
    rho_opt = solver.solve()
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/optimization/ipopt_solver.py", line 199, in solve
    results = self.ipopt_problem.solve(guess)
  File "src/cyipopt.pyx", line 581, in cyipopt.problem.solve
  File "src/cyipopt.pyx", line 617, in cyipopt.objective_cb
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/reduced_functional_numpy.py", line 36, in __call__
    return self.rf.__call__(self.set_local(m_copies, m_array))
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/tape.py", line 46, in wrapper
    return function(*args, **kwargs)
  File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/reduced_functional.py", line 137, in __call__
    func_value = self.scale * self.functional.block_variable.checkpoint
AttributeError: 'float' object has no attribute 'block_variable'

I believe the problem is in this error:

File "/Users/rafaelferr0/anaconda3/envs/fenics2/lib/python3.8/site-packages/pyadjoint/reduced_functional.py", line 137, in __call__
    func_value = self.scale * self.functional.block_variable.checkpoint

Does anyone know what this error means?

I’m generating a minimal code to improve the discussion, but if anyone identifies the problem I’d be grateful.

@sebastkm @dokken @bleyerj

3

Hi, I don’t use dolfin-adjoint so I cannot help you a lot but I think that the functional J must be a form or some kind of expression. Here r is just a float corresponding to your numerically computed eigenvalue and dolfin-adjoint cannot optimize over this single float. This is why you get the error.

4

Thanks for your return. I can see this situation and I’m trying to change the eigenvalues output to some kind of eigenvalue dependent function, something like an eigenfunction of type J = assemble(F(Lambda)*dx) (I don’t know if that’s right). I’m studying some works but it’s still not clear.

In that case, would you, or anyone else, have any ideas on how to build an eigenvalue dependent function?

5

As @bleyerj pointed out, r is just a float representing your eigenvalue and can’t be passed into dolfin adjoint as it is not an adjoint problem nor does it contain any of the information needed to calculate the gradients. Therefore you cannot use an adjoint method to find the sensitivites.

You will need to calculate the sensitivities explicitly using the equation you provided in your original post.

\frac{\partial \lambda_{i} }{\partial \rho_{e} } =\Phi^{T}_{i} \left[ \frac{\partial K}{\partial \rho_{e} } -\lambda_{i} \frac{\partial M}{\partial \rho_{e} } \right] \Phi_{i}

This involves locally assembling the two matricies \frac{\partial K}{\partial \rho_{e} } and \frac{\partial M}{\partial \rho_{e} }. It is much more efficient to assemble locally and then perform the matrix algebra to find \frac{\partial \lambda_{i} }{\partial \rho_{e} } for each element.

Also note that this can only be done on an elemental formulation, i.e using DG0 as your Function Space for design variables. This is due to the fact that the sensitivity matrices are assembled for each element and would therefore be unable to distinguish the effects of your design variables if they were defined as nodal values.

This paper explains the process: (PDF) Multiscale optimisation of resonant frequencies for lattice-based additive manufactured structures

1 Like
6

@bleyerj or @dokken, can you help on how to perform this calculation of explicit sensitivities?
\frac{\partial \lambda_{i} }{\partial \rho_{e} } =\Phi^{T}_{i} \left[ \frac{\partial K}{\partial \rho_{e} } -\lambda_{i} \frac{\partial M}{\partial \rho_{e} } \right] \Phi_{I}

7

Hello guys!

@Rafael_Ferro, did you manage to solve this problem?

@conpierce8 or @mNightingale could you help with this problem?

Could anyone help us?

8

Dear Silvia,

First, sorry for the late reply. Well, I’m glad you’re working on this issue, but unfortunately I still haven’t found an effective solution. I would love if someone from the community could stop for a while and help us. As I’ve practically been working alone, it’s very difficult to make progress on some things.

I’m working on a minimal code, I believe this week I’ll post it here.

9

Dear,

Below is a MWE. The idea is to find the sensitivity in the eigenvalue and then insert it into the optimization with a constraint (to be discussed later).

from fenics import *
import numpy as np

# Input data:
E, nu, rho_d = 200e9, 0.3, 7860 # E=N/m2 and rho_d=kg/m3
L, H, B = 3.0, 1.0, 0.02 # Geometry of the design domain

# penalisation and SIMP constants
p, eps1 = Constant(3.0), Constant(1.0e-3)

# Mesh constants
mesh = BoxMesh(Point(0., 0., 0.), Point(L, H, B), 30, 10, 2)

# Define function space and base functions
V, W = VectorFunctionSpace(mesh, "CG", 1), FunctionSpace(mesh, "CG", 1)

# Boundary Condition
left = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14) 
bcs = [DirichletBC(V, Constant((0.0, 0.0, 0.0)), left)]

# Forward function adapted to find the eigenvalues:
def forward_FindEigenValues():
    v = TestFunction(V)
    du = TrialFunction(V) 
    mu = E/2./(1+nu)
    lmbda = E*nu/(1+nu)/(1-2*nu)
    def eps(v):
        return sym(grad(v))
    def sigma(v):
        dim = v.geometric_dimension()
        return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)
    """
    a dummy l_form and right-hand side vector have been introduced 
    to call this function
    """
    # Linear Elasticity Bilinear Form:
    k_form = inner(sigma(du), eps(v))*dx 

    l_form = Constant(1.)*v[0]*dx  
    
    K = PETScMatrix()
    b = PETScVector()
    assemble_system(k_form, l_form, bcs, A_tensor=K, b_tensor=b)
    
    m_form = rho_d*inner(du, v)*dx
    M = PETScMatrix()
    assemble(m_form, tensor=M)
    # To Calculate the largest eigenvalue
    #bcs.zero(M)
    # Zero out the rows to enforce boundary condition:
    for bc in bcs:
        bc.zero(M)
    
    # Solve eigenmodes
    eigensolver = SLEPcEigenSolver(K, M)
    eigensolver.parameters['problem_type'] = 'gen_hermitian'
    eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
    eigensolver.parameters['spectral_shift'] = 0.
    
    N_eig = 1   # number of eigenvalues
    eigensolver.solve(N_eig)
    
    eigenvector = []
    eigenvalue = []
   
    # Extraction:
    #for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(0)
    # r = eigenvalue
    # rx = eigenvector
    
    # Initialize function and assign eigenvector:
    eigenvector = Function(V, name = "Eigenvectors")
    eigenvector.vector()[:] = rx
    eigenvalue = Function(V, name = "Eigenvalues")
    eigenvalue.vector()[:] = r
   
    # eigenfrequency
    freq = sqrt(r)/2/pi
    
    return freq, eigenvector


# First run the the eigenvalue analysis
eigenvalue, Eigenmode = forward_FindEigenValues()

rho_n = interpolate(Constant(float(eigenvalue)), W)

def gradient(rho_n):

    # Initilise a zero matrix for the gradients of dimension: num of cells
    num_local_cells = mesh.num_cells()
    dJdsF= np.zeros((num_local_cells,1))

    v = TestFunction(V)  # Test CG_Vector Function Space function
    u = TrialFunction(V)  # Trial CG_Vector Function Space function
    
    #mu = E/2./(1+nu)
    #lmbda = E*nu/(1+nu)/(1-2*nu)
    E_rho = eps1 + (E-eps1)*(rho_n**p)
    lmbda = nu*E_rho/((1+nu)*(1-2*nu))
    mu = E_rho/(2*(1+nu))
    
    def epsg(v):
        return sym(grad(v))
    
    def sigmag(v):
        dim = v.geometric_dimension()
        return lmbda*tr(epsg(v))*Identity(dim) + 2.0*mu*epsg(v)
    
    # pre-assemble the mass matrix
    dm = rho_d*dot(u,v)*dx
    
    # pre-assemble the global stiffness grad matrix 
    dk = inner(sigmag(u), epsg(v))*dx
    
    for cell in cells(mesh): ##Only calculate for the cells in your design space.
        
        # Extract differentiated stiffness and Mass matrix on cell
        cell_id = cell.index()
        cell = Cell(mesh,cell_id)            
        dofs = V.dofmap().cell_dofs(cell_id)
        
        # Extract differentiated stiffness and Mass matrix on cell
        dK = assemble_local(dk, cell)
        dM = assemble_local(dm, cell)
        
        eigvect_e = np.array(Eigenmode[dofs])
        
        Kgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dK),  eigvect_e)))
        Mgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dM),  eigvect_e)))
        gradg =  (Kgrad - eigenvalue*Mgrad)
        
        dJdsF[cell_id] = gradg  # add this gradient into the temp matrix
        
    return dJdsF
    
dJdsF = gradient(rho_n)

See what generates the error below:

Traceback (most recent call last):

  File "/Users/.../OTE_modal_00_discourse.py", line 142, in <module>
    dJdsF = gradient(rho_n)

  File "/Users/.../OTE_modal_00_discourse.py", line 132, in gradient
    eigvect_e = np.array(Eigenmode[dofs])

  File "/Users/r.../site-packages/ufl/exproperators.py", line 449, in _getitem
    all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)

  File "/Users/.../site-packages/ufl/index_combination_utils.py", line 172, in create_slice_indices
    elif ind == Ellipsis:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

However, all dJdsF values are set to zero.

I would really like someone to help me.

10

@dokken, @bleyerj or someone else, could you check if the above code is on the way to finding the solution to the equation below:
\frac{\partial \lambda_{i} }{\partial \rho_{e} } =\Phi^{T}_{i} \left[ \frac{\partial K}{\partial \rho_{e} } -\lambda_{i} \frac{\partial M}{\partial \rho_{e} } \right] \Phi_{I}

11

You can avoid the error message by doing

eigvect_e = np.array(Eigenmode.vector()[dofs])

However, I don’t know if this correctly achieves what you are trying to do.

12

Thanks a lot @conpierce8, it really worked and no more error. However, I would like some fenics specialist to know if the code solves the aforementioned equation of the derivative of the eigenvalue in relation to a variable, which in this case is the density of the material interpolated in the SIMP method.

13

Dear @dokken, could you help me by checking if the code matches the solution of the above equation.

14

Dear, I ended up exchanging the eigenvalue entry for the frequency. Now considering the fix, the following error appears:

  File "/Users/.../OTE_modal_00_discourse.py", line 146, in <module>
    dJdsF = gradient(rho_0)

  File "/Users/.../OTE_modal_00_discourse.py", line 142, in gradient
    dJdsF[cell_id] = gradg  # add this gradient into the temp matrix

ValueError: could not broadcast input array from shape (3,) into shape (1,)

Here is the corrected code:

from fenics import *
import numpy as np

# Input data:
E, nu, rho_d = 200e9, 0.3, 7860 # E=N/m2 and rho_d=kg/m3
L, H, B = 3.0, 1.0, 0.02 # Geometry of the design domain

# penalisation and SIMP constants
p, eps1 = Constant(3.0), Constant(1.0e-3)

rho_0= Constant(0.5) # Initial guess 

# Mesh constants
mesh = BoxMesh(Point(0., 0., 0.), Point(L, H, B), 30, 10, 2)

# Define function space and base functions
V, W = VectorFunctionSpace(mesh, "CG", 1), FunctionSpace(mesh, "CG", 1)

# Boundary Condition
left = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14) 
bcs = [DirichletBC(V, Constant((0.0, 0.0, 0.0)), left)]

# Forward function adapted to find the eigenvalues:
def forward_FindEigenValues():
    v = TestFunction(V)
    du = TrialFunction(V) 
    mu = E/2./(1+nu)
    lmbda = E*nu/(1+nu)/(1-2*nu)
    def eps(v):
        return sym(grad(v))
    def sigma(v):
        dim = v.geometric_dimension()
        return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)
    """
    a dummy l_form and right-hand side vector have been introduced 
    to call this function
    """
    # Linear Elasticity Bilinear Form:
    k_form = inner(sigma(du), eps(v))*dx

    l_form = Constant(1.)*v[0]*dx  
    
    K = PETScMatrix()
    b = PETScVector()
    assemble_system(k_form, l_form, bcs, A_tensor=K, b_tensor=b)
    
    m_form = rho_d*inner(du, v)*dx
    M = PETScMatrix()
    assemble(m_form, tensor=M)
    # To Calculate the largest eigenvalue
    #bcs.zero(M)
    # Zero out the rows to enforce boundary condition:
    for bc in bcs:
        bc.zero(M)
    
    # Solve the beam eigenmodes
    eigensolver = SLEPcEigenSolver(K, M)
    eigensolver.parameters['problem_type'] = 'gen_hermitian'
    eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
    eigensolver.parameters['spectral_shift'] = 0.
    
    N_eig = 1   # number of eigenvalues
    eigensolver.solve(N_eig)
    
    eigenvector = []
    eigenvalue = []
   
    # Extraction:
    #for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(0)
    # r = eigenvalue
    # rx = eigenvector
    
    # Initialize function and assign eigenvector:
    eigenvector = Function(V, name = "Eigenvectors")
    eigenvector.vector()[:] = rx
    eigenvalue = Function(V, name = "Eigenvalues")
    eigenvalue.vector()[:] = r
   
    # eigenfrequency
    freq = sqrt(r)/2/pi
    
    #return freq, eigenvalue, eigenvector
    return eigenvalue, eigenvector


# First run the the eigenvalue analysis
eigenvalue, Eigenmode = forward_FindEigenValues()
#freq, eigenvalue, Eigenmode = forward_FindEigenValues()

#rho_n = interpolate(Constant(float(eigenvalue)), W)
#rho_n = interpolate(Constant(float(freq)), W)
rho_n = interpolate(rho_0, W)

def gradient(rho_n):
#def gradient():

    # Initilise a zero matrix for the gradients of dimension: num of cells
    num_local_cells = mesh.num_cells()
    dJdsF= np.zeros((num_local_cells,1))

    v = TestFunction(V)  # Test CG_Vector Function Space function
    u = TrialFunction(V)  # Trial CG_Vector Function Space function
    
    #mu = E/2./(1+nu)
    #lmbda = E*nu/(1+nu)/(1-2*nu)
    E_rho = eps1 + (E-eps1)*(rho_n**p)
    lmbda = nu*E_rho/((1+nu)*(1-2*nu))
    mu = E_rho/(2*(1+nu))
    
    def epsg(v):
        return sym(grad(v))
    
    def sigmag(v):
        dim = v.geometric_dimension()
        return lmbda*tr(epsg(v))*Identity(dim) + 2.0*mu*epsg(v)
    
    # pre-assemble the mass matrix
    dm = rho_d*dot(u,v)*dx
    
    # pre-assemble the global stiffness grad matrix 
    dk = inner(sigmag(u), epsg(v))*dx
    
    for cell in cells(mesh): ##Only calculate for the cells in your design space.
        
        # Extract differentiated stiffness and Mass matrix on cell
        cell_id = cell.index()
        cell = Cell(mesh,cell_id)            
        dofs = V.dofmap().cell_dofs(cell_id)
        
        # Extract differentiated stiffness and Mass matrix on cell
        dK = assemble_local(dk, cell)
        dM = assemble_local(dm, cell)
        
        #eigvect_e = np.array(Eigenmode[dofs])
        eigvect_e = np.array(Eigenmode.vector()[dofs])
        
        Kgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dK),  eigvect_e)))
        Mgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dM),  eigvect_e)))
        gradg =  (Kgrad - eigenvalue*Mgrad)
        
        dJdsF[cell_id] = gradg  # add this gradient into the temp matrix
        
    return dJdsF
    
dJdsF = gradient(rho_n)

Dear colleagues, please help me.

15

For each eigenpair, you get a scalar eigenvalue and a (n \times 1) eigenvector, where n is the number of degrees of freedom of the system. It is correct to put the real component of the eigenvector (rx) into the function eigenvector, but there is no need to make a separate function for eigenvalues. They are scalar and you can keep them in a list.

Here is the corrected code

from fenics import *
import numpy as np

def epsilon(v):
    return sym(grad(v))
def sigma(v):
    dim = v.geometric_dimension()
    return lmbda*tr(epsilon(v))*Identity(dim) + 2.0*mu*epsilon(v)

# Input data:
E, nu, rho_d = 200e9, 0.3, 7860 # E=N/m2 and rho_d=kg/m3
L, H, B = 3.0, 1.0, 0.02 # Geometry of the design domain
n_eig = 3   # number of eigenvalues


mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# penalisation and SIMP constants
p, eps1 = Constant(3.0), Constant(1.0e-3)

volfrac, penal= Constant(0.5), 3 # Initial guess 

# Mesh constants
mesh = BoxMesh(Point(0., 0., 0.), Point(L, H, B), 30, 10, 2)

# Define function space and base functions
U = VectorFunctionSpace(mesh, "CG", 1)
D = FunctionSpace(mesh, "DG", 0)

u, v = TestFunction(U), TrialFunction(U)
rho = interpolate(volfrac, D)


# Boundary Condition
left = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14) 
bcs = [DirichletBC(U, Constant((0.0, 0.0, 0.0)), left)]

# Linear Elasticity Bilinear Form:
k_form = inner(rho**penal*sigma(u), epsilon(v))*dx
m_form = rho_d*inner(rho*u, v)*dx
l_form = Constant(1.)*v[0]*dx  

dK_form = diff(k_form, rho)
dM_form = diff(m_form, rho)

# Forward function adapted to find the eigenvalues:
def forward_get_eigenpairs(rho, n_eig):
    
    K = PETScMatrix()
    M = PETScMatrix()
    b = PETScVector()
    assemble_system(k_form, l_form, bcs, A_tensor=K, b_tensor=b)
    assemble(m_form, tensor=M)

    # Zero out the rows to enforce boundary condition:
    for bc in bcs:
        bc.zero(M)
    
    # Solve the beam eigenmodes
    eigensolver = SLEPcEigenSolver(K, M)
    eigensolver.parameters['problem_type'] = 'gen_hermitian'
    eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
    eigensolver.parameters['spectral_shift'] = 0.
    
    eigensolver.solve(n_eig)
    
    eigenvectors = np.zeros((U.dim(),n_eig))
    eigenvalues = []
   
    # Extraction:
    for i in range(n_eig):
        # Extract eigenpair
        r, c, rx, cx = eigensolver.get_eigenpair(i)
        eigenvectors[:,i] = rx        
        eigenvalues.append(r)
    return eigenvalues, eigenvectors


# First run the the eigenvalue analysis
eigenvalue, eigenmode = forward_get_eigenpairs(rho, n_eig)

def gradient(rho,n_eig):
    # Initilise a zero matrix for the gradients of dimension: num of cells
    num_local_cells = mesh.num_cells()
    dJdsF= np.zeros((num_local_cells,n_eig))

    for i in range(n_eig):
        for cell in cells(mesh): ##Only calculate for the cells in your design space.
            
            # Extract differentiated stiffness and Mass matrix on cell
            cell_id = cell.index()
            cell = Cell(mesh,cell_id)            
            dofs = U.dofmap().cell_dofs(cell_id)
            
            # Extract differentiated stiffness and Mass matrix on cell
            dK = assemble_local(dK_form, cell)
            dM = assemble_local(dM_form, cell)
            
            eigvect_e = eigenmode[dofs,i]
            
            Kgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dK),  eigvect_e)))
            Mgrad = (np.dot(eigvect_e.T,  np.dot(np.array(dM),  eigvect_e)))
            gradg =  (Kgrad - eigenvalue[i]*Mgrad)
            
            dJdsF[cell_id,i] = gradg  # add this gradient into the temp matrix
        
    return dJdsF
    
dJdsF = gradient(rho,n_eig)

Can you please have a look at it @mNightingale?