adds bilinear-form for a penalty term with arguments L2 projected on the constants

src/hyteg/hdgp1functionspace/HDGP1LocalAssembly.hpp

@@ -259,51 +259,105 @@ inline void jump( const std::array< Point3D, 6 >& points,
    const uint_t difference  = ( levelUp > levelDown ) ? ( levelUp - levelDown ) : ( levelDown - levelUp );
    const auto   numSubEdges = real_t( 1u << difference );
 
-   real_t leftPtUp, rightPtUp, leftPtDown, rightPtDown, edgeLength;
+   real_t leftPtUp, rightPtUp, leftPtDown, rightPtDown, alphaTimesEdgeLength;
    if ( levelUp > levelDown )
    {
-      leftPtUp    = 0.;
-      rightPtUp   = 1.;
-      leftPtDown  = real_t( offset ) / numSubEdges;
-      rightPtDown = real_t( offset + 1 ) / numSubEdges;
-      edgeLength  = alpha * ( points[0] - points[1] ).norm();
+      // Situation:
+      //     Up-triangle
+      //         0   1
+      //         +---+
+      // +---o---o---o---+
+      // 0  1/4 2/4 3/4  1
+      //     Down-triangle
+      leftPtUp             = 0.;
+      rightPtUp            = 1.;
+      leftPtDown           = real_t( offset ) / numSubEdges;
+      rightPtDown          = real_t( offset + 1 ) / numSubEdges;
+      alphaTimesEdgeLength = alpha * ( points[0] - points[1] ).norm();
    }
    else
    {
-      leftPtUp    = real_t( offset ) / numSubEdges;
-      rightPtUp   = real_t( offset + 1 ) / numSubEdges;
-      leftPtDown  = 0.;
-      rightPtDown = 1.;
-      edgeLength  = alpha * ( points[3] - points[4] ).norm();
+      // Reverse situation:
+      //     Up-triangle
+      // 0  1/4 2/4 3/4  1
+      // +---o---o---o---+
+      //         +---+
+      //         0   1
+      //     Down-triangle
+      leftPtUp             = real_t( offset ) / numSubEdges;
+      rightPtUp            = real_t( offset + 1 ) / numSubEdges;
+      leftPtDown           = 0.;
+      rightPtDown          = 1.;
+      alphaTimesEdgeLength = alpha * ( points[3] - points[4] ).norm();
    }
 
    const real_t middlePtUp   = 0.5 * ( leftPtUp + rightPtUp );
    const real_t middlePtDown = 0.5 * ( leftPtDown + rightPtDown );
 
-   // Simpson's rule
-   tensor[index2x4( 0, 0 )] += +edgeLength / 6 *
-                               ( ( 1 - leftPtUp ) * ( 1 - leftPtUp ) + 4. * ( 1 - middlePtUp ) * ( 1 - middlePtUp ) +
-                                 ( 1 - rightPtUp ) * ( 1 - rightPtUp ) );
-   tensor[index2x4( 0, 1 )] +=
-       +edgeLength / 6 *
-       ( ( 1 - leftPtUp ) * ( leftPtUp ) + 4. * ( 1 - middlePtUp ) * ( middlePtUp ) + ( 1 - rightPtUp ) * ( rightPtUp ) );
-   tensor[index2x4( 0, 2 )] += -edgeLength / 6 *
-                               ( ( 1 - leftPtUp ) * ( 1 - leftPtDown ) + 4. * ( 1 - middlePtUp ) * ( 1 - middlePtDown ) +
-                                 ( 1 - rightPtUp ) * ( 1 - rightPtDown ) );
-   tensor[index2x4( 0, 3 )] +=
-       -edgeLength / 6 *
-       ( ( 1 - leftPtUp ) * ( leftPtDown ) + 4. * ( 1 - middlePtUp ) * ( middlePtDown ) + ( 1 - rightPtUp ) * ( rightPtDown ) );
-   tensor[index2x4( 1, 0 )] +=
-       +edgeLength / 6 *
-       ( ( leftPtUp ) * ( 1 - leftPtUp ) + 4. * ( middlePtUp ) * ( 1 - middlePtUp ) + ( rightPtUp ) * ( 1 - rightPtUp ) );
-   tensor[index2x4( 1, 1 )] +=
-       +edgeLength / 6 * ( ( leftPtUp ) * ( leftPtUp ) + 4. * ( middlePtUp ) * ( middlePtUp ) + ( rightPtUp ) * ( rightPtUp ) );
-   tensor[index2x4( 1, 2 )] +=
-       -edgeLength / 6 *
-       ( ( leftPtUp ) * ( 1 - leftPtDown ) + 4. * ( middlePtUp ) * ( 1 - middlePtDown ) + ( rightPtUp ) * ( 1 - rightPtDown ) );
-   tensor[index2x4( 1, 3 )] +=
-       -edgeLength / 6 *
-       ( ( leftPtUp ) * ( leftPtDown ) + 4. * ( middlePtUp ) * ( middlePtDown ) + ( rightPtUp ) * ( rightPtDown ) );
+   // Applying Simpson's rule to integrate the hat functions:
+   // (the comments below are just here to retain useful formatting.
+   tensor[index2x4( 0, 0 )] += +alphaTimesEdgeLength / 6 *
+                               ( +( 1 - leftPtUp ) * ( 1 - leftPtUp )           //
+                                 + 4. * ( 1 - middlePtUp ) * ( 1 - middlePtUp ) //
+                                 + ( 1 - rightPtUp ) * ( 1 - rightPtUp ) );
+   tensor[index2x4( 0, 1 )] += +alphaTimesEdgeLength / 6 *
+                               ( +( 1 - leftPtUp ) * ( leftPtUp )           //
+                                 + 4. * ( 1 - middlePtUp ) * ( middlePtUp ) //
+                                 + ( 1 - rightPtUp ) * ( rightPtUp ) );
+   tensor[index2x4( 0, 2 )] += -alphaTimesEdgeLength / 6 *
+                               ( +( 1 - leftPtUp ) * ( 1 - leftPtDown )           //
+                                 + 4. * ( 1 - middlePtUp ) * ( 1 - middlePtDown ) //
+                                 + ( 1 - rightPtUp ) * ( 1 - rightPtDown ) );
+   tensor[index2x4( 0, 3 )] += -alphaTimesEdgeLength / 6 *
+                               ( +( 1 - leftPtUp ) * ( leftPtDown )           //
+                                 + 4. * ( 1 - middlePtUp ) * ( middlePtDown ) //
+                                 + ( 1 - rightPtUp ) * ( rightPtDown ) );
+   tensor[index2x4( 1, 0 )] += +alphaTimesEdgeLength / 6 *
+                               ( +( leftPtUp ) * ( 1 - leftPtUp )           //
+                                 + 4. * ( middlePtUp ) * ( 1 - middlePtUp ) //
+                                 + ( rightPtUp ) * ( 1 - rightPtUp ) );
+   tensor[index2x4( 1, 1 )] += +alphaTimesEdgeLength / 6 *
+                               ( ( leftPtUp ) * ( leftPtUp )            //
+                                 + 4. * ( middlePtUp ) * ( middlePtUp ) //
+                                 + ( rightPtUp ) * ( rightPtUp ) );
+   tensor[index2x4( 1, 2 )] += -alphaTimesEdgeLength / 6 *
+                               ( +( leftPtUp ) * ( 1 - leftPtDown )           //
+                                 + 4. * ( middlePtUp ) * ( 1 - middlePtDown ) //
+                                 + ( rightPtUp ) * ( 1 - rightPtDown ) );
+   tensor[index2x4( 1, 3 )] += -alphaTimesEdgeLength / 6 *
+                               ( +( leftPtUp ) * ( leftPtDown )           //
+                                 + 4. * ( middlePtUp ) * ( middlePtDown ) //
+                                 + ( rightPtUp ) * ( rightPtDown ) );
+}
+
+/// This function assembles the terms
+///    alpha * inner(jump(P(u)), jump(P(v)))*dS
+/// over the interior boundary, were alpha is some coefficient and 
+/// P the L2-projection on the space of constant functions on the edges.
+/// It assumes a vertex ordering of the form
+///     0             1
+///     +-------------+
+///     3             4
+/// where the indices above (0 and 1) are the local indices of the upper cell,
+/// and the indices below (3 and 4) are the local indices of the lower cell.
+inline void jumpProjected( const std::array< Point3D, 6 >& points,
+                           const uint_t&                   levelUp,
+                           const uint_t&                   levelDown,
+                           const uint_t,
+                           const real_t& alpha,
+                           EdgeTensor&   tensor )
+{
+   // we integrate over the smaller edge which belongs to a face with a higher refinement level.
+   auto const edgeLength           = levelUp > levelDown ? ( points[0] - points[1] ).norm() : ( points[3] - points[4] ).norm();
+   auto const alphaTimesEdgeLength = alpha * edgeLength;
+   tensor[index2x4( 0, 0 )] += +alphaTimesEdgeLength / 4 * ( +1 );
+   tensor[index2x4( 0, 1 )] += +alphaTimesEdgeLength / 4 * ( +1 );
+   tensor[index2x4( 0, 2 )] += +alphaTimesEdgeLength / 4 * ( -1 );
+   tensor[index2x4( 0, 3 )] += +alphaTimesEdgeLength / 4 * ( -1 );
+   tensor[index2x4( 1, 0 )] += +alphaTimesEdgeLength / 4 * ( +1 );
+   tensor[index2x4( 1, 1 )] += +alphaTimesEdgeLength / 4 * ( +1 );
+   tensor[index2x4( 1, 2 )] += +alphaTimesEdgeLength / 4 * ( -1 );
+   tensor[index2x4( 1, 3 )] += +alphaTimesEdgeLength / 4 * ( -1 );
 }
 
 /// Creates the matrix for the affine transformation to our counter-clockwise reference triangle.
@@ -735,6 +789,7 @@ class VariableMacroEdgeForm
       const real_t delta = Transpose ? eps_ : delta_;
       const real_t eps   = Transpose ? delta_ : eps_;
       jump( points, levelUp, levelDown, offset, sigma_ * hInv, tensor );
+      // jumpProjected( points, levelUp, levelDown, offset, sigma_ * hInv, tensor );
       jumpAvg< UpperFaceOrientatedCW >( points, levelUp, levelDown, offset, eps, delta, kappa_, tensor );
    }