introduces more general boundary conditions

apps/CMakeLists.txt

@@ -87,3 +87,5 @@ add_subdirectory(SphericalHarmonics)
 add_subdirectory(Surrogates)
 
 add_subdirectory(blockOperators)
+
+add_subdirectory(RobinLaplace)

apps/RobinLaplace/CMakeLists.txt

+if(NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/output)
+    file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/output")
+endif()
+
+waLBerla_link_files_to_builddir( *.prm )
+
+waLBerla_add_executable( NAME RobinLaplace
+        FILES  RobinLaplace.cpp
+        DEPENDS hyteg)
\ No newline at end of file

apps/RobinLaplace/RobinLaplace.cpp

+#include "core/mpi/Environment.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/mpi/MPIManager.h"
+ 
+#include "hyteg/mesh/MeshInfo.hpp"
+#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
+#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
+#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
+#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
+
+#include "hyteg/solvers/CGSolver.hpp"
+#include "hyteg/solvers/MinresSolver.hpp"
+
+#include "hyteg/dataexport/VTKOutput.hpp"
+#include "hyteg/primitivestorage/Visualization.hpp"
+#include "hyteg/gridtransferoperators/P1toP1LinearProlongation.hpp"
+#include "hyteg/gridtransferoperators/P1toP1LinearRestriction.hpp"
+#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
+#include "hyteg/p1functionspace/P1Function.hpp"
+#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
+#include "hyteg/solvers/GaussSeidelSmoother.hpp"
+#include "hyteg/solvers/GeometricMultigridSolver.hpp"
+
+#include "helpers/P1OperatorWithBC.hpp"
+#include "helpers/P1OperatorWithBC.cpp"
+#include "helpers/P1CustomBoundaryOperator.hpp"
+#include "helpers/P1CustomBoundaryOperator.cpp"
+ 
+using walberla::real_c;
+using walberla::real_t;
+using namespace hyteg;
+
+void generateRightRobinSide(  const std::shared_ptr< PrimitiveStorage >&        storage,  
+                              hyteg::P1Function< real_t >&                      rhsFunction,
+                              std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
+                              uint_t                                            level,
+                              const DoFType&                                    flag )
+{
+  std::vector< real_t > quadratureNodes_ = { 0.5 };
+  std::vector< real_t > quadratureWeights_ = { 1.0 };
+  communication::syncFunctionBetweenPrimitives( rhsFunction, level );
+  rhsFunction.communicate< Vertex, Edge >( level );
+  rhsFunction.communicate< Edge, Vertex >( level );
+
+  std::vector< hyteg::PrimitiveID > edgeIDs = storage->getEdgeIDs();
+  for ( int i = 0; i < int_c( edgeIDs.size() ); i++ )
+  {
+    hyteg::Edge& edge = *storage->getEdge( edgeIDs[uint_c( i )] );
+
+    const DoFType edgeBC = rhsFunction.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
+    if ( testFlag( edgeBC, flag ) )
+    {   
+      for( const auto& microEdgeIndices: hyteg::indexing::MicroEdgeIterator( levelinfo::num_microedges_per_edge( level ) ))
+      {  
+        real_t edgeNum = real_t( levelinfo::num_microedges_per_edge( level ) ); 
+        Point3D microEdgeLength = (edge.getCoordinates()[1] - edge.getCoordinates()[0]) / edgeNum;
+        Point3D coord0 = edge.getCoordinates()[0] + real_t( microEdgeIndices[0].col() ) * microEdgeLength;
+        Point3D coord1 = edge.getCoordinates()[0] + real_t( microEdgeIndices[1].col() ) * microEdgeLength;
+
+        real_t integralMicro0 = 0;
+        real_t integralMicro1 = 0;
+
+        for(uint_t j = 0; j < quadratureNodes_.size(); j++) 
+        {
+          Point3D currentNodeCoord({ coord0[0] * quadratureNodes_[j] + coord1[0] * (1 - quadratureNodes_[j]),
+            coord0[1] * quadratureNodes_[j] + coord1[1] * (1 - quadratureNodes_[j])});
+          real_t boundaryValue = boundaryFunction( currentNodeCoord );
+
+          integralMicro0 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[j];
+          integralMicro1 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[j]);   
+        }
+
+        PrimitiveDataID< FunctionMemory< real_t >, Edge > rhsEdgeDoFIdx = rhsFunction.getEdgeDataID();
+        auto rhsData = edge.getData( rhsEdgeDoFIdx )->getPointer( level );
+        rhsData[ microEdgeIndices[0].col() ] += integralMicro0;
+        rhsData[ microEdgeIndices[1].col() ] += integralMicro1;
+      }
+      rhsFunction.communicateAdditively< Edge, Vertex >( level, hyteg::NeumannBoundary, *storage, false );
+      rhsFunction.communicateAdditively< Edge, Face >( level, hyteg::NeumannBoundary, *storage, false );
+    }
+  }
+}
+
+int main( int argc, char** argv ) 
+{
+  walberla::Environment env( argc, argv );
+  walberla::mpi::MPIManager::instance()->useWorldComm();
+  uint_t numProcesses = uint_c( walberla::mpi::MPIManager::instance()->numProcesses() );
+
+  walberla::shared_ptr< walberla::config::Config > cfg( new walberla::config::Config );
+  cfg->readParameterFile( "./RobinLaplace.prm" );
+  walberla::Config::BlockHandle parameters = cfg->getOneBlock( "Parameters" );
+  parameters.listParameters();
+
+  const uint_t    minLevel        = parameters.getParameter< uint_t >( "minLevel" );
+  const uint_t    maxLevel        = parameters.getParameter< uint_t >( "maxLevel" );
+  const uint_t    maxKrylowDim    = parameters.getParameter< uint_t >( "maxKrylowDim" );
+  const real_t    robinFactor     = parameters.getParameter< real_t >( "robinFactor" );
+  const bool      useVTK          = parameters.getParameter< bool >( "vtkOutput" ); 
+    
+  hyteg::MeshInfo meshInfo = hyteg::MeshInfo::meshRectangle( hyteg::Point2D( {0.0, 0.0} ), hyteg::Point2D( {1.0, 1.0} ), hyteg::MeshInfo::CRISS, 2, 2 );
+  hyteg::SetupPrimitiveStorage setupStorage( meshInfo, numProcesses );
+  setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, false );
+  hyteg::loadbalancing::roundRobin( setupStorage );
+  std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage ); 
+
+  if( useVTK )
+  {
+    hyteg::writeDomainPartitioningVTK( storage, "./output", "RobinLaplace_domain" );
+  }
+
+  std::function< real_t( const hyteg::Point3D& ) > exactSolutionInterpolate = [] ( const hyteg::Point3D& x ) 
+  {
+    return std::pow( x[0], 2 ) - x[0]*x[1] - std::pow( x[1], 2 ) + x[1] + 1.0;
+  };
+
+  std::function< real_t( const hyteg::Point3D& ) > boundaryFunction = [ robinFactor ]( const hyteg::Point3D& x ) 
+  {   
+    real_t boundaryTOL = 0.05;
+    real_t boundaryEvaluate = 0.0;
+    if( x[0] >= 1-boundaryTOL ) 
+    {
+      boundaryEvaluate =  2.0 - x[1] + robinFactor * (2.0 - x[1]*x[1]);
+    } 
+    else if( x[1] >= 1-boundaryTOL ) 
+    {
+      boundaryEvaluate = -x[0] - 1.0 + robinFactor * (x[0]*x[0] - x[0] + 1.0);
+    } 
+    else if( x[0] <= boundaryTOL )
+    {
+      boundaryEvaluate = x[1] + robinFactor * (-x[1]*x[1] + x[1] + 1.0);
+    } 
+    else if( x[1] <= boundaryTOL )
+    {
+      boundaryEvaluate = x[0] - 1.0 + robinFactor * (x[0]*x[0] + 1.0);
+    } 
+    return boundaryEvaluate;
+  };
+
+  std::vector< real_t > boundaryQuadNodes = { 0.5 };
+  std::vector< real_t > boundaryQuadWeights = { 1.0 };
+
+  hyteg::P1ConstantLaplaceOperator innerOperator( storage, minLevel, maxLevel );
+  P1CustomBoundaryOperator boundaryOperator( storage, minLevel, maxLevel, robinFactor, boundaryFunction, boundaryQuadNodes, boundaryQuadWeights );
+
+  P1OperatorWithBC< P1ConstantLaplaceOperator, P1CustomBoundaryOperator > CustomLaplaceOperator( storage, minLevel, maxLevel, innerOperator, boundaryOperator );
+  hyteg::P1Function< real_t > f( "f", storage, minLevel, maxLevel ); 
+  f.interpolate( 0.0, maxLevel, hyteg::NeumannBoundary );
+
+  std::vector< PrimitiveID > vertexIDs = storage->getVertexIDs();
+  for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
+  {
+    Vertex& vertex = *storage->getVertex( vertexIDs[uint_c( i )] );
+    PrimitiveDataID< FunctionMemory< real_t >, Vertex > fVertexDoFIdx = f.getVertexDataID();
+    auto fData = vertex.getData( fVertexDoFIdx )->getPointer( maxLevel );
+    real_t fAtVertex = fData[0];
+
+    std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " : f( vertex.getCoordinates() ) = " << fAtVertex << std::endl;
+  } 
+
+  hyteg::P1Function< real_t > u( "u", storage, minLevel, maxLevel ); 
+  u.interpolate( 0.5, maxLevel, hyteg::DirichletBoundary );
+
+  hyteg::VTKOutput vtkOutput("./output", "RobinLaplace", storage);
+  if( useVTK )
+  {
+    vtkOutput.add( u );
+    vtkOutput.add( f );
+    vtkOutput.write( maxLevel, 0 );
+  } 
+
+  hyteg::MinResSolver minresSolver = hyteg::MinResSolver< P1OperatorWithBC< hyteg::P1ConstantLaplaceOperator, P1CustomBoundaryOperator > >( storage, minLevel, maxLevel, maxKrylowDim );
+  minresSolver.solve( CustomLaplaceOperator, u, f, maxLevel );
+
+  if( useVTK )
+  {
+    vtkOutput.write( maxLevel, 1 );
+  }
+
+  hyteg::P1Function< real_t > exactSolution( "exactSolution", storage, minLevel, maxLevel );
+  exactSolution.interpolate( exactSolutionInterpolate, maxLevel, hyteg::All );
+  std::cout << "sup( exactSolution ) = " << exactSolution.getMaxMagnitude( maxLevel, hyteg::All, false ) << std::endl;
+  std::cout << "sup( u|complete ) = " << u.getMaxMagnitude( maxLevel, hyteg::All, false ) << std::endl;
+  std::cout << "sup( u|boundary ) = " << u.getMaxMagnitude( maxLevel, hyteg::NeumannBoundary, false ) << std::endl;
+  hyteg::P1Function< real_t > residual( "residual", storage, minLevel, maxLevel );
+  residual.assign( {1.0, -1.0}, {u, exactSolution}, maxLevel, hyteg::All );
+  std::cout << " sup norm of residual (complete)        : " << residual.getMaxMagnitude( maxLevel, hyteg::All, false ) << std::endl;
+  std::cout << " sup norm of residual (boundary)        : " << residual.getMaxMagnitude( maxLevel, hyteg::NeumannBoundary, false ) << std::endl;
+  real_t numDoFs = real_t( hyteg::numberOfGlobalDoFs<P1FunctionTag>( *storage, maxLevel));
+  std::cout << " p1-norm per DoF of residual (complete) : " << (residual.dotGlobal( residual, maxLevel, hyteg::All ) / numDoFs) << std::endl;
+
+  std::vector< PrimitiveID > vertexIDsDebug = storage->getVertexIDs();
+  for ( int i = 0; i < int_c( vertexIDsDebug.size() ); i++ )
+  {
+    Vertex& vertex = *storage->getVertex( vertexIDsDebug[uint_c( i )] );
+    PrimitiveDataID< FunctionMemory< real_t >, Vertex > exactVertexDoFIdx = exactSolution.getVertexDataID();
+    PrimitiveDataID< FunctionMemory< real_t >, Vertex > uVertexDoFIdx = u.getVertexDataID();
+    auto exactData = vertex.getData( exactVertexDoFIdx )->getPointer( maxLevel );
+    auto uData = vertex.getData( uVertexDoFIdx )->getPointer( maxLevel );
+    real_t exactAtVertex = exactData[0];
+    real_t uAtVertex = uData[0];
+
+    std::cout << "vertexID = " << vertexIDsDebug[uint_c( i )] << " , vertex.getCoordinates() = " << vertex.getCoordinates() << 
+      " : exactSolution( vertex.getCoordinates() ) = " << exactAtVertex << " , u( vertex.getCoordinates() ) = " << uAtVertex << std::endl; 
+  }  
+
+}
+
+

