Derivatives at the quadrature points

Question

Given a Function u in a FunctionSpace, I would like to obtain the associated derivative as a Function defined on the quadrature points. This may currently be done by using project. However it is far too expensive with respect to my needs, especially when u is a vector with tensor-valued derivatives.

As an example, consider

mesh = UnitSquareMesh(10,10)
Vu = FunctionSpace(mesh, 'CG', 1) 
Vq = VectorFunctionSpace(mesh, 'Quadrature', 0) 
u = Function(Vu)

u = interpolate(Expression("pow(x[0],2)"),Vu)
eps = grad(u)

The currently working method to retrieve a function with eps at the quadrature point is

project(eps,Vq)

Something like interpolate(eps,Vq) does not work because eps is not a GenericFunction or Expression. Is there any other efficient method to get the value of eps by simple and unprocessed evaluation of the derivative of the shape functions at the Gauss points? Of course, I want to avoid iteration loops though elements and Gauss points in python.

I known that there exists the recently added LocalSolver to have a cheap version of 'project'. I am currently using this approach. But it is still somewhat unnecessary complex and undirect.

P.S.: I need this because I am looking at this because I am trying to develop a python version of a solver for plasticity (I known there exists https://bitbucket.org/fenics-apps/fenics-solid-mechanics in C++).

Answer 1

Not sure about the LocalSolver, but here is a suggestion that is very fast if you want repeated evaluations of project(grad(u), Vq).

from dolfin import *

mesh = UnitSquareMesh(100,100)
Vu = FunctionSpace(mesh, 'CG', 1) 
Vq = VectorFunctionSpace(mesh, 'Quadrature', 0) 
u = Function(Vu)

u = interpolate(Expression("pow(x[0],2)"),Vu)
eps = grad(u)

# This is the result you're after
v0 = project(eps, Vq)

# Now we'll do it more efficiently
# You want to solve this variational form
# inner(uq, vq)*dx = inner(vq, grad(u))*dx
uq = TrialFunction(Vq)
vq = TestFunction(Vq)

# The mass matrix is diagonal. Get the diagonal 
M = assemble(inner(uq, vq)*dx)
ones = Function(Vq)
ones.vector()[:] = 1.
Md = M*ones.vector()

# If you only want to assemble right hand side once:
grad_eps = assemble(inner(vq, grad(u))*dx)

# For repeated projects it is much faster to compute the rhs 
# like this, where P can be preassembled just once
#P = assemble(inner(vq, grad(TrialFunction(Vu)))*dx)
#grad_eps = P*u.vector()

# solve
v1 = Function(Vq)
v1.vector().set_local(grad_eps.array() / Md.array())

# Check that solutions match
print all(abs(v0.vector().array() - v1.vector().array()) < 1e-8)

Comments

Thank you for the very nice and smart trick. It works even in parallel.
This was not exactly what I was looking for but surely an efficient solution. Do you have an idea why this is not the default implementation of project in FEniCS?

---

It requires a preassembled matrix and as such it is quite demanding in terms of memory. It is very fast, though, so if you can afford it use it:-) And the mass matrix is only diagonal for some FunctionSpaces.

Answer 2

LocalSolver will be much more efficient for this than a projection. You can use a UFC function for computing derivatives at arbitrary points, but it's not very fast.

In either case, you should not build a Function that holds all the derivatives at all points. It uses too much memory to be useful. A development branch of FEniCS Solid Mechanics (https://bitbucket.org/fenics-apps/fenics-solid-mechanics/branch/garth/restructure) that avoids global storage will be merged into the Solid Mechanics master branch soon. The approach in the Solid Mechanics Apps master branch stores the derivatives, stress update and tangent globally, but is very slow and can only solve small 2D problems and very small 3D problems because the of excessive memory requirements.

Comments

Thanks Garth, I look forward for the new version of fenics-solid-mechanics.

My final goal is to couple plasticity with a gradient damage model. I am hesitating in developing it on the top of fenics-solid-mechanics in c++ or independently in python.

Do you think that in python would be possible to do something efficient, eventually using some compile_extension_module? Or do you definitively suggest using c++ interface?

Just for info, here a minimal demonstrative code for the gradient damage model

Regarding the global storage of the tensor valued field, do you not need to have them somehow if you want to save the plastic deformations to a file? It sounds odd to me that storing six values for gauss point with linear elements would be a problem for memory in a cluster with 2Gb/core memory ratio.

---

I don't think it's possible to make an efficient solver for plasticity, even more so gradient plasticity, in pure Python. What I think is possible and desirable is a Python interface that allows a user to provide certain functions, e.g. the yield function. This could be done with some Python by wrapping a few of the Solid Mechanics virtual functions. Gradient damage is far more tractable than plasticity in Python because the stress update is much simpler.

The new branch of the Solid Mechanics app stores the minimum required, which for classical plasticity is the plastic strain and history variables. The old implementation computed strain etc, and computed the stress update and tangent for all points and stored them - this was too expensive memory-wise for large problems. The new branch computes the stress update and tangent cell-by-cell inside the assembly loop. Using Eigen for the return mapping algorithm, the code is expressive and very fast.

Linear elements generally perform poorly for plastic problems. For a P2 element, say you have 4 quadrature points and you store the strain, stress and tangent, then you are storing 466(66) = 5184 doubles per cell. It starts to hurt.

Kristian Oelgaard has done some gradient plasticity with FEniCS. It's in his thesis, which will be available at the end of the year. You could email him to ask for a copy.

---

Thanks a lot for you precious suggestions.