# 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)[0]

# 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!

Hi, your discussion is very helpful to me! But Iâ€™m not sure how to achieve it in the new version of dolfinx. For example, how to update the following line

``````psi_max_old.assign(project(psi_max,Q,
form_compiler_parameters
``````

Thank you!

hi,
you can now directly interpolate UFL expressions, see here:
https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html#stress-computation

``````psi_max_old.interpolate(psi_max)
``````

But it generated an error:

``````RecursionError: maximum recursion depth exceeded in __instancecheck__
``````

Do you have any suggestions? Thanks again!

Please provide a code that can reproduce this error, as it is impossible to guide you without the context of the error message.

Hi dokken, here is a code reproducing this recursion error.

``````import numpy as np
import dolfinx.fem as fem
import dolfinx.mesh as mesh
import ufl
from mpi4py import MPI

domain = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=[[0,0], [1,1]],
ndim = domain.geometry.dim
fdim = ndim-1

element_u = ufl.VectorElement('Lagrange', domain.ufl_cell(), degree=1, dim=2)
V_u = fem.FunctionSpace(domain, element_u)
element_alpha = ufl.FiniteElement('Lagrange', domain.ufl_cell(), degree=1)
V_alpha = fem.FunctionSpace(domain, element_alpha)

# Measures

u = fem.Function(V_u)
# Function to store previous max value
psi_max_old = fem.Function(V_alpha)
# Use this in the variational problem where psi_max is needed
psi_max = ufl.max_value(psi, psi_max_old)

for i_t, t in enumerate(np.linspace(0, 0.01, 1)):
# ... solve the problem ...
psi_max_old.interpolate(psi_max)
``````

I guess that I made a mistake with the function space of psi_max_old. But I have no idea how to get through.

Thanks!

As seen in the example pointed out by @bleyerj, you need to wrap the `psi_max` as a `dolfinx.fem.Expression` before interpolating

Thanks for your help! I wrapped it up as follows and it worked.

``````V_psi_max = fem.FunctionSpace(domain, ("Lagrange", 1))