Enhance VectorTools::project

[ghost]

Since there is now matrix-free+hp support, the VectorTools::project functions might be updated in order to increase the capabilities. As I use a similar function to project the gradients of a given finite element function I tried to revise the current implementation of VectorTools::project. In the following a short list of changes:

• Since we pass in VectorTools::internal::project_matrix_free() the finite element degree as -1 to the MatrixFreeOperators::MassOperator the documentation, which says 'the degree of FiniteElement is less than nine', can be removed.
• Personally I consider multi-physical problems, where I get immediately more than four components. Thus I increased the number of supported components to ten and updated the documentation accordingly.
• Use range-based loops for the loop over all cells or all locally_owned_elements
• Adapt the matrix-free functions to handle the hp case as well

I have also some questions:

• In the function VectorTools::internal::project we check if we could use the matrix-free implementation by asking is_supported() and n_base_elements()==1. Looking into deal.II/matrix_free/shape_info.templates.h I assume that the lines

const FiniteElement<dim, spacedim> *fe_ptr = &(fe.base_element(base));
if (fe_ptr->n_components() != 1)
  return false;

just ask

dof.get_fe(0).n_base_elements() == 1

Am I wrong?

• Is there a more elegant way to copy the passed quadrature to hp::QCollection<dim> quadrature_mf in project_matrix_free()?
• I am coupling different domains with different physics like it is proposed in step-46. Hence I use the hp capabilities together with the FE_Nothing element. As I see it correctly in deal.II/matrix_free/shape_info.templates.h, this element is not supported and thus the project function would switch to the generic matrix-based implementation. However, then it is not possible to use a p::d::Tria? Wouldn't this be a desirable feature?
• In the current state there is no implementation for the parallel generic matrix-based case? So if we use a p::d::Tria and more than four components an error is thrown?

And some possible remaining TODOs:

• Maybe some tests would be help full
• Changelog entry

[masterleinad]

Overall, this looks pretty good to me. Some remarks:

• We clearly need some tests.
• We deliberately decide to only support a small number of components since compiling all the template instantiations that are pretty expensive for the matrix-free case. We should critically evaluate how much more expensive supporting up to 10 components is. In the worst case, users should be able to call the internal functions using the matrix-free framework directly if a larger number of components is required.
• We decided to only support parallel Triangulations for the matrix-free case since the support for matrix-based approaches, in this case, would require external packages and the implementation must depend on the configurations provided. We could certainly do that but there doesn't seem to have been a large demand (also because the current implementation should cover the most common use cases). This should be a separate discussion IMHO.

[ghost]

Thanks for your remarks!

• I thought about the tests. Since there was hp support previously I would assume that there are already some serial hp tests. The new tests then sould just cover the parallel hp case, right? Is it straight foreward extending them?

• OK, the evaluation might exceed my knowledge of the dealii library. I can live with four components in the original dealii and work with my fork. For sake of simplicity I revert the changes and go back to four components.

• OK, fine, I can live with that, just asked out of interest.

Unfortunately, the library is not yet compiling due to some undefined references to VectorTools::create_right_hand_side. As I see it correctly the one inst for create_right_hand_side in the hp + LA::d::Vector case is missing. In order to save compile time I have explicitly defined only the missing inst. Just an additionaly question: Don't we need a compress() call in the create_right_hand_side functions for the parallel cases?

[kronbichler]

We deliberately decide to only support a small number of components since compiling all the template instantiations that are pretty expensive for the matrix-free case. We should critically evaluate how much more expensive supporting up to 10 components is. In the worst case, users should be able to call the internal functions using the matrix-free framework directly if a larger number of components is required.

Just to comment on this aspect as I don't have time at the moment to go deeper into the other comments: I believe it would be good if could create a variant of our matrix-free loops (and possibly operator classes) that can handle arbitrary numbers of components by concatenating scalar evaluations one after the other. Our class FEEvaluation supports to hook into an arbitrary vector component, as long as we use an FESystem. What we therefore need to do is to create a loop over the number of desired components in the operator inside the cell_operation (the one that gets called within MatrixFree::cell_loop or MatrixFree::loop) and selects one component after the other. We would need to check how we can integrate this into our operators, maybe by a template value -1 in the number of components, and adapt the code accordingly. @gfcas would you be interested in such a variant?

[ghost]

@kronbichler In its generality I would like such a variant as for the fe_degree, but I can also live with applying my changes locally. For me an more interesting issue would be (if it is one) allowing FE_Nothing in the parallel use_matrix_free==true case.

[ghost]

* We clearly need some tests.

@masterleinad If I see it correctly we can maybe just adapt the test from the simplex-case?

[masterleinad]

/rebuild

[masterleinad]

It should be straightforward to create MPI parallel versions of

• tests/hp/project_01.cc
• tests/hp/project_01_curved_boundary.cc
• tests/hp/project_02.cc

[marcfehling]

One of the simplex tests triggered an assertion:

The following tests FAILED:
	 11 - simplex/data_out_write_vtk_02.debug (Failed)

An error occurred in line <1910> of file </home/runner/work/dealii/dealii/include/deal.II/lac/vector_operations_internal.h> in function
    static Number dealii::internal::VectorOperations::functions<Number, Number2, dealii::MemorySpace::Host>::dot(const std::shared_ptr<dealii::parallel::internal::TBBPartitioner>&, dealii::internal::VectorOperations::size_type, const dealii::MemorySpace::MemorySpaceData<Number2, dealii::MemorySpace::Host>&, dealii::MemorySpace::MemorySpaceData<Number, dealii::MemorySpace::Host>&) [with Number = double; Number2 = double; dealii::internal::VectorOperations::size_type = unsigned int]
The violated condition was: 
    dealii::numbers::is_finite(sum)

Can you investigate?

The jenkins CI ran out of virtual memory (virtual memory exhausted). The next iteration of checks should run through.

[ghost]

Of course I will try to find it out. However, I am out of office for the next few week.

[bangerth]

Any news?

[ghost]

Any news?

Sorry! I haven't had time for it in the last few weeks because I still have some deadlines to meet. But I would like to improve that.

[peterrum]

@gfcas ?

[ghost]

@gfcas ?

I know this is still open, but I don't have time for it currently, unfortunately. I would still like to contribute.
I'll see what I can do.

[ghost]

This PR is outdated. I will try to review the current version of VectorTools and contribute maybe some single elements of this PR in a new one.