Computing element-wise spatial derivatives of the solution

General
1

Currently, I am solving a scalar linear PDE using first-order basis functions. I need to compute element-wise 2nd order derivatives to generate the hessian of the solution. To this end, I currently project the solution to 2nd order basis functions and then would like to compute the derivates. The projection problem that I solve is given by:

#msh represents the mesh and using fenicsx.

wh = fem.FunctionSpace(msh, ("Lagrange", 2)) # higher order function space 
bcp = fem.dirichletbc(value=ScalarType(0.0), dofs=dofs, V=wh) #boundary condition of the problem

up = ufl.TrialFunction(wh)
vp = ufl.TestFunction(wh)

ap =  inner(up,vp) * dx #L2 projection bilinear form
Lp =  uh* vp * ufl.dx # load vector computed using first-order solution uh

proj_problem = fem.petsc.LinearProblem(ap, Lp, bcs=[bcp], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
up = proj_problem.solve()

#computing the hessian
ux1p = ufl.diff(up,x1)
ux1x1p = ufl.diff(ux1p,x1)

ux2p = ufl.diff(up,x2)
ux2x2p = ufl.diff(ux2p,x2)
ux1x2p = ufl.diff(ux1p,x2)

#to compute ux1x1
W = fem.FunctionSpace(msh, (“DG”, 0))
vdg = ufl.TestFunction(W)

ap = inner(ux1p,vdg) * dx
Mp = fem.form(ap)
dx1x1 = fem.petsc.assemble_vector(Mp) # I plan to divide this by the volume of the element as 
                                                                   # element wise hessian will be constant.

but when I run this exactly at ap = inner(ux1p,vdg) * dx

I get the following error:

This integral is missing an integration domain.
ERROR:UFL:This integral is missing an integration domain.
Traceback (most recent call last):
File “/Users/ankit/Desktop/fenicscodes/Exponential-integrators/conv_diff_steady_state/conv_diff_clean.py”, line 287, in
fenics_solver(filename,connect,coord)
File “/Users/ankit/Desktop/fenicscodes/Exponential-integrators/conv_diff_steady_state/conv_diff_clean.py”, line 255, in fenics_solver
ap = inner(ux1p,vdg) * dx
File “/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ufl/measure.py”, line 432, in rmul
error(“This integral is missing an integration domain.”)
File “/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ufl/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.

can someone kindly tell me how to compute spatial derivatives of a solution and its projection to a richer space?

Regards
Ankit

2

You haven’t made this example fully reproducible, as several definitions are missing. Please make sure that the error is reproducible.

If I assume x1,x2 = ufl.SpatialCoordinate(msh)
then you are saying that you want to solve the problem
v_i = \frac{\partial^2 u_p}{\partial^2 x_i} \text{ in } \Omega, i=1,\dots,d, \Omega\in\mathbb{R}^d, which should be integrated by parts (on the rhs, as the double derivative of a linear function is 0).

3

Since u_p is a projection of u_h to second order basis obtained by solving a global projection problem, Doesn’t it make the derivative non-zero? I want to solve the problem as v_i = \frac{\partial^2 u_p}{\partial x^2} but rather than solving a linear system, I was looking more towards computing the functional \int \frac{\partial^2 u_p}{\partial x^2} w_{0} \, dx where w_0 is DG-0 basis function for every element in the mesh which I want to follow by dividing it with the volume.

Is there any other way to compute the second derivatives of a projection of a first-order approximate solution for every element?

I need these derivatives to drive an anisotropic Taylor series-based mesh adaptation and hence cant use second-order elements because, in that case, I will need to approximate third-order derivatives for the mesh adaptation.

Regards
Ankit

Blockquote

.

4

Projecting a linear function into a quadratic function space will still render the function linear in time.

Imagine that you have uh = x
Then up = x, and thus the double derivative is still zero.
This can be illustrated with the following:

import dolfinx
from mpi4py import MPI
import numpy as np
import ufl

cell = ufl.Cell("triangle", geometric_dimension=2)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, 1))
x = np.array([[0., 0.], [1., 0.], [0., 1.]])
cells = np.array([[0, 1, 2]], dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)

V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
W = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 2))

u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0]+3*x[1])

w = dolfinx.fem.Function(W)

uh = ufl.TrialFunction(W)
vh = ufl.TestFunction(W)
a = ufl.inner(uh, vh) * ufl.dx
L = ufl.inner(u, vh) * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, u=w)
problem.solve()

print(w.x.array)

when we get

[-2.22478286e-14  1.00000000e+00  3.00000000e+00  2.00000000e+00
  1.50000000e+00  5.00000000e-01]

which matches the ordering of

i.e. the function can be expressed as:

\begin{align} u &= \sum_{i=0}^{5} c_i \phi_i\\ & = 0 \cdot (2x^2+4xy-3x+2y^2-3y+1) \\ &+ 1 \cdot (x (2x-1)) + 3 \cdot ( y (2y-1))\\ &+ 2\cdot (4xy)\\ &+ 1.5 \cdot (4y (-x-y+1))\\ &+ 0.5 \cdot ( 4x(-x-y+1))\\ &= 0\\ &+2x^2-x + 6y^2 - 3y\\ &+8xy\\ &-6xy- 6y^2 + 6y\\ & - 2x^2 - 2xy +2x\\ &= x+3y \end{align}

and thus also the double derivative of this function will be zero.