Constructing operators in HOG

The construction of operators in the HOG is currently very rigid and cannot mirror more complex cases.

Let's look at an example:

I'd like to implement bilinear forms with both volume and boundary integrals:

\int_\Omega F dx + \int_\Gamma G ds

where F and G are different integrands and either integral integrates over a different subdomain (domain \Omega and boundary \Gamma).

Just as one example of rigidness in the HOG, the element matrix is tied to the HyTeGElementwiseOperator class such that there is no way to have multiple different element matrices being executed after the other in the same operator.

I am using this issue to document/discuss a structure that is more flexible. I have already started refactoring with ce2c107c, but that was not sufficient.

Just some pseudocode to show what I think could be a better structure (incomplete braindump): essentially, most information should be passed to the actual kernels (I may want to be able to optimize them individually).

class KernelType:
    """
    Represents the actual low level kernel over a macro-volume.
    """
    volume_element,
    integration_domain,  # volume, boundary (facets in the future maybe for DG)
    integrand,           # aka form, aka element_matrix
    blending_map,
    quadrature_rule,
    loop_strategy,       # possible choices tied to the integration domain
    optimizations,

class KernelWrapperType:
    """
    Handles pre- and post-communication, switching between 2D/3D, the macro-loop, and calls the actual kernels.

    It also shall detect when to apply kernels that are flagged with integration_domain "boundary" and skip them.
    Possibly we need to supply the BCs for each passed kernel.

    Generated (pseudo)code e.g. looks like this:
 
    void apply():                                        // <- specified by this class
      if domain_is_2D:
        for macro_face in macro_faces:
            pre_communication( ... )                     // <- specified by this class
            some_2D_volume_kernel( ... )                 // <- passed kernel
            if macro.boundary_at_facet(0) == NeumannBC:  // <- conditional generated automatically by this class
                some_2D_boundary_kernel_facet_0( ... )   // <- passed kernel
            if macro.boundary_at_facet(1) == NeumannBC:
                some_2D_boundary_kernel_facet_1( ... )
            ..
            post_communication( ... )                    // <- specified by this class
      else if domain_is_3D:
        for macro_tet in macro_tets:
            pre_communication( ... )                     // <- specified by this class
            some_3D_volume_kernel( ... )                 // <- passed kernel
            ..
     
    """
    operation,     # e.g., apply, assemble, gemv, ... passed to the kernel(?)
    kernels_2D,    # list of kernels - each one corresponding to one integral in the weak form - all of them are added up
    kernels_3D,    # we may have operators for 2D or 3D only - so we set them explicitly
    

class HyTeGElementwiseOperator:
    """
    Represents a single bilinear form (possibly composed of multiple integrals, like the example above) 
    and exposes one or multiple methods that feature that bilinear form in some sense 
    (apply, assemble, diagonal assemle, GEMV, ...)
    """
    trial_space,     # already inside the integrand information currently, 
                     # and likely needed there(?) but should be the same for the entire operator(?) 
                     # -> we could pass this down to the form from here ...
    test_space,      #
    data_type        # real_t, float, double, could actually move somewhere else, ...
    kernel_wrappers  # list of KernelWrapperType

Feedback is more than welcome because I might miss/confuse stuff. This is certainly not the best solution, but got to start somewhere..