Apply boundary integral to diagonal in MatrixFree

[yy.wayne]

Hi there,

I've been applying MatrixFree method to Step-29 where face integrals are involved, different from step-37, step-75 and so on. Firstly I think loop() instead of cell_loop() should be used. Secondly, similar to step-75, Chebyshev smoother and a coarsest level matrix are needed in my code. In step-75 we get diagonal vector by MatrixFreeTools::compute_diagonal and MatrixFreeTools::compute_matrix.

However, the face integral makes them different. Both two functions only take cell_operations(FEEvaluation), not face_operations(FEFaceEvaluation). My question is how to compute diagonal and matrix when my weak form has face integrals terms (only boundary integral, no interface integral). 

[yy.wayne]

To write the boundary_integral in MatrixFree for step-29 is kind of confusing as well...

For FEFaceEvaluation with 2 components(real and imaginary parts), whether a index correspond to src(alpha and beta) or basis function(phi_0 for phi and phi_1 for psi) is a problem. I cannot test this code yet due to the former question.

[yy.wayne]

Maybe I should describe my question again.

In short, to use Chebyshev smoother and solve mg coarse grid straightly, we need MatrixFreeTools::compute_diagonal and MatrixFreeTools::compute_matrix, respectively.

However, for problems involve Robin boundary condition where integral is required, those 2 function will not work(I think). That's because FEFaceEvaluation should be included as well.

I think I will write new functions if there's no appropriate one. Any hint, especially on FEFaceEvaluation and multi-component MatrixFree formulation is welcomed.

[Peter Munch]

Hi,

you are right such a function is unfortunately missing in deal.II. It has been on my TODO list for a long time but I didn't get around to implement it, since I haven't had to use DG for a while and in most of my applications I can move the BCs to the RHS.

Setting up the matrix is not particular difficult. You can take a look at our CFD code (https://github.com/exadg/exadg/blob/66a38ee6c0a03bbcc003ff4c19d94354face492b/include/exadg/operators/operator_base.cpp#L801-L827, https://github.com/exadg/exadg/blob/66a38ee6c0a03bbcc003ff4c19d94354face492b/include/exadg/operators/operator_base.cpp#L1681-L1967). Computing the diagonal is a bit more difficult. What the function in deal.II does is to compute a column of the element stiffness matrix and to apply right away the constraints without ever having to assemble the system matrix. I am not sure how much work is to extend the function so that it considers also face integrals; I think the logic should be similar.

>  think I will write new functions if there's no appropriate one.

That would be great!

Hope this help,
Peter