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 Teach FiniteElement how to fill MF::ShapeInfo #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.
• Teach the matrix-free evaluators the case with anisotropic polynomials.
• 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

  dealii/include/deal.II/matrix_free/mapping_info.templates.h

  Lines 86 to 88 in 7dd7f14

  	this->update_flags_cells =
  	MappingInfoStorage<dim, dim, VectorizedArrayType>::compute_update_flags(
  	update_flags_cells, quad);

  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

  dealii/source/fe/fe_raviart_thomas_nodal.cc

  Lines 528 to 572 in f541dd4

  	const unsigned int n = this->degree;
  	const unsigned int face_no = 0;
  	Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
  	for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
  	// face support points are in lexicographic ordering with x running
  	// fastest. invert that (y running fastest)
  	{
  	unsigned int i = local % n, j = local / n;
  	// face_orientation=false, face_flip=false, face_rotation=false
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	0) =
  	j + i * n - local;
  	// face_orientation=false, face_flip=false, face_rotation=true
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	1) =
  	i + (n - 1 - j) * n - local;
  	// face_orientation=false, face_flip=true, face_rotation=false
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	2) =
  	(n - 1 - j) + (n - 1 - i) * n - local;
  	// face_orientation=false, face_flip=true, face_rotation=true
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	3) =
  	(n - 1 - i) + j * n - local;
  	// face_orientation=true, face_flip=false, face_rotation=false
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	4) = 0;
  	// face_orientation=true, face_flip=false, face_rotation=true
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	5) =
  	j + (n - 1 - i) * n - local;
  	// face_orientation=true, face_flip=true, face_rotation=false
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	6) =
  	(n - 1 - i) + (n - 1 - j) * n - local;
  	// face_orientation=true, face_flip=true, face_rotation=true
  	this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
  	7) =
  	(n - 1 - j) + i * n - local;
  	// for face_orientation == false, we need to switch the sign
  	for (unsigned int i = 0; i < 4; ++i)
  	this->adjust_quad_dof_sign_for_face_orientation_table[face_no](local,
  	i) = 1;

  (in particular quad_dof_sign_for_face_orienation_table) and then apply the operations from

  dealii/source/fe/fe_poly_tensor.cc

  Lines 531 to 537 in f541dd4

  	if (adjust_quad_dof_sign_for_face_orientation(
  	local_quad_dof_index,
  	face_index_from_dof_index,
  	cell->face_orientation(face_index_from_dof_index),
  	cell->face_flip(face_index_from_dof_index),
  	cell->face_rotation(face_index_from_dof_index)))
  	dof_sign = -1.0;

  and

  dealii/source/fe/fe_poly_tensor.cc

  Lines 595 to 598 in f541dd4

  	for (unsigned int k = 0; k < n_q_points; ++k)
  	for (unsigned int d = 0; d < dim; ++d)
  	output_data.shape_values(first + d, k) =
  	dof_sign * fe_data.transformed_shape_values[k][d];

  .
• 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

  dealii/source/fe/fe_raviart_thomas_nodal.cc

  Lines 41 to 42 in f541dd4

  	template <int dim>
  	class PolynomialsRaviartThomasNodal : public TensorPolynomialsBase<dim>

  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):

• For the evaluation of RT polynomials in evaluation_kernels.h, one has to (a) branch here

  dealii/include/deal.II/matrix_free/evaluation_kernels.h

  Lines 1772 to 1775 in f541dd4

  	run(const unsigned int n_components,
  	const EvaluationFlags::EvaluationFlags evaluation_flag,
  	const Number * values_dofs,
  	FEEvaluationData<dim, Number, false> & fe_eval)

  as another case with another ElementType. This needs to be made in sync with step 2 that fills ShapeInfo
• Select the right anisotropic polynomials at FEEvaluationImpl

  dealii/include/deal.II/matrix_free/evaluation_kernels.h

  Lines 211 to 213 in f541dd4

  	Eval eval0 = create_evaluator_tensor_product(univariate_shape_data[0]);
  	Eval eval1 = create_evaluator_tensor_product(univariate_shape_data[1]);
  	Eval eval2 = create_evaluator_tensor_product(univariate_shape_data[2]);

• Implement the actual interpolation with anisotropic polynomials and the corresponding offsets with a new class analogous to

  dealii/include/deal.II/matrix_free/tensor_product_kernels.h

  Lines 1764 to 1771 in f541dd4

  	EvaluatorTensorProduct<evaluate_evenodd,
  	dim,
  	n_rows,
  	n_columns,
  	Number,
  	Number2>::apply(const Number2 *DEAL_II_RESTRICT shapes,
  	const Number * in,
  	Number * out)

• 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

  dealii/include/deal.II/matrix_free/tensor_product_kernels.h

  Lines 1801 to 1804 in 7dd7f14

  	for (int i2 = 0; i2 < n_blocks2; ++i2)
  	{
  	for (int i1 = 0; i1 < n_blocks1; ++i1)
  	{

  or

  dealii/include/deal.II/matrix_free/tensor_product_kernels.h

  Lines 1441 to 1444 in 7dd7f14

  	for (int i2 = 0; i2 < n_blocks2; ++i2)
  	{
  	for (int i1 = 0; i1 < n_blocks1; ++i1)
  	{

  or

  dealii/include/deal.II/matrix_free/tensor_product_kernels.h

  Lines 360 to 363 in 7dd7f14

  	for (int i2 = 0; i2 < n_blocks2; ++i2)
  	{
  	for (int i1 = 0; i1 < n_blocks1; ++i1)
  	{

  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 #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.