Dear FreeFem++ community,
I am using H(div) finite elements for Maxwell’s equations and I have a few questions regarding the accuracy of the implementations. For some cases, I have an analytical expression of the vector field from which to calculate the degrees of freedom. When I compare the interpolated values to the analytical expression, I find that the divergence is off by a factor of around 3 when using RT0 finite elements. When using BDM1 finite elements, they agree quite well. Now, the documentation states in regard to RT0 finite elements:
“To compute the flux, we use a quadrature with one Gauss point, the mid-point of the edge.”
Nothing is said on BDM1 elements, but I assume two quadrature points are used since there are also two DoFs per edge. Since I use FreeFem++ in conjunction with Fortran, I tried to calculate the RT0 DoFs with two (or more) quadrature points in Fortran and pass the values to FreeFem++. This way, the resultant value for the divergence is on par with the value obtained from BDM1 elements. This raises a few questions:
Is there an option to specify the quadrature formula to be used when projecting analytical expressions to the FE function? As far as I understand it, this is already possible for int2d etc.
Could you provide the formulas used to calculate these projections, and the derivatives (dx, dy etc.)? This could be included in the documentation but I would already be happy with a reference to the relevant C++ sources. This would also be interesting for the RT1 finite elements; the results I get suggest some other kind of bug.