Assemble forms with test functions

dolfinx
,
1

Hello,

I’m trying to compute the integral of fv, where f is a function defined in a finite element function space V, and v is the test function of the same space.

In Python, this is very straightforward and works correctly for any degree I specify. Here is a minimal working example:

from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
from ufl import dx

degree = 3
msh = mesh.create_interval(comm=MPI.COMM_WORLD, nx=10, points=(0.0, 1.0))
V = fem.functionspace(msh, ("Lagrange", degree))

v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2) / 0.01)
L = fem.form(f * v * dx)

print(fem.assemble_scalar(L))

No matter what degree I choose, this implementation runs smoothly.

However, when I attempt a similar approach using the C++ interface, I encounter issues. I construct the form using UFL and Basix as follows:

from basix.ufl import element
from basix import LagrangeVariant
from ufl import *

cell = "interval"
degree_cg = 3
e = element("Lagrange", cell, degree_cg, LagrangeVariant.equispaced)

coord = element("Lagrange", cell, 1, shape=(1,))
mesh = Mesh(coord)

V = FunctionSpace(mesh, e)

v = TestFunction(V)
f = Coefficient(V)
a = f * v * dx

Then, in C++, I create and assemble the form as:

auto a = std::make_shared<Form<T,U>>(fem::create_form<T,U>(*form_test_1D_P3_gauss1, {V_}, {{"f", Ff_}}, {}, {}));
T g = dolfinx::fem::assemble_scalar(*a);

This works fine for polynomial degrees 1 and 2. However, when I use degree 3, I encounter the following error:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

It seems like the issue is specific to higher-order elements, possibly related to quadrature or DOF mapping. I’d appreciate any insights into what might be causing this or suggestions for resolving it.

2

This may show up my utter ignorance… but I’m actually puzzled as to why your original example assembles a scalar value

from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
from ufl import dx

degree = 3
msh = mesh.create_interval(comm=MPI.COMM_WORLD, nx=10, points=(0.0, 1.0))
V = fem.functionspace(msh, ("Lagrange", degree))

v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2) / 0.01)
L = fem.form(f * v * dx)

print(fem.assemble_scalar(L))

The linear formulation in L should be defined as the vector

\int_{k} f \, v \, \mathrm{d}\vec{x} \quad \forall k \in \mathcal{T}^h(\Omega)

where \mathcal{T}^h is the mesh of your domain \Omega and k is an element in the mesh.

I suspect DOLFINx is not performing an arity check in your simulation. I’d be careful as I think the “scalar quantity” you’re computing is meaningless.

1 Like
3

I’m sorry I didn’t clearly state my purpose.
Actually, I’m checking the error of the Gauss’s law, for example, in 1D it should be:
\partial_xE = \rho
because I solved it in the weak form, so the above conclusion also only hold in the weak form:
-(E_h,\partial_x\phi) = (\rho_h,\phi),\qquad\forall\phi\in V
And what I wanted to do is compute the values of the two terms on the two sides of equation, and compare the difference, that’s why I have the test function in the form.

Actually, it used to work well, but recently I updated my DOLFINx from version v0.6.0 to v.7.3, then the problems happen. I have attached my codes before:

import basix
from basix.ufl_wrapper import create_element, create_vector_element
from ufl import Coefficient, TestFunction, TrialFunction, dx, CellVolume, Mesh, div, grad, dot

degree_cg = 3
cell = "interval"

element = create_element("Lagrange",cell, degree_cg, basix.LagrangeVariant.equispaced)
coord_element = create_vector_element("DG", cell, 3)
mesh = Mesh(coord_element)

v = TestFunction(element)
rho = Coefficient(element)
gauss = -rho * v * dx

and

auto gauss = std::make_shared<Form<T>>(fem::create_form<T>(*form_Maxwell_1D_Q3_gauss, {V}, {{"rho", rho_}}, {}, {}));

I can see there are difference between, but not sure how can I change it, since I followed the demos.

4

Acutally the scalar quantity is just what I need, since I want to compute the error, but I don’t know why it doesn’t work for P3 elements.

5

The problem is that the equation above does not yield a scalar quantity.
As you need to test for all \phi_i\in V_h, you get a vector b_i with entries
\begin{pmatrix}(\rho, \phi_0)\\ (\rho, \phi_1)\\ \vdots \\(\rho, \phi_{N-1})\end{pmatrix}

When you mean computing the error between the two, you should compute the two vectors, and then compare them element by element, or equivalently assemble the difference into a single vector, and then perform an l_2 norm.

2 Likes
6

Ok I see, yes you are right, the test function is not a single function, so I should get a vector instead of a scalar.