apps/RobinLaplace/RobinLaplace.prm

+Parameters
+{
+    minLevel 0;
+    maxLevel 0;
+    maxKrylowDim 5;
+    robinFactor -1;
+    vtkOutput true;
+}
\ No newline at end of file

apps/RobinLaplace/helpers/MicroEdgeIterator.hpp

+//Add header
+
+#pragma once
+
+#include "core/debug/Debug.h"
+#include "hyteg/indexing/Common.hpp"
+
+namespace hyteg {
+namespace indexing {
+
+class MicroEdgeIterator
+{
+public:
+
+  using iterator_category = std::input_iterator_tag;
+  using value_type        = std::array< Index, 2 >;
+  using reference         = value_type const&;
+  using pointer           = value_type const*;
+  using difference_type   = ptrdiff_t;
+
+  MicroEdgeIterator( const uint_t & width, 
+                     const uint_t & offsetToCenter = 0, 
+                     const bool &   backwards = false, 
+                     const bool &   end = false ) 
+    : width_( width )
+    , offsetToCenter_( offsetToCenter )
+    , backwards_( backwards )
+    , totalNumberOfDoFs_( width - 2 * offsetToCenter )
+    , step_( 0 )
+  {
+    WALBERLA_ASSERT_GREATER( width, 0, "Size of edge must be larger than zero!" );
+    WALBERLA_ASSERT_LESS( offsetToCenter - 1, width, "Offset to center is beyond edge width!" );
+
+    coordinates_[0].dep() = 0;
+    coordinates_[0].row() = 0;
+    coordinates_[1].dep() = 0;
+    coordinates_[1].row() = 0;
+
+    if ( backwards )
+    {
+      coordinates_[0].col() = width - 1 - offsetToCenter;
+      coordinates_[1].col() = width - 2 - offsetToCenter;
+    }
+    else
+    {
+      coordinates_[0].col() = offsetToCenter;
+      coordinates_[1].col() = offsetToCenter + 1;
+    }
+
+    if ( end )
+    {
+      step_ = totalNumberOfDoFs_;
+    }
+  }
+
+  MicroEdgeIterator begin() { return MicroEdgeIterator( width_, offsetToCenter_, backwards_ ); }
+  MicroEdgeIterator end()   { return MicroEdgeIterator( width_, offsetToCenter_, backwards_, true ); }
+
+  bool operator==( const MicroEdgeIterator & other ) const
+  {
+    return other.step_ == step_;
+  }
+
+  bool operator!=( const MicroEdgeIterator & other ) const
+  {
+    return other.step_ != step_;
+  }
+
+  reference operator*()  const { return  coordinates_; };
+  pointer   operator->() const { return &coordinates_; };
+
+  MicroEdgeIterator & operator++() // prefix
+  {
+    WALBERLA_ASSERT_LESS_EQUAL( step_, totalNumberOfDoFs_, "Incrementing iterator beyond end!" );
+
+    step_++;
+    if ( backwards_ )
+    {
+      coordinates_[0].col()--;
+      coordinates_[1].col()--;
+    }
+    else
+    {
+      coordinates_[0].col()++;
+      coordinates_[1].col()++;
+    }
+
+    return *this;
+  }
+
+  MicroEdgeIterator operator++( int ) // postfix
+  {
+    const MicroEdgeIterator tmp( *this );
+    ++*this;
+    return tmp;
+  }
+
+
+private:
+
+  uint_t                                    width_;
+  uint_t                                    offsetToCenter_;
+  bool                                      backwards_;
+  uint_t                                    totalNumberOfDoFs_;
+  uint_t                                    step_;
+  std::array< hyteg::indexing::Index, 2 >   coordinates_;
+
+};
+
+}
+}
\ No newline at end of file

