Inconsistent derivatives at the quadrature points

dolfinx
1

Dear community, I found the inconsistent derivatives at the quadrature points (in Dolfinx v0.8.0) and would like to know how to fix it. Here is the MWE:

from mpi4py import MPI
from dolfinx.mesh import create_rectangle, CellType
import scipy
import ufl

from dolfinx.fem import functionspace, Function, Expression, form
from dolfinx.fem.petsc import create_vector, assemble_vector, create_matrix, assemble_matrix
from basix.ufl import quadrature_element


mesh = create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [2, 5]],
                        [1, 1], CellType.quadrilateral)
mesh.geometry.x[0, 1] += 0.5

S0 = functionspace(mesh, ("DG", 0))
S = functionspace(mesh, ("CG", 1))
quad_element = quadrature_element(mesh.topology.cell_name(), degree=2)

SQ = functionspace(mesh, quad_element)
u = Function(S)
uq = Function(SQ)
v = ufl.TestFunction(SQ)

# Compute element volumes
v0 = ufl.TestFunction(S0)
volume_form = form(v0*ufl.dx)
volume_vec = create_vector(volume_form)
assemble_vector(volume_vec, volume_form)
vol = Function(S0)
vol.vector.array = volume_vec.array

# FEniCSx solution
num_pts = 4
dudu_form = form(ufl.derivative(num_pts*u*v/vol*ufl.dx, u))
dudu_mat = create_matrix(dudu_form)
assemble_matrix(dudu_mat, dudu_form)
dudu_mat.assemble()
ai, aj, av = dudu_mat.getValuesCSR()
mat = scipy.sparse.csr_matrix((av, aj, ai)).todense()

# Finite difference
pts = uq.function_space.element.interpolation_points()
u_expr = Expression(u, pts)
uq.interpolate(u_expr)
value1 = uq.vector.array.copy()

idx, eps = 0, 1e-6
u.vector.array[idx] += eps
uq.interpolate(u_expr)
value2 = uq.vector.array.copy()
fenicsx_sol = mat[:, idx].ravel()
finite_difference = (value2-value1) / eps

print(fenicsx_sol)
print(finite_difference)

The results are

[[0.60310759 0.17173114 0.16160219 0.04601522]]
[0.62200847 0.16666667 0.16666667 0.0446582 ]

and we can see the inconsistencies.

If we use the degree-1 quadrature element or use a rectangular element by commenting

mesh.geometry.x[0, 1] += 0.5

the derivatives are consistent again:

[[0.62200847 0.16666667 0.16666667 0.0446582 ]]
[0.62200847 0.16666667 0.16666667 0.0446582 ]
2

It would help if you did specify what the two methods are trying to compute (mathematically), and why they should be equivalent.

As far as i can tell by looking at this on my phone:

makes the element non-affine, as you are perturbing a single coordinate in your geometry.

Then you might be underintegrating the form:

(Note that it is unclear why you call the derivative here, as the form is linear in u.)

3

Thanks for the reply, @dokken! My initial goal is to compute a Jacobian matrix, \partial p_i / \partial q_j, with p_i defined at the quadrature point (such as the strain component) and q_j defined at the node (such as the displacement). Indices i and j are the DOF counts of p and q, respectively.

In the current simplified setup, I set p = u defined at the quadrature point and q = u defined at the node. The first method is something artificial. It is less explainable and may not be correct, although it seems to work for rectangular elements. The second method is the finite difference to approximate \partial p_i / \partial q_j.

Currently, I am working on the unstructured mesh, and the element might be non-affine. So I wonder if there is a way to evaluate \partial p_i / \partial q_j for the unstructured mesh.

4

I have fixed this problem now and would like to share the revised code here for future readers:

from mpi4py import MPI
from dolfinx.mesh import create_rectangle, CellType
import scipy
import ufl
import numpy as np

from dolfinx.fem import functionspace, Function, Expression, form
from dolfinx.fem.petsc import create_matrix, assemble_matrix
from basix.ufl import quadrature_element
import basix


mesh = create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [3.5, 4]],
                        [1, 2], CellType.quadrilateral)
mesh.geometry.x[0, 1] += 0.2

quad_degree = 1
quad_element = quadrature_element(mesh.topology.cell_name(), degree=quad_degree)

S0 = functionspace(mesh, ("DG", 0))
S = functionspace(mesh, ("CG", 1))
SQ = functionspace(mesh, quad_element)

u = Function(S)
w = Function(SQ)
uq = Function(SQ)
v = ufl.TestFunction(SQ)

quad_points, quad_weights = basix.make_quadrature(
    basix.cell.string_to_type(mesh.topology.cell_name()), degree=quad_degree)
w.vector.array = np.tile(quad_weights, int(w.vector.array.size/quad_weights.size))

# Direct computation
factor = w * abs(ufl.JacobianDeterminant(mesh))
expr = ufl.derivative(u*v/factor*ufl.dx, u)
dudu_form = form(expr)
dudu_mat = create_matrix(dudu_form)
assemble_matrix(dudu_mat, dudu_form)
dudu_mat.assemble()
ai, aj, av = dudu_mat.getValuesCSR()
mat = scipy.sparse.csr_matrix((av, aj, ai)).todense()

# Finite difference
u_expr = Expression(u, quad_points)
uq.interpolate(u_expr)
value1 = uq.vector.array.copy()

idx, eps = 0, 1e-6
u.vector.array[idx] += eps
uq.interpolate(u_expr)
value2 = uq.vector.array.copy()

# Check the values
fenicsx_sol = mat[:, idx].ravel()
finite_difference = (value2-value1) / eps
print(fenicsx_sol)
print(finite_difference)

The key line is factor = w * abs(ufl.JacobianDeterminant(mesh)). For a variable, p, defined at a quadrature point, it is evaluated as m = w \cdot j \cdot p in FEniCSx. Here, w is the weight at that quadrature point, and j is the Jacobian determinant.

To get the p value, we simply evaluate p = m / f with f = w \cdot j. In my original post, it accidentally holds N = 1/w and V \equiv j for an affine quad element. Here, N is the number of integration points per element, and V is the element volume. For the non-affine element, we can use the code above.

1 Like