Hello, I would like to be able to access the adjoint variable (solution of the adjoint system).

I calculate the gradient \frac{dJ}{dm} by using the command `compute_gradient(J, m)`

. I’m assuming that the gradient calculation is done using the adjoint approach, so following the notation in the dolfin-adjoint mathematical background of adjoints, the adjoint equation is solved (\frac{\partial F(u,m)}{\partial u})^* \lambda=\frac{\partial J}{\partial u}^* and the result is used to calculate \frac{dJ}{dm}^* = -\frac{\partial F}{\partial m}^* \lambda + \frac{\partial J}{\partial m}^*. I’d like to be able to get the variable \lambda.

I’ve tried using `J.block_variable.adj_value`

, but I’m not sure that I’m using it correctly because I get back a value of 1, when it should be a vector.

In this simple example, I’m using dolfin-adjoint 2019.1.0 and solving a linear elasticity pde.

```
from fenics import *
from fenics_adjoint import *
import numpy as np
from ufl import nabla_div
# Define parameters
nelx = 30
nely = 15
nelz = 15
lx = float(nelx)
ly = float(nely)
lz = float(nelz)
E = 1
v = 0.3
mu = Constant(E / (2*(1+v)))
lmbda = Constant(E*v / ((1+v)*(1-2*v)))
# Create a mesh
mesh = BoxMesh.create(
[Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
[nelx, nely, nelz], # number of elements in each direction
CellType.Type.hexahedron)
mesh = Mesh(mesh)
# FE function space and functions
V = VectorFunctionSpace(mesh, "CG", 1)
D = FunctionSpace(mesh, "DG", 0)
u = TrialFunction(V)
v = TestFunction(V)
u_sol = Function(V) # function for storing the solution
density = Function(D) # density field
# Bounday condition
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.)
left = Left()
bc = DirichletBC(V, Constant((0., 0., 0.)), left)
# Stress and strain functions
def sigma(u):
n = len(u) # size of u
return lmbda*nabla_div(u)*Identity(n) + 2*mu*epsilon(u)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# Define density field
density.vector()[:] = np.ones(nelx*nely*nelz)
# Define forms
f = Constant((0., 0., -1.))
L = dot(f, v) * dx
a = inner(density*sigma(u), nabla_grad(v)) * dx
# Solve
A, b = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u_sol.vector(), b)
# Get derivative
J = assemble(action(L, u_sol)) # define compliance
m = Control(density) # define density as the design variable
dJ_dm = compute_gradient(J, m)
# Get adjoint solution
adjoint_sol = J.block_variable.adj_value
print(type(adjoint_sol))
print('adjoint value:', adjoint_sol)
```