 Obtaining the maximum strain energy over the history

Hello,

I have the following weak form of the pde:

`R_z = y*2*z*(psi_max)*dx+ k*(y*(z-1)/epsilon + epsilon*inner(grad(z),grad(y)))*dx`

Here z is the trial function, y is the test function and k,epsilon are material parameters. psi_max is the maximum strain energy over the entire simulation history (I am applying load in steps). Does anyone have any ideas on how I can keep track of the maximum strain energy (at quadrature points essentially)?

I suppose I need to update the values of psi_max after each load step, but I am not quite sure how to do that.

Thank you!

Quadrature elements are currently a bit buggy, but can be used by using some deprecated functionality as follows:

``````from dolfin import *
from ufl import Max

# The following shuts off a deprecation warning for quadrature representation:
import warnings

# A consistent quadrature degree must be used throughout, to avoid mismatch

mesh = UnitSquareMesh(10,10)

# Function space with DoFs at quadrature points:
Q = FunctionSpace(mesh, Qelement)

# Function to store previous max strain energy:
psi_max_old = Function(Q)

# Placeholder, to be replaced with whatever the expression for strain
# energy is in the problem:
psi = SpatialCoordinate(mesh)

# Use this in the variational problem where psi_max is needed:
psi_max = Max(psi,psi_max_old)

# After solving in each load step, update psi_max_old.  When the unknowns of
# a variational problem are in a Quadrature-type space, you need to use the
psi_max_old.assign(project(psi_max,Q,
form_compiler_parameters

# Plot (projection onto linears) to verify:
from matplotlib import pyplot as plt
plot(project(psi_max_old,FunctionSpace(mesh,"Lagrange",1),
plt.show()
``````

Also check out Jeremy Bleyer’s plasticity tutorial, which uses quadrature spaces:

1 Like

Thank you David for your reply. I tried what you suggested. The updating of psi_max_old was really slow, so I tried to do the local projection as mentioned in the Jeremy Bleyer’s tutorial:

``````def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return

#Call local_project function to update psi_max_old
local_project(psi_max, Q, psi_max_old)
``````

but this gives me a compilation error. Do you have any ideas on what could be going wrong?

Thanks again!

``````parameters["form_compiler"]["representation"] = 'quadrature'
``````

after importing `dolfin` it should work. The disadvantage of this is that quadrature representation will be used for all forms, and is less robust than the default (`"uflacs"`) for very complicated forms. I usually try to use quadrature representation only where necessary.

Another trick I’ve found useful is to use an iterative method with a Jacobi preconditioner in a global L^2 projection to quadrature spaces, such that the preconditioner is the inverse of the (diagonal) matrix.

1 Like

Hi kamensky!

would you please tell me how you “usually try to use quadrature representation only
where necessary”?
For instance I want to use

``````import dolfin as df
df.parameters["form_compiler"]["representation"] = "uflacs"
``````

for my global problem. But for any projection (involving quadrature spaces) I would like to
specify

``````df.parameters["form_compiler"]["representation"] = "quadrature"
``````

in the `df.LocalSolver` of the `local_project` function defined as above.

Do you know how to do that?

Thanks,
Philipp

Maybe try something like this:

``````def local_project(v, V, u):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dx
b_proj = inner(v, v_)*dx
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
solver.solve_local_rhs(u)
parameters["form_compiler"]["representation"] = 'uflacs'

local_project(psi_max, Q, psi_max_old)
``````

Ah okay. I thought there might be an option to specify form compiler parameters to
the LocalSolver directly, but I guess this should also work.

Thanks!