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.
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
from ffc.quadrature.deprecation \
import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("ignore", QuadratureRepresentationDeprecationWarning)
# A consistent quadrature degree must be used throughout, to avoid mismatch
# between the quadrature space and quadature of the variational form:
QUAD_DEG = 4
dx = dx(metadata={"quadrature_degree":QUAD_DEG})
mesh = UnitSquareMesh(10,10)
# Function space with DoFs at quadrature points:
Qelement = FiniteElement("Quadrature", mesh.ufl_cell(),
degree=QUAD_DEG,
quad_scheme="default")
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
# deprecated quadrature representation.
psi_max_old.assign(project(psi_max,Q,
form_compiler_parameters
={"quadrature_degree":QUAD_DEG,
"representation":"quadrature"}))
# Plot (projection onto linears) to verify:
from matplotlib import pyplot as plt
plot(project(psi_max_old,FunctionSpace(mesh,"Lagrange",1),
form_compiler_parameters={"quadrature_degree":QUAD_DEG}))
plt.show()
Also check out Jeremy Bleyer’s plasticity tutorial, which uses quadrature spaces:
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?
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.
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