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.
Best regards,
Patrick Lainer