### postprocessing on refined mesh

I've pasted an MWE that shows me the initial coarse mesh on Paraview and the post-processing to find the energy density $psi$ is also done on the coarse space.

""""

MWE for elasticity

How to plot and work with solution on adaptively-refined mesh??

"""

from fenics import *

# Optimization options for the form compiler

parameters["form_compiler"]["cpp_optimize"] = True

ffc_options = {"optimize": True, \

"eliminate_zeros": True, \

"precompute_basis_const": True, \

"precompute_ip_const": True}

# Create mesh and define function space

mesh = RectangleMesh(Point(-1,-1),Point(1,1), 4, 4)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

W = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary conditions

F11=1.0

F12=0.0

F21=0.0

F22=1.0

y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"),F11=F11,F12=F12,F21=F21,F22=F22,degree=2)

def x0_boundary(x, on_boundary):

return on_boundary

bc = DirichletBC(V, y_BC, x0_boundary)

# Define functions

w = TrialFunction(V)

v = TestFunction(V)

y = Function(V,name="deformation")

#initial guess

y_reference = Expression(("x[0]","x[1]"),degree=2)

y=project(y_reference,V)

# Kinematics

d = y.geometric_dimension()

I = Identity(d) # Identity tensor

F = I + grad(y) # Deformation gradient

C = F.T*F # Right Cauchy-Green tensor

Ic = tr(C)

J = det(F)

# Elasticity model parameters

E, nu = 10.0, 0.3

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

b=Expression(("10*exp((-x[0]*x[0]-x[1]*x[1])/0.1)","0"),degree=2) # body force

psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 # strain energy density (compressible neo-Hookean)

# Total potential energy

Pi = psi*dx - inner(b,y)*dx

F = derivative(Pi, y, v) # Compute first variation of Pi

J = derivative(F, y, w) # Compute Jacobian of F

# Refinement parameters

M = inner(y,y)*dx()

tol = 1.e-5

# Solve variational problem

solve(F == 0, y, bcs=bc, J=J,tol=tol,M=M)

# Save solution in PVD format

fileYpvd = File("deformation.pvd")

fileYpvd << y

psi_field=Function(W,name="energy_density")

psi_field = project(psi,W)

filePsipvd = File("Energy_Density.pvd")

filePsipvd << psi_field

### 1 Answer

When I run your MWE (after fixing the indentation in the function def), it outputs both fields (deformation and energy density) on the coarse grid, and that's what I see in ParaView.

When I change `fileYpvd << y`

to `fileYpvd << y.leaf_node()`

, then this outputs the deformation field on the adapted grid, which I see in ParaView.

It is unclear to me why there is another FunctionSpace `W`

onto which you project the energy density field. I suppose that your intention is indeed to project the solution back onto the coarse grid; so I suppose the example already works as expected in this case.

Thanks for the information about deformation on the adapted grid. That's very helpful. However, the other part of my difficulty is to now compute $psi$ on the adapted grid.

```
solve(S == 0, y, bcs=bc, J=J,tol=tol,M=M)
y=y.leaf_node()
mesh=mesh.leaf_node()
W_refined=FunctionSpace(mesh.leaf_node(),"Lagrange",1)
psi_field = project((mu/2)*(tr((I + grad(y)).T*(I + grad(y))) - 3) - mu*ln(det(I+grad(y))) + (lmbda/2)*(ln(det(I+grad(y))))**2,W_refined)
```

I'm sure there is a way to define a function $psi$ to avoid these long expressions.

Very good. So is your question answered? If the answer we're discussing here doesn't suffice, but you know the answer now, then it'd be good to post the correct answer (yourself) and accept the correct answer.

I'm glad you're making progress!

It's possible that I've fixed an issue similar to yours before, but it's really hard for me to tell without an example that runs and demonstrates the issue.

As a wild guess, I would suggest you might want to try swapping

`W_refined`

with`W_refined.leaf_node()`

and`psi_field`

with`psi_field.leaf_node()`

. At least I've often been perplexed by what I see in ParaView, only to realize that I missed a`.leaf_node()`

call somewhere.Thanks for the MWE!

Side note: I am pretty sure that you can remove the line

`psi_field=Function(W,name="energy_density")`

before the line`psi_field = project(psi,W)`

, because the latter returns a new instance (breaking the reference to the old one) anyway.