Dear all,
I am interesting in the evaluation of an expression at certain points of a facet. More precisely, I want to evaluate it at the interpolation points in a Discontinuos Galerkin (DG) space.
To do so, I have followed the approach in Cohesive zone modeling of debonding and bulk fracture — Computational Mechanics Numerical Tours with FEniCSx, section Facet expressions interpolation. In that tutorial, a basix quadrature element is defined with integration points equal to the interpolation points, see:
\int_F eq \, dS = \sum_{g=1}^n |F| e \, \omega_g (x_g)q(x_g)
with x_g the integration points being equal to the interpolation points, \omega the weights, e the expression to be evaluated, F the facet and |F| the area measure, and q a test function.
Then, the form obtained after quadrature integration is divided by the area of the facet to end up with an assembled vector that contains exactly the values of the expression (and not their integration in the area). For a DG of degree 1, in that tutorial the weights are set to \omega_g=1.
My question is why the weights are set to 1 and not to \omega_g=0.5, since there is two integration points and not just one, as far as I understand. I might be missing some details about how the basix.ufl.quadrature_element works to justify the value of 1 for the weights.
Thank you in advance.
We use here a “trick” to form a vector containing the value of the expression to evaluate at quadrature points, we do not want to do a correct evaluation of the integral. By computing
\int_F \dfrac{e}{|F|}q dS
with \omega_g=1, we obtain for this form with test function q in the quadrature space:
\sum_{g=1}^{n}e(x_g)q(x_g)
Hence, we get a vector with values exactly e(x_g). This also works for a larger value of quadrature points than 2 (DG 2 space for instance).
Thank you for the response. Does any special consideration need to be taken for 2D faces rather than 1D ones?
I ask because I implemented this “trick” for 2D facets and, for a toy problem with an expression e equal to a constant value of e=1, the resulting value–in Paraview, interpolated into a DG, 1–was exactly doubled, i.e., 2, and not 1. Thus, I thought that the quadrature scheme in ufl may consider a quadrature rule based on an isoparametric element in the spirit of \xi\in[-1,1], and that the weights would be referred to such a reference element with Jacobian of the transformation of coordinates half the length of the Facet |F|.
However, from what you describe, I understand that the parameters of the quadrature rule defined inside the ufl.Measure directly stands for:
\sum_{g=1}^n |F| \omega_g e(x_g)q(x_g),
with |F| the total surface of the face (global domain). Is this right? If so, I will keep searching for the cause not in the quadrature rule.
Thank you.
In basix, the reference interval is [0;1], so yes in 2D, you have to adapt the formula to the reference triangle area which is 1/2. Hence, the factor 2.