apps/RobinLaplace/helpers/P1CustomBoundaryOperator.cpp

+#include "P1CustomBoundaryOperator.hpp"
+
+namespace hyteg {
+
+using walberla::real_t;
+
+P1CustomBoundaryOperator::P1CustomBoundaryOperator(   const std::shared_ptr< PrimitiveStorage >&        storage,
+                                                      size_t                                            minLevel,
+                                                      size_t                                            maxLevel,
+                                                      real_t                                            robinFactor,
+                                                      std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction )
+: Operator( storage, minLevel, maxLevel )
+, robinFactor_( robinFactor )
+, boundaryFunction_( boundaryFunction )
+{
+   quadratureNodes_.push_back( 0.5 );
+   quadratureWeights_.push_back( 1.0 );
+}
+
+P1CustomBoundaryOperator::P1CustomBoundaryOperator(   const std::shared_ptr< PrimitiveStorage >&        storage,
+                                                      size_t                                            minLevel,
+                                                      size_t                                            maxLevel,
+                                                      real_t                                            robinFactor,
+                                                      std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
+                                                      std::vector< real_t >                             quadratureNodes,
+                                                      std::vector< real_t >                             quadratureWeights )
+: Operator( storage, minLevel, maxLevel )
+, robinFactor_( robinFactor )
+, boundaryFunction_( boundaryFunction )
+, quadratureNodes_( quadratureNodes )
+, quadratureWeights_( quadratureWeights )
+{}
+
+void P1CustomBoundaryOperator::apply(  const P1Function< real_t >& src,
+                                       const P1Function< real_t >& dst,
+                                       size_t                      level,
+                                       DoFType                     flag,
+                                       UpdateType                  updateType ) const
+{
+   if ( this->storage_->hasGlobalCells() )
+   {
+      WALBERLA_ABORT( "Operator with custom boundary conditions is not available for global cells yet" );
+   }
+   else
+   {
+      std::cout << "\n\nBoundary operator was called..." << std::endl;
+      communication::syncFunctionBetweenPrimitives( src, level );
+
+      std::vector< hyteg::PrimitiveID > edgeIDs = this->getStorage()->getEdgeIDs();
+      for ( int i = 0; i < int_c( edgeIDs.size() ); i++ )
+      {
+         hyteg::Edge& edge = *this->getStorage()->getEdge( edgeIDs[uint_c( i )] );
+
+         const DoFType edgeBC = dst.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
+         if ( testFlag( edgeBC, flag ) )
+         {   
+            for( const auto& microEdgeIndices: hyteg::indexing::MicroEdgeIterator( levelinfo::num_microedges_per_edge( level ) ))
+            {  
+               quadratureOnBoundaryEdge( level, edge, microEdgeIndices, src, dst, updateType );
+            }
+         }
+      }
+      dst.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage_, false );
+      dst.communicate< Vertex, Edge >( level );
+      dst.communicate< Edge, Face >( level ); 
+   }
+}
+
+void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t                                 level,
+                                                         const Edge&                            edge,
+                                                         const std::array< indexing::Index, 2 > microedgeIdx,
+                                                         const P1Function< real_t >&            src,
+                                                         const P1Function< real_t >&            dst,
+                                                         UpdateType                             update ) const
+{  
+   std::cout << "\nquadratureOnBoundaryEdge() was called..." << std::endl;
+   real_t edgeNum = real_t( levelinfo::num_microedges_per_edge( level ) ); 
+   Point3D microEdgeLength = (edge.getCoordinates()[1] - edge.getCoordinates()[0]) / edgeNum;
+   Point3D coord0 = edge.getCoordinates()[0] + real_t( microedgeIdx[0].col() ) * microEdgeLength;
+   Point3D coord1 = edge.getCoordinates()[0] + real_t( microedgeIdx[1].col() ) * microEdgeLength;
+
+   std::cout << "Quadrature on boundary with parameters..." << std::endl;
+   std::cout << "     edge.getCoordinates()[0] = " << edge.getCoordinates()[0] << std::endl;
+   std::cout << "     edge.getCoordinates()[1] = " << edge.getCoordinates()[1] << std::endl;
+   std::cout << "     microedgeIdx[0].col()    = " << microedgeIdx[0].col() << std::endl;
+   std::cout << "     microedgeIdx[1].col()    = " << microedgeIdx[1].col() << std::endl;
+   std::cout << "     edgeNum                  = " << edgeNum << std::endl;
+   std::cout << "     coord0                   = " << coord0 << std::endl;
+   std::cout << "     coord1                   = " << coord1 << std::endl;
+   std::cout << "     microEdgeLength          = " << microEdgeLength << std::endl; 
+
+   real_t integralUU = 0;
+   real_t integralUV = 0;
+   real_t integralVU = 0;
+   real_t integralVV = 0;
+   real_t integralMicro0 = 0;
+   real_t integralMicro1 = 0;
+
+   for(uint_t i = 0; i < quadratureNodes_.size(); i++) 
+   {
+      Point3D currentNodeCoord({ coord0[0] * quadratureNodes_[i] + coord1[0] * (1 - quadratureNodes_[i]),
+                                 coord0[1] * quadratureNodes_[i] + coord1[1] * (1 - quadratureNodes_[i])});
+      real_t boundaryValue = boundaryFunction_( currentNodeCoord );
+
+      integralUU += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * quadratureNodes_[i] * quadratureNodes_[i];
+      integralUV += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * quadratureNodes_[i] * (1 - quadratureNodes_[i]);
+      integralVU += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * (1 - quadratureNodes_[i]) * quadratureNodes_[i];
+      integralVV += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * (1 - quadratureNodes_[i]) * (1 - quadratureNodes_[i]);
+      integralMicro0 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[i];
+      integralMicro1 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[i]);
+   }
+   
+   std::cout << "Boundary operator came up with the following integrals:" << std::endl;
+   std::cout << "     integralUU     = " << integralUU << std::endl;
+   std::cout << "     integralUV     = " << integralUV << std::endl;
+   std::cout << "     integralVU     = " << integralVU << std::endl;
+   std::cout << "     integralVV     = " << integralVV << std::endl;
+   std::cout << "     integralMicro0 = " << integralMicro0 << std::endl;
+   std::cout << "     integralMicro1 = " << integralMicro1 << std::endl;  
+
+   PrimitiveDataID< FunctionMemory< real_t >, Edge > dstEdgeDoFIdx = dst.getEdgeDataID();
+   PrimitiveDataID< FunctionMemory< real_t >, Edge > srcEdgeDoFIdx = src.getEdgeDataID();
+   auto dstData = edge.getData( dstEdgeDoFIdx )->getPointer( level );
+   auto srcData = edge.getData( srcEdgeDoFIdx )->getPointer( level );
+
+   std::cout << "OLD: dstData[" << microedgeIdx[0].col() << "] = " << dstData[microedgeIdx[0].col()] << ", dstData[" << microedgeIdx[1].col() << "] = " << dstData[microedgeIdx[1].col()] << std::endl;
+
+   if (update == Replace) 
+   {
+      if( microedgeIdx[0].col() == 0 ) 
+      {
+         dstData[ microedgeIdx[0].col() ] = integralUU * srcData[ microedgeIdx[0].col() ] + integralUV * srcData[ microedgeIdx[1].col() ];
+      } else  
+      {
+         dstData[ microedgeIdx[0].col() ] += integralUU * srcData[ microedgeIdx[0].col() ] + integralUV * srcData[ microedgeIdx[1].col() ];
+      }
+      dstData[ microedgeIdx[1].col() ] = integralVU * srcData[ microedgeIdx[0].col() ] + integralVV * srcData[ microedgeIdx[1].col() ];
+   } else if (update == Add) 
+   {
+      dstData[ microedgeIdx[0].col() ] += integralUU * srcData[ microedgeIdx[0].col() ] + integralUV * srcData[ microedgeIdx[1].col() ];
+      dstData[ microedgeIdx[1].col() ] += integralVU * srcData[ microedgeIdx[0].col() ] + integralVV * srcData[ microedgeIdx[1].col() ];
+   } 
+   dstData[ microedgeIdx[0].col() ] += integralMicro0;
+   dstData[ microedgeIdx[1].col() ] += integralMicro1;
+   std::cout << "NEW: dstData[" << microedgeIdx[0].col() << "] = " << dstData[microedgeIdx[0].col()] << ", dstData[" << microedgeIdx[1].col() << "] = " << dstData[microedgeIdx[1].col()] << std::endl;
+   
+}
+
+} // namespace hyteg
+
+

apps/RobinLaplace/helpers/P1CustomBoundaryOperator.hpp

+#pragma once
+
+#include "hyteg/operators/Operator.hpp"
+#include "hyteg/p1functionspace/P1Function.hpp"
+#include "hyteg/p1functionspace/P1Operator.hpp"
+#include "hyteg/forms/P1Form.hpp"
+#include "hyteg/communication/Syncing.hpp"
+
+#include "MicroEdgeIterator.hpp"
+
+// #include "P1CustomBoundaryOperator.cpp"
+
+namespace hyteg {
+
+using walberla::real_t;
+
+class P1CustomBoundaryOperator : public Operator< P1Function< real_t >, P1Function< real_t > >
+{
+public:
+   P1CustomBoundaryOperator(  const std::shared_ptr< hyteg::PrimitiveStorage >&       storage,
+                              size_t                                                  minLevel,
+                              size_t                                                  maxLevel,
+                              real_t                                                  robinFactor,
+                              std::function< real_t( const hyteg::Point3D& ) >&       boundaryFunction );
+
+   P1CustomBoundaryOperator(  const std::shared_ptr< hyteg::PrimitiveStorage >&       storage,
+                              size_t                                                  minLevel,
+                              size_t                                                  maxLevel,
+                              real_t                                                  robinFactor,
+                              std::function< real_t( const hyteg::Point3D& )>&        boundaryFunction,
+                              std::vector< real_t >                                   quadratureNodes,
+                              std::vector< real_t >                                   quadratureWeights );
+
+   virtual ~P1CustomBoundaryOperator() = default;
+
+   void apply( const P1Function< real_t >& src,
+               const P1Function< real_t >& dst,
+               size_t                      level,
+               DoFType                     flag,
+               UpdateType                  updateType = Replace ) const override final;
+
+   void quadratureOnBoundaryEdge(   uint_t                                 level,
+                                    const Edge&                            edge,
+                                    const std::array< indexing::Index, 2 > microedgeIdx, 
+                                    const P1Function< real_t >&            src,
+                                    const P1Function< real_t >&            dst,
+                                    UpdateType                             update ) const;
+
+private:
+    real_t robinFactor_;
+    std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction_;
+    std::vector< real_t > quadratureNodes_;
+    std::vector< real_t > quadratureWeights_;
+};
+
+}
\ No newline at end of file

apps/RobinLaplace/helpers/P1OperatorWithBC.cpp

+#include "P1OperatorWithBC.hpp"
+
+namespace hyteg {
+
+using walberla::real_t;
+