MatrixFree & vector valued FE : block-wise structure on the input/output vectors is not allowed

[njase]

Hello,

I am doing some experiments on MatrixFree with vector valued FE and have an observation+debug finding. Would appreciate a response to know the underlying reason.

I found this while working on tests/matrix_vector_stokes.cc which works on vectors with block-wise structure (vector of vectors) and test using scalar valued FE. Other test cases like inverse_mass_02.cc use vector-valued FE but use one single vector mf_test.cc.zip

Details:

If used with vector-valued FE, MatrixFree does not allow block-wise structure on the input and output vectors given to cell_loop.

A small test (attached) shows this behaviour, with result as:

  An error occurred in line <2839> of file <xxx/dealii_build/include  /deal.II/matrix_free/fe_evaluation.h> in function
    void dealii::internal::check_vector_compatibility(const VectorType &, const   internal::MatrixFreeFunctions::DoFInfo &) [VectorType = dealii::Vector<double>]
  The violated condition was: 
    (vec.size()) == (dof_info.vector_partitioner->size())
  Additional information: 
    Dimension 125 not equal to 250.

During debugging I find that the function read_dof_values_plain (const VectorType *src[]) in fe_evaluation.h has a comment:

// case with vector-valued finite elements where all components are
// included in one single vector. Assumption: first come all entries to
// the first component, then all entries to the second one, and so
// on. This is ensured by the way MatrixFree reads out the indices.

What is the reason for such an assumption? - Is this feature unsupported, practically undesired or a bug (or)?

Thanks

[kronbichler]

What is the reason for such an assumption? - Is this feature unsupported,

The short answer is yes, this feature is currently unsupported. The reason is simplicity of implementation, because for supporting to split the degrees of freedom in the function read_write_operation requires some additional offsets that we do not currently compute.

The good news is: We already have code that does fix this, but it is not yet committed because it sits inside a big project, the face evaluation part of MatrixFree that lets one evaluate DG operators of all kind. The project is bigger than I would like and there are some parts that do not work yet, which needs to be removed/properly handled, which in turn is why my co-worker Katharina Kormann and I have not yet committed the code. But we plan to do that in the next month. I will set up another issue to track the steps so you can have a look how close it is.

[kronbichler]

For your information, the issue is #5667. I will keep you updated once we have implemented the feature.

[kronbichler]

A workaround now (you maybe already have figured it out from the tests): Create two DoFHandler objects and pass them as a std::vector<const DoFHandler<dim>*> to the MatrixFree::reinit() call, then you can access both components inside a loop.