Changed extract_dofs so that it works in parallel

[mathsen]

Current implementation of extract_dofs for hp::DoFHandler doesn't work in parallel, since the call to "get_fe()" throws an exception:

The violated condition was:
(dynamic_cast<const dealii::DoFHandler<DoFHandlerType::dimension, DoFHandlerType::space_dimension> *>( this->dof_handler) != nullptr) || (this->is_locally_owned() || this->is_ghost())
Additional information:
You can only query active_fe_index information on cells that are either locally owned or (after distributing degrees of freedom) are ghost cells.

I propose the small attached changes, which make the function run in parallel.

[kronbichler]

Thank you very much for your contribution - this is greatly appreciated. I agree with the patch. I would very much like if we could add a small test case to cover this new feature. Would you be interested in that? Since you appear to be new to the deal.II repository, we would be happy to help you through the infrastructure (or help with inserting it to the test suite in case you want to provide a simple program). Some general info about our tests is found here: https://dealii.org/developer/developers/testsuite.html

[mathsen]

Hello kronbichler,
I have had a look into the test suite and really like it very much. In my opinion the new test would fit best into the "dofs" folder. I wrote one, following some other tests in the dofs and mpi folder, and updated the pull request here with the proposal.
I have another question which came up, when I was looking at the test suite. On my system many tests in the dofs folder didn't pass. I had a closer look at the cause and realized, that it's just because of some "spaces" missing in my locally compiled output files. I have attached an example, the test for dof_tools_output_12 as source files (zipped) and a little screenshot with a relevant extraction.
The relevant code parts in dofs/dof_tools_12.cc 12 are:

check(const FiniteElement<dim> &fe, const std::string &name)
{
  deallog << "Checking " << name << " in " << dim << "d:" << std::endl;
[...]
}

[...]

#define CHECK(EL, deg, dim) \
  {                         \
    FE_##EL<dim> EL(deg);   \
    check(EL, #EL #deg);    \
  }

[...]

#define CHECK_SYS3(sub1, N1, sub2, N2, sub3, N3, dim) \
  {                                                   \
    FESystem<dim> q(sub1, N1, sub2, N2, sub3, N3);    \
    check(q, #sub1 #N1 #sub2 #N2 #sub3 #N3);          \
  }

[...]

// systems of systems
CHECK_SYS3(
  (FESystem<2>(FE_Q<2>(1), 3)), 3, FE_DGQ<2>(0), 1, FE_Q<2>(1), 3, 2);
CHECK_SYS3(FE_DGQ<2>(3),
           1,
           FESystem<2>(FE_DGQ<2>(0), 3),
           1,
           FESystem<2>(FE_Q<2>(2), 1, FE_DGQ<2>(0), 1),
           2,
           2);

To my understanding the spaces should be there, and I don't get why in the reference file the spaces are missing?
Nevertheless,
https://cdash.kyomu.43-1.org/testSummary.php?project=1&name=dofs%2Fdof_tools_12.debug&date=2019-06-04
shows that the test passes on the test suite. On my local machine (same compiler, GCC 8.3) it fails.
Greetings,

Mathias

dof_tools_output_12_comparison.zip

[masterleinad]

It would be helpful to clarify that the output has to be interpreted with respect to

dof_handler.locally_owned_dofs().index_within_set(indices[i])

Otherwise, returning an IndexSet object would probably be a better choice.

[masterleinad]

Currently, you are only ever checking FE_Q elements, i.e. scalar ones with only one component. Thus, the size of mask is one and you are testing the same thing twice.

[mathsen]

Thanks for the review! You are absolutely right, I added vector valued elements and (most) of the elements of dof_tools_common.
I also use now "MPILogInitAll" to get the output of all processes, but, as you said, the selected_dofs vectors are just the locally owned vectors, which are not mapped to global DoF indices. I don't think that this is a problem in this testing scenario, or am I wrong there? Therefore I wouldn't change this, if there is no objection.

Greetings
Mathias

[masterleinad]

No, it's not a problem that you are considering the locally owned degrees of freedom. At least for me, it is just unintuitive that the vector returned doesn't correspond to actual degrees of freedoms but has to be mapped using the IndexSet of locally owned DoFs to make sense (with a distributed mesh). We often define a separate function returning an IndexSet instead of a std::vector<bool> when porting functions from serial Triangulations to distributed ones. Of course, this behavior was pre-existing for the non-hp version and it's fine to extend that behavior to the hp case. Nevertheless, I would appreciate a comment in the documentation on how to interpret the returned object.

[mathsen]

masterleinad: OK, I understand. I explained this fact in the comment of the beginning of the test, or is there another place, where this should go?

[masterleinad]

Sorry for not being precise. Such a comment should go to

dealii/include/deal.II/dofs/dof_tools.h

Lines 1318 to 1350 in 6e6a2bf

	/**
	* Extract the indices of the degrees of freedom belonging to certain vector
	* components of a vector-valued finite element. The @p component_mask
	* defines which components or blocks of an FESystem are to be extracted
	* from the DoFHandler @p dof. The entries in the output array @p
	* selected_dofs corresponding to degrees of freedom belonging to these
	* components are then flagged @p true, while all others are set to @p
	* false.
	*
	* If the finite element under consideration is not primitive, i.e., some or
	* all of its shape functions are non-zero in more than one vector component
	* (which holds, for example, for FE_Nedelec or FE_RaviartThomas elements),
	* then shape functions cannot be associated with a single vector component.
	* In this case, if <em>one</em> shape vector component of this element is
	* flagged in @p component_mask (see
	* @ref GlossComponentMask),
	* then this is equivalent to selecting <em>all</em> vector components
	* corresponding to this non-primitive base element.
	*
	* @param[in] dof_handler The DoFHandler whose enumerated degrees of freedom
	* are to be filtered by this function.
	* @param[in] component_mask A mask that states which components you want
	* to select. The size of this mask must be compatible with the number of
	* components in the FiniteElement used by the @p dof_handler. See
	* @ref GlossComponentMask "the glossary entry on component masks"
	* for more information.
	* @param[out] selected_dofs A vector that will hold @p true or @p false
	* values for each degree of freedom depending on whether or not it
	* corresponds to a vector component selected by the mask above. The size
	* of this array must equal DoFHandler::n_locally_owned_dofs(), which for
	* sequential computations of course equals DoFHandler::n_dofs(). The
	* previous contents of this array are overwritten.
	*/

[mathsen]

Ah, I didn't understand this, it was just "obvious" to me, that this is, what the function does. I have added a little note in the description and hope this is OK.

[masterleinad]

Yes, that looks OK. How did you verify the correctness of the test output? 😉

[mathsen]

I didn't... correct until someone else proofs the opposite ;)

//Edit: To be more precise, the main code change (extract_dofs()) works in a convergence test for me for a real world incompressible flow problem with a vector valued FE. So I strongly assume that the main change here does what it should do. Nevertheless I didn't check the output of the test file.

[masterleinad]

You could probably also unify the hp and the non-hp version here. Would you mind to give it a shot?

[mathsen]

That was what I first wanted to do, because to be honest, the separation between the hp and normal DoFHandler is awkward.
You can have a look at this (failed) try on this github branch:
extract_dofs_parallel_templated_test
My attempt was to generally template the check_this() function, like I did in my 12a file:

template <int dim, typename DoFHandlerType>
void
check_this(const DoFHandlerType &dof_handler)

I edited all the dof_tools... files and changed the implementation of the method to this function header.

And basically add a call in "dof_tools_common" as in 12a:

check_this<dim, DoFHandler<dim>>(dof_handler);
check_this<dim, hp::DoFHandler<dim>>(hp_dof_handler);

So if every check_this() function could just be called for both DoFHandlers, this would be a great solution. Unfortunately not all function calls are yet adjusted to the hp:: case, so my idea was to create a "fake_hp.h" file with a partial template specialization of check_this():

template <int dim>
void
check_this(const hp::DoFHandler<dim> &dof_handler)
{
	// nothing to do here
}

which I would include in all the tests, which are not yet compatible to the hp:: case, so that this "special" do nothing function is being called which just does nothing.
Unfortunately the C++ standard doesn't implement partial function specialization (yet). The situation is well explained here:
function templates partial specialization
and unfortunately function overloading here is not really a nice situation, since you then have again to implement every check_this() method twice (for the hp and non-hp case), which is what I have now done.

If you have any other idea on how to achieve the union I would be very interested in it.

Greetings,
Mathias

//EDIT: Nevertheless my test 12a unifies the hp and non-hp case, actually as I wrote above in my file both DoFHandlers are checked, which more or less doubles test 12. But since I had to remove the 1D elements, it doesn't completely replace test 12, that's why I named it 12a and didn't just replace 12.

[masterleinad]

Ah, I was not talking about the tests, but rather DoFTools::extract_dofs. You should be able to use the hp-implementation also for the non-hp case.

[mathsen]

Ah, I see, this should indeed be possible. I will give it a try in the afternoon :)

[mathsen]

And I implemented a "united version" for both DoFHandlers with templates. As I wrote in the commit auto deduction now is with this solution a little problem and I had to change other locations in the library, too that call extract_dofs().
But this was the first idea how to merge these two functions into one.
Greetings

Mathias

[masterleinad]

You can use DoFHandlerType::dimension and DoFHandlerType::space_dimension here and in the other places.

[mathsen]

Thanks for this valuable hint! Indeed this makes life so much easier.
I have (again...) a new proposal now and could merge all the ideas I had at the beginning. Let me explain in a nutshell what I tried to achieve and I begin with the testsuite (and I like it really very much):
Basically I splitted the 12a test into a dof_tools_common_parallel.h file, which is similar to the non parallel version, just with the parallel triangulation and not using (the non supported) 1d elements. So basically if someone wants to make a test for MPI in the future, he can just include this file, similar to dof_tools_common.
Further I made the check_this() template dependent of the DoFHandler type - also this means to change many code lines. The benefit is, that dof_cools_common.h can now test for both DoFHandlers (hp and not hp) without code deduplicating. I demonstrate this with the updated test dof_tools_12: it's basically the same code as before, but now tests against the normal DoFHandler and the hp::DoFHandler. The other tests could just do the same, if every function implemented the full functionality for DoFHandler and hp::DoFHandler. Until this they just have to include the dof_tools_fake_hp.inst.in so that basically just the normal DoFHandler is being tested. So in my opinion this will make testing in the future easier, when more hp::DoFHandler stuff gets merged (maybe some tests already work with the hp::DoFHandler, I didn't really look into this).

Now to the extract_dofs() function: as you proposed I changed it, that it doesn't depend on dim or spacedim template arguments and therefore auto deduction works as expected.

I still feature the problem with the different output files on my local machine and some reference output files, what I describe above - I really don't know where the different outputs come from and I don't know, why the whitespaces aren't there in the reference output files. This could be a problem with the output file of test 12a, since it includes the whitespaces in the FESystems, also the updated test 12 output files from my system include these whitespaces.
Greetings,
Mathias

[masterleinad]

/rebuild

[masterleinad]

This looks mostly good to me. Thanks for following my suggestions!
I have just a few more things:

• Please rename tests/dofs/dof_tools_12a.output to tests/dofs/dof_tools_12a.mpirun=1.output. Otherwise, the test doesn't pass if deal.II is built without MPI support.
• Please write a change log entry in doc/news/changes/minor/ describing your changes.
• We use the file extension inst.in only for generating inst files for explicit instantiations. Can you rename tests/dofs/dof_tools_fake_hp.inst.in to tests/dofs/dof_tools_common_fake_hp.h?
• The newly created files should just have 2019 as copyright year.

After that, we should be good to go.

[mathsen]

Next try :)
I changed your proposals.
I'm still not getting the cause for the different *.output files and think that could be a showbreaker when testing the changes.
Greetings,
Mathias

[masterleinad]

I'm still not getting the cause for the different *.output files and think that could be a showbreaker when testing the changes.

The tester is fine with it. 🙂

[mathsen]

The tester is fine with it. slightly_smiling_face

I skipped through the output of the tester and at least I didn't see that the test for dof_tools_12a was executed (neither in the mpi nor in the serial tester). I also didn't find the output for the other parallel test files in the dofs folder, e.g. the dof_tools_04 test. Do I just miss this or do you have to do something else that the tests are being executed?

Greetings

Mathias

[masterleinad]

I skipped through the output of the tester and at least I didn't see that the test for dof_tools_12a was executed (neither in the mpi nor in the serial tester). I also didn't find the output for the other parallel test files in the dofs folder, e.g. the dof_tools_04 test. Do I just miss this or do you have to do something else that the tests are being executed?

I ran the test locally and it passed. We don't run all the tests on the CI tester.

[masterleinad]

Looks good to me! Let's see if someone else has something to say.

[masterleinad]

I am having problems with this on older compilers. Can you please replace it here and in all the other places?

Suggested change

	std::map<types::boundary_id, const Function<dof_handler.dimension> *>
	std::map<types::boundary_id, const Function<DoFHandlerType::dimension> *>

[mathsen]

Fixed, I hope I greped all occurrences :)