Question about using curved surface elements with a complex HCurl discretisation

Hi All,

I have a quick question concerning curved surface elements obtained via mesh.Curve(curve) for a complex HCurl discretisation. We want to repeatably compute integrals of the form $\int_\Omega E \cdot \overline{E} d \Omega$ where $E$ is known solution to a vectorial PDE. To do this, we have been computing x @ M @ conj(x), where x is the vector of solution coefficients for a grid function Ehp approximating E and M is the mass matrix.

To fix ideas, we consider a simple BVP Find $E$ s.t.

$\nabla \times \nabla \times E + E = f$ in $\Omega$
$ n \times E = n \times E_0$  on $\partial \Omega$

Where $E_0 = (sin(x+y) * (2+1j), -sin(x+y)*(2+1j), 0)$ (and is also the exact solution to the problem) and $f$ is a known source term. In this example, we solve using HCurl conforming order $p=0,1,2,...$ elements over a unit radius sphere $\Omega$. As the geometry is curved, we consider different mesh.Curve(curve), $curve\ge 1$. Recall, for p=0, the H(curl) elements approx the tangential component of E as a constant on edges, but the basis functions are linear.

First, we have tried Intval=Integrate(InnerProduct(Ehp,Ehp), mesh, order=2*(p+1)) Where Ehp is a grid function and we fix the number of integration points based on the degree of the integrand. Our understanding is that we can use Strang's second lemma to fix the number integration points independent of the order of the geometry approximation. We specify order as the default value for the max polynomial degree of the integrand is and we might want to consider higher order elements and also to make comparisons with the following.

Second, we set up a bilinear form for the inner product of a test and trial function and export the resulting mass matrix as a sparse scipy matrix M. Then compute Scival = x @ M @ conj(x)

We have found the following

1. If we use curve=1 (ie linear geometry, Intval and Scival agree to a high level of precision circa 1e-14 for tetrahedral and boundary layer meshes.
2. If we use curve > 1 and a boundary layer mesh (with prisms close to the sphere's surface) Intval and Scival agree to a high level of precision circa 1e-14
3. But if we use curve>1 and just a tetrahedral mesh we see differences in Intval and Scival. e.g. for curve=3 and p=1 of the order of 1e-6.

We think that the number of integration points used to compute the exported mass matrix in 3) for the BilinearForm may not be correct. Our understanding is the correct number of points should be obtained (automatically) from using integration points based on a polynomial of order 2*(p+1) (by applying Strang's lemma), but possibly 2*p for p\ge 1 is being used in case 3) only. Indeed, if we specify bonus_intorder=2 for the BilinearForm, this appears to fix the issues.

We think there is also a similar issue for LinearForm in 3) and a similar fix appears to work.

Note that in 1) and 2) there are no issues and so think the correct number of integration points are being automatically assigned for the Bilinear forms in these cases.

github.com/JamesElgy/NGSolve_HCurl_Example_Problem

The link above has a short Python script that recreates this problem using Python 3.10.6, NGSolve 6.2.2204 installed via the pip installer, and Ubuntu 22.04.

Does this sound correct? Is there indeed a possible issue with the number of integration points and exported mass matrix in use 3)?


Regards,
James and Paul