postprocessing on refined mesh


205
views
0
3 months ago by
K D  
I solve a nonlinear elastic problem with mesh refinement.  I'm trying to then postprocess the solution, in particular I want to compute the elastic energy density using the computed solution.  I'm able to do this easily on the original mesh but am not sure how to do that on the refined mesh.  Also, when I save and open in Paraview, I only see the original 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
Community: FEniCS Project
1

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.

written 3 months ago by Alexander G. Zimmerman  
Thanks Alexander.  I've now added an MWE.
written 3 months ago by K D  

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.

written 3 months ago by Alexander G. Zimmerman  

1 Answer


3
3 months ago by

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 << yto 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 Alexander.  About function spaces: V is a VectorFunctionSpace, and W is a (scalar) FunctionSpace.  The deformation lives in V and the energy density lives in W.  It's not my intention to project the energy density back to the coarse grid.  I just don't know how to compute it on the fine grid.

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.
written 3 months ago by K D  
I figured out my problem.  I was thinking too much like a mathematician and not like a programmer.  I assumed that $psi$ was *defined* as a function of $y$ when I first wrote it out "# Elasticity model parameters".  I guess it was only evaluated there.  So I had to re-evaluate it after I got the final converged solution.  This code snippet after the solve step does it:
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.
written 3 months ago by K D  

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!

written 3 months ago by Alexander G. Zimmerman  
Yes, your answer led me to the correct version which I have posted.  Thanks very much!
written 3 months ago by K D  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »