I am solving some linear second order PDE using first order Lagrangian elements. Now, I would like to compute some physical quantities that involve second order derivatives, the Laplacian for instance.
Consider the following example:
from dolfin import *
from mshr import *
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Expression('-logf(sqrt(pow(x[0]+1, 2) + pow(x[1]+1, 2)))', degree = 1)
u = project(u,V)
Delta_u = project(div(grad(u)),V)
Here, u is chosen to be a harmonic function (apart from the point (-1,-1)).
However, this won’t work. Aparently, fenics is working on each triangle separately, so any second order derivative will vanish for this piecewise linear function u.
Is there a canonical way in fenics to compute derivatives which are of higher order than the underlying function space?
the behaviour you describe is apparently due to the fact, that the Expression has degree = 1, thus is interpolated into piecewise linear function (with vanishing second derivative).
You have basically two options:
Increase degree of expression
This way you will make the polynomial approximation more accurate and also have second derivative non zero. But keep in mind, that many functions are very badly approximated with polynomials (as the one with singularities…)
Rewite the Expression into UFL language
x = SpatialCoordinate(mesh)
u = log(sqrt((x[0] + 1) ** 2 + (x[1] + 1) ** 2))
Delta_u = project(div(grad(u)), V)
Please note the difference. The second approach is much better in most cases, since you will be evaluating the symbolic true nonlinear function (log and power…) at quadrature points in order to assemble the projection matrix. These symbols log, sqrt, ... are coming from FEniCS symbolic language - UFL.
The first approach might be favourable if for some reason you really need to interpolate expression into polynomial finite element space - but whenever you have an expression in form or projection - go with symbolic UFL language.
thank you for your answer! Indeed, with the functions from UFL this is working. Remark: In my version the (natural) logarithm is implemented in ufl.ln, not in ufl.log.
However, I must admit that my example was a bit misleading. The point is that my function u is the result of a finite element computation based on first order Lagrangian elements. So, u is given as a piecewise linear function, and this is my starting point.
It seems that both of your suggestions won’t help me here?