Teach FiniteElement how to fill MF::ShapeInfo

[kronbichler]

The information for the tensor product information of certain FiniteElement classes is extracted in the struct internal::MatrixFreeFunctions::ShapeInfo from a given finite element object and a 1D quadrature formula. We do this by some dynamic cast magic in

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

Lines 103 to 128 in 1f3f4df

	const FE_Poly<TensorProductPolynomials<dim>, dim, dim> *fe_poly =
	dynamic_cast<
	const FE_Poly<TensorProductPolynomials<dim>, dim, dim> *>(fe);
	const FE_Poly<
	TensorProductPolynomials<dim,
	Polynomials::PiecewisePolynomial<double>>,
	dim,
	dim> *fe_poly_piece =
	dynamic_cast<const FE_Poly<
	TensorProductPolynomials<dim,
	Polynomials::PiecewisePolynomial<double>>,
	dim,
	dim> *>(fe);
	const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim> *>(fe);
	const FE_Q_DG0<dim> *fe_q_dg0 = dynamic_cast<const FE_Q_DG0<dim> *>(fe);
	element_type = tensor_general;
	if (fe_poly != nullptr)
	scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
	else if (fe_poly_piece != nullptr)
	scalar_lexicographic =
	fe_poly_piece->get_poly_space_numbering_inverse();
	else if (fe_dgp != nullptr)

but this is not a real solution and does not scale if we want to support more elements (and I want to add some tensor product evaluators for RT and Nedelec with appropriate basis functions).
The real solution is to add this as an interface to the FiniteElement classes and let them return tensor product information. I plan to do this in three steps:

• Move the include/deal.II/matrix_free/shape_info.{h,templates.h} files to the fe/ subfolder, extract it from the internal namespace and give it a proper name. I was thinking along the lines FiniteElementTensorProductInformation. I would put a template for the number type on this class but not the dimension as a tensor product element is in principle not restricted to a particular dimension.
• Create a virtual function FiniteElement<dim,spacedim>::extract_tensor_product_information(const Quadrature<1> &quadrature, const unsigned int base_element_index = 0) that returns a FiniteElementTensorProductInformation<double>. For FESystem, the data is filled partly by the children.
• Use this information in MatrixFree by copying to the chosen Number type there (e.g. VectorizedArray<float>) and create an alias at the existing ShapeInfo place.

What I want to discuss here before I start is the naming (any suggestions?) and how exactly we would like to set up the interface. One can think of this class/struct as a collection of arrays with evaluation of 1D shape functions along a 1D quadrature and for various combinations of faces/subfaces, so it is pretty low-level and a bit similar to the InternalDataBase objects in FiniteElement, but without the inheritance within the class.

[peterrum]

This is a good idea!

I guess then we can also eliminate this mess:

dealii/source/multigrid/mg_transfer_internal.cc

Lines 713 to 722 in 1f3f4df

	std::string fe_name = dof_handler.get_fe().base_element(0).get_name();
	{
	const std::size_t template_starts = fe_name.find_first_of('<');
	Assert(fe_name[template_starts + 1] ==
	(dim == 1 ? '1' : (dim == 2 ? '2' : '3')),
	ExcInternalError());
	fe_name[template_starts + 1] = '1';
	}
	const std::unique_ptr<FiniteElement<1>> fe(
	FETools::get_fe_by_name<1, 1>(fe_name));

and maybe part of this:

dealii/source/multigrid/mg_transfer_internal.cc

Lines 582 to 640 in 1f3f4df

	template <int dim, typename Number>
	void
	setup_element_info(ElementInfo<Number> & elem_info,
	const FiniteElement<1> & fe,
	const dealii::DoFHandler<dim> &dof_handler)
	{
	// currently, we have only FE_Q and FE_DGQ type elements implemented
	elem_info.n_components = dof_handler.get_fe().element_multiplicity(0);
	AssertDimension(Utilities::fixed_power<dim>(fe.dofs_per_cell) *
	elem_info.n_components,
	dof_handler.get_fe().dofs_per_cell);
	AssertDimension(fe.degree, dof_handler.get_fe().degree);
	elem_info.fe_degree = fe.degree;
	elem_info.element_is_continuous = fe.dofs_per_vertex > 0;
	Assert(fe.dofs_per_vertex < 2, ExcNotImplemented());
	// step 1.2: get renumbering of 1D basis functions to lexicographic
	// numbers. The distinction according to fe.dofs_per_vertex is to support
	// both continuous and discontinuous bases.
	std::vector<unsigned int> renumbering(fe.dofs_per_cell);
	{
	AssertIndexRange(fe.dofs_per_vertex, 2);
	renumbering[0] = 0;
	for (unsigned int i = 0; i < fe.dofs_per_line; ++i)
	renumbering[i + fe.dofs_per_vertex] =
	GeometryInfo<1>::vertices_per_cell * fe.dofs_per_vertex + i;
	if (fe.dofs_per_vertex > 0)
	renumbering[fe.dofs_per_cell - fe.dofs_per_vertex] =
	fe.dofs_per_vertex;
	}
	// step 1.3: create a dummy 1D quadrature formula to extract the
	// lexicographic numbering for the elements
	Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
	ExcNotImplemented());
	const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
	const unsigned int n_child_dofs_1d =
	(fe.dofs_per_vertex > 0 ? (2 * fe.dofs_per_cell - 1) :
	(2 * fe.dofs_per_cell));
	elem_info.n_child_cell_dofs =
	elem_info.n_components * Utilities::fixed_power<dim>(n_child_dofs_1d);
	const Quadrature<1> dummy_quadrature(
	std::vector<Point<1>>(1, Point<1>()));
	internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
	shape_info.reinit(dummy_quadrature, dof_handler.get_fe(), 0);
	elem_info.lexicographic_numbering = shape_info.lexicographic_numbering;
	// step 1.4: get the 1d prolongation matrix and combine from both children
	elem_info.prolongation_matrix_1d.resize(fe.dofs_per_cell *
	n_child_dofs_1d);
	for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
	for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
	for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
	elem_info
	.prolongation_matrix_1d[i * n_child_dofs_1d + j + c * shift] =
	fe.get_prolongation_matrix(c)(renumbering[j], renumbering[i]);
	}

[peterrum]

Move the include/deal.II/matrix_free/shape_info.{h,templates.h} files to the fe/ subfolder, extract it from the internal namespace and give it a proper name. I was thinking along the lines FiniteElementTensorProductInformation. I would put a template for the number type on this class but not the dimension as a tensor product element is in principle not restricted to a particular dimension.

I just wanted to suggest the same!

[kronbichler]

I guess then we can also eliminate this mess: [some ugly code in the multigrid transfer]

I hope that we can do; however, this means we need to add another field to the ShapeInfo class for the parent/child relation we use in GMG (which is useful information and cheap to obtain anyway).

[peterrum]

@kronbichler Let me know if I can do something or try something out. I have simply copied those two code snippets to my new transfer operators, but I am not happy with that...

I will have a look at my code maybe I find some additional requirements!

[kronbichler]

The evaluator side also relates to #7039 that should be generalized/restructured in one go.

[peterrum]

This PR will involve some time we don't have at the time. Adjusting the milestone.

[kronbichler]

In #14385 we noticed that the cases where we need special treatments, to be addressed by this PR, just keep increasing. I think we should aim for this effort before the next release.

[peterrum]

This is a place where this would be useful:

dealii/source/dofs/dof_accessor_set.cc

Lines 253 to 254 in e8f5713

	fe.get_prolongation_matrix(child, cell.refinement_case())
	.vmult(tmp, local_values);

since we could replace the full matrix-vector product by sum factorization.