Teach FEEvaluation RT and Nedelec evaluation

[kronbichler]

This issue is to keep track of the work to support Raviart-Thomas and Nedelec elements in the context of MatrixFree and FEEvaluation. We need the following steps:

• [ ] We need suitable polynomial spaces compatible with tensor products and ideally with matrix-free evaluators. I am thinking of nodal variants in terms of Gauss/Gauss-Lobatto bases or, for the case we want to do face integrals (incompressible flow with velocity represented by RT), we would probably even prefer Polynomials::HermiteLikeInterpolation. I have yet to review what polynomials we have exactly in the various RT and Nedelec classes; at least at first sight they do not seem to be easily decomposable into tensor products.
• [ ] Teach internal::MatrixFreeFunctions::ShapeInfo to read out the information for RT and Nedelec elements. Once we have #9544 in place, it should be straight-forward to do so; one extension we still need is support for anisotropic polynomials, i.e., the tensor product of multiple 1D formulas because we have degree k in some directions and k+1 in some other.
• [x] Teach the matrix-free evaluators the case with anisotropic polynomials.
• [x] Implement the Piola and Nedelec mappings in FEEvaluation. Given that this is a switch to a different formula in all get_value(), submit_value(), get_curl functions etc., we could either add it by an additional template parameter MappingKind that selects the right formula, or by a simple if statement. The latter would be less intrusive but require additional if clauses in inner loops. The compiler is typically pretty good at figuring them out, but if we put too many it might at some point refuse to optimize what we want to be optimized. We can start with that option. Note that for gradients with RT we will need the gradient of the Jacobian of the mapping, so the places where we construct the mapping in https://github.com/dealii/dealii/blob/7dd7f14b46ab44e2c460e1fbcf159147e463e91a/include/deal.II/matrix_free/mapping_info.templates.h#L86-L88 need to get aware of this fact.
• [ ] Regarding the case of faces/edges not in the standard orientation: We need to adjust the basis functions so that Piola transform is still correct. Probably we need https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/source/fe/fe_raviart_thomas_nodal.cc#L528-L572 (in particular quad_dof_sign_for_face_orienation_table) and then apply the operations from https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/source/fe/fe_poly_tensor.cc#L531-L537 and https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/source/fe/fe_poly_tensor.cc#L595-L598.
• [ ] Check what else needs to be done with adaptive grids. I'm undecided if we need to introduce a factor of 0.5 or 2. somewhere, so that the Piola transformation looks the same on both sides.

[jwitte08]

With a student, I am working on RT elements build up as a tensor product of 1D elements.

Today, we had a look into FE_RaviartThomas and FE_RaviartThomasNodal's implementations. The first is close to the mathematical theory and assumes less regularity from the ansatz space. However, the latter (based on a nodal, primitive ansatz) has clear advantages seen from a technical perspective.

For example, in 2D one constructs an anisotropic quadrature formula Q_{k+2,k+1} for the x-component and Q_{k+1,k+2} for the y-component using QProjector. The node functionals with a one-to-one correspondence to these quadrature points uniquely determine the continuity in normal direction. Therefore, each component is readily decomposed into 1D node functionals.

Thus, we decided to start with the FE_RaviartThomasNodal filling the ShapeInfo having #9544 in mind. Let me know if you need some help with the new interface provided by FiniteElement.

[bangerth]

On 2/20/20 2:05 PM, jwitte08 wrote:

Today, we had a look into |FE_RaviartThomas| and |FE_RaviartThomasNodal|'s implementations. The first is close to the mathematical theory and assumes less regularity from the ansatz space. However, the latter (based on a nodal, primitive ansatz) has clear advantages seen from a technical perspective.

If one reads through old code, there are always places where you now know that things could be written in a cleaner way, differently, or with better comments. Feel free to make these changes and submit pull requests for them -- code cleanups are always welcome!

[kronbichler]

After discussion with @nfehn, we came up with a more detailed plan for approaching the problem.

Regarding the first point:

• [ ] Extend the PolynomialsRaviartThomas class https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/source/fe/fe_raviart_thomas_nodal.cc#L41-L42 so that also the d * (k+1)^d space can be represented besides the RT - this is not strictly needed for RT, but makes it easier to create other H(div) conforming spaces and just test out non-anisotropic spaces can be done. I have an idea how to make that work.

Regarding the third bullet point (make the evaluators capable of doing anisotropic polynomials):

• [x] For the evaluation of RT polynomials in evaluation_kernels.h, one has to (a) branch here https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/include/deal.II/matrix_free/evaluation_kernels.h#L1772-L1775 as another case with another ElementType. This needs to be made in sync with step 2 that fills ShapeInfo
• [x] Select the right anisotropic polynomials at FEEvaluationImpl https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/include/deal.II/matrix_free/evaluation_kernels.h#L211-L213
• [x] Implement the actual interpolation with anisotropic polynomials and the corresponding offsets with a new class analogous to https://github.com/dealii/dealii/blob/f541dd4710b82663a96a59ae395609792a91c52d/include/deal.II/matrix_free/tensor_product_kernels.h#L1764-L1771
• [x] It might be that the previous point is more easily addressed by breaking out the 1D interpolation routines from the loop over the other coordinate dimensions like https://github.com/dealii/dealii/blob/7dd7f14b46ab44e2c460e1fbcf159147e463e91a/include/deal.II/matrix_free/tensor_product_kernels.h#L1801-L1804 or https://github.com/dealii/dealii/blob/7dd7f14b46ab44e2c460e1fbcf159147e463e91a/include/deal.II/matrix_free/tensor_product_kernels.h#L1441-L1444 or https://github.com/dealii/dealii/blob/7dd7f14b46ab44e2c460e1fbcf159147e463e91a/include/deal.II/matrix_free/tensor_product_kernels.h#L360-L363 that are repetitive anyway.

Regarding the fourth and fifth point I made some more specific comments in the original post.

[kronbichler]

Apart from the cleanup in https://github.com/dealii/dealii/issues/12728, I think we need to postpone the additional topics in this issue to after the 9.4 release. I think we have a nice "experimental" state in terms of affine transformations and basic stuff, and can fix hanging nodes (#13671) as well as non-affine meshes to the next release. We would not have the time to properly test and implement all missing features to have full support for RT (let alone Nedelec, that was also part of the topic).

[kronbichler]

We made progress with #13990 and some preceding commits, but the rest needs to be postponed.