After some discussions with @nfehn, we decided to start with a general implementation for the Laplace operator involving coefficients. The operator should be capable of handling:

• Scalar and vector valued solution fields (n_components),
• Non-coupling coefficients, i.e. one or n_components many scalars,
• Coupling coefficients, i.e.
  - dyadic coefficients for scalar-valued problems, and
  - tetradic coefficients for vector-valued problems, (although not too likely to be needed soon)
• Coefficients that are variable in space.
• Coefficients that are variable/constant in time.

The implementation is expected to be as generic as possible in terms of rank, but for example differentiate and optimize variable/constant coefficients.

The special case of constant coefficients in space (like unit coefficients) is planned to maintain its own implementation, as some performance overhead is expected due to the flexibility of this generalized operator. Hopefully, all other use cases present or planned in ExaDG can be covered by this operator.

Sounds like a good plan, I am in favor. As discussed offline with @nfehn (and with @peterrum in other occasions), for the base operator (constant coefficient) we need to figure out how precisely realize problem-specific optimizations for Poisson operators that are performance critical and that include our recent state of the art in terms of data locality optimizations, overhead-reductions, etc. and other things that were developed via the CEED BP benchmarks and other contexts. It is yet to be decided whether the right place for these optimizations is here in ExaDG or deal.II, but in either case this fits well with the given strategy of the operator.

As we've agreed, I've been spending some time to create a generalized Laplace operator. I'd like to get your opinion on the matter of BoundaryDescriptor before going on with the implementation.

So far, operators have been using the BoundaryDescriptor in the user interface of their own module. The two "operators" which are not bound to a module are the mass operator and the rhs operator, but they don't have boundary terms.

So my question is where would be the best place to define the new BoundaryDescriptor for the new operator. Should I put it in the header file of the operator, or make a new folder, create a new file and also put everything related to the operator in a new folder, or do something else?

Also, does it make sense at all to define boundary descriptor for a single operator, instead of a problem? I think this is also a design question, but currently the operators need to know about the boundary descriptor.

If you have any opinions / suggestions about this, please let me know. If you want to check out the current state of the implemenation to get a better picture, you can find the branch here.

I also thought about the problem with the BoundaryDescriptor during the last days. I did not come up with a clean idea/solution so far, though. Maybe we could discuss this in person.

The derivation of the DG discretization does not take into account in my opinion that the coefficient can be discontinuous between elements, which needs some averaging/flux computation also for the diffusivity (harmonic mean vs. arithmetic mean). This can be extracted from the ExaDG code for the viscous operator in the IncNS module. Unfortunately, I did not describe this in detail in my PhD thesis, but you can probably find information and pointers to the literature in the papers or the thesis by B. Krank.