Enable hanging nodes for continuous Galerkin discretizations

[nfehn]

This PR aims to realize local refinements / hanging-nodes for continuous Galerkin discretizations for nonlinear problems (e.g. nonlinear structural mechanisms) with special requirements in terms of constraints / inhomogeneous BC as discussed in issue #386.

closes #386

closes #446 (@bugrahantemur please take a look at the changes and let me know if you do not agree)

TODOs:

• introduce second AffineConstraints object for Poisson spatial discretization (which is the second case of continuous Galerkin discretizations in ExaDG apart from the Structure module; other modules use continuous Galerkin only inside MG, which is unproblematic since the operators used in MG preconditioning are always homogeneous)
• set_inhomogeneous_boundary_values() will be removed from OperatorBase::evaluate() functions. Make sure that all modules currently calling these functions are adapted accordingly.
• initialize dof_index_inhomogeneous with an invalid number and make sure that a correct value has been set whenever using dof_index_inhomogeneous

[nfehn]

@kronbichler the lengthy comment in commit f02f2c1 is kind of a summary of our discussion from #386 (... or what I took home from this discussion). Could you have a look at this commit?

[nfehn]

I think I have to resolve #507 first, because we need the separate dof-index and integrator also in OperatorBase for functions of type rhs_add() for example.

[nfehn]

The code seems to crash in OperatorBase::add_diagonal() with the error message

An error occurred in line <7847> of file </home/code/dealii/include/deal.II/matrix_free/fe_evaluation.h> in function
    void dealii::FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number, VectorizedArrayType>::evaluate(dealii::EvaluationFlags::EvaluationFlags) [with int dim = 3; int fe_degree = -1; int n_q_points_1d = 0; int n_components_ = 1; Number = float; VectorizedArrayType = dealii::VectorizedArray<float, 8>]
The violated condition was: 
    this->dof_values_initialized == true
Additional information: 
    (none)

What about the TODO here:

exadg/include/exadg/operators/operator_base.cpp

Lines 432 to 448 in 980cc10

	dealii::MatrixFreeTools::
	compute_diagonal<dim, -1, 0, n_components, Number, dealii::VectorizedArray<Number>>(
	*matrix_free,
	diagonal,
	[&](auto & integrator) -> void {
	// TODO this line is currently needed as bugfix, but should be
	// removed because reinit is now done twice
	this->reinit_cell(integrator, integrator.get_current_cell_index());
	integrator.evaluate(integrator_flags.cell_evaluate);
	this->do_cell_integral(integrator);
	integrator.integrate(integrator_flags.cell_integrate);
	},
	data.dof_index,
	data.quad_index);

The line with this->reinit_cell() in the above code seems to be problematic. Without that line, the code does not crash (for Poisson). However, for general operators, this function is really necessary, as it may be overwritten by derived classes, see e.g. ElasticityOperatorBase::reinit_cell(), where the material handler is set in this function. Hence, we cannot simply remove this line. One can see this as design problem on the deal.II side, or on the ExaDG side. I think regarding dealii, the function dealii::MatrixFreeTools::compute_diagonal should receive a separate, user-specified lambda for the reinit part, in order to have more flexibility. I might be able to fix the problem on the ExaDG side, but this makes the code more ugly and not easier to understand I fear. Another option would be to simply implement this code for diagonal computation in ExaDG for the continuous Galerkin case, just as we do it in ExaDG for discontinuous Galerkin. @kronbichler, what do you think?

[kronbichler]

The code inside deal.II is complicated, because it tracks down the constraints and computes their influence on the local degrees of freedom, doing a long preparation https://github.com/dealii/dealii/blob/3bba74760bff6724e1766bf22e4ee06eaf98bfc1/include/deal.II/matrix_free/tools.h#L422-L854 and then https://github.com/dealii/dealii/blob/3bba74760bff6724e1766bf22e4ee06eaf98bfc1/include/deal.II/matrix_free/tools.h#L872-L908 to extract the entries. I recently rewrote around 300 of the 500 lines of that code in dealii/dealii#15135 and dealii/dealii#15147 because the initial implementation was performing poorly - the details are not so important, but it shows that re-writing this in ExaDG is probably not a good idea.

I see two options here:

• We equip deal.II with a lambda that triggers once per cell next to the phi->reinit(cell) inside the first link I posted. This leads to a second lambda function next to the operation to be provided or where we can provide a default. The problem here is that if we want to support both deal.II-9.5 and master, we need an ifdef and still option 2.
• We keep track of what cell has been used before,

  unsigned int last_cell = numbers::invalid_unsigned_int;

 dealii::MatrixFreeTools:: 
   compute_diagonal<dim, -1, 0, n_components, Number, dealii::VectorizedArray<Number>>( 
     *matrix_free, 
     diagonal, 
     [&](auto & integrator) -> void { 
       if (last_cell != integrator.get_current_cell_index())
         this->reinit_cell(integrator, integrator.get_current_cell_index()); 
       ...

and this way only call it once per cell, not once per column. We still use another evaluator here, but that one is hidden inside compute_diagonal and nothing we can do something about.

I usually do not like lambda functions, because they are also complicated for a user to understand. But here it might make sense to go down that way because it is a bit cleaner for ExaDG nonetheless in my opinion.

[nfehn]

We should only do it once per cell and not once per column. So the lambda function I was thinking of would be called by the deal.II code only once per cell, just as integrator.reinit(cell).

Overall this is in line with my general design ideas/wishes regarding generic marix-free operators: The user wants to specify (i) reinit() functions and (ii) what needs to be done on a quadrature point (the "kernel" operation). The evaluate/integrate steps as well as the quadrature point loop should be generic code not visible in user code ... At least this is what I would formulate as a design goal for a generic matrix-free code.

[kronbichler]

Overall this is in line with my general design ideas/wishes regarding generic marix-free operators: The user wants to specify (i) reinit() functions and (ii) what needs to be done on a quadrature point (the "kernel" operation).

Yes indeed, this is what one would like to have, fully agree.

[nfehn]

I would like to avoid the above variant with last_cell, because the code is hard to understand from the user perspective if one does not know the implementation in deal.II. The problem is that integrator comes from the deal.II code, while last_cell is injected/catched from the ExaDG code and the body of the lambda does not only depend on both of them, it may also change both variables (i.e. they are not const). I think this is exactly what makes lambdas hard-to-understand to error-prone, so I agree with the general comment on lambdas. One should not mis- or over-use them.

[nfehn]

@kronbichler I wonder why the tests have passed in PR #508 given that the line with reinit_cell() has been changed in PR #508.

[nfehn]

@kronbichler I think I resolved the reinit_cell() problem for now in ExaDG. I also updated the multigrid preconditioner for Structure (adding another AffineConstraints object). Ready for review from my side ...

[kronbichler]

I mostly agree with this PR. Below, I list some cosmetic changes, plus some code design issues regarding the setup of constraints, plus one discussion as a point of documentation.

[nfehn]

Note that the utility functions introduced in c3dc80b could also be used to avoid code duplications in the multigrid implementation. However, I suggest to do this as a follow-up PR (since this might touch many files).

[nfehn]

this line seems to have introduced a bug in debug mode, see PR #535, since we later call copy_from and then add additional constraints.

[nfehn]

The solution is to call close() in the end once all the constraints have been added (which somewhat hinders writing modular code ...)