wip adds normal (wrong) normal field for blended geometry

src/hyteg/elementwiseoperators/P1ToP2ElementwiseOperator.cpp

@@ -573,6 +573,7 @@ template class P1ToP2ElementwiseOperator<
 
 template class P1ToP2ElementwiseOperator< P1ToP2Form_divt< 0 > >;
 template class P1ToP2ElementwiseOperator< P1ToP2Form_divt< 1 > >;
+template class P1ToP2ElementwiseOperator< P1ToP2Form_divt< 2 > >;
 
 
 } // namespace hyteg

src/hyteg/elementwiseoperators/P1ToP2ElementwiseOperator.hpp

@@ -142,6 +142,7 @@ typedef P1ToP2ElementwiseOperator< P1ToP2FenicsForm< fenics::NoAssemble, p1_to_p
 
 typedef P1ToP2ElementwiseOperator< P1ToP2Form_divt< 0 > > P1ToP2ElementwiseBlendingDivTxOperator;
 typedef P1ToP2ElementwiseOperator< P1ToP2Form_divt< 1 > > P1ToP2ElementwiseBlendingDivTyOperator;
+typedef P1ToP2ElementwiseOperator< P1ToP2Form_divt< 2 > > P1ToP2ElementwiseBlendingDivTzOperator;
 
 
 } // namespace hyteg

src/hyteg/elementwiseoperators/P2NormalFieldWithBlendingEstimator.hpp

+/*
+ * Copyright (c) 2017-2020 Daniel Drzisga, Dominik Thoennes, Nils Kohl.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+#pragma once
+
+#include "hyteg/composites/P1StokesFunction.hpp"
+#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
+#include "hyteg/composites/P2P1TaylorHoodStokesBlockPreconditioner.hpp"
+#include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
+#include "hyteg/elementwiseoperators/P1ToP2ElementwiseOperator.hpp"
+#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
+#include "hyteg/elementwiseoperators/P2ToP1ElementwiseOperator.hpp"
+
+namespace hyteg {
+
+/// Estimates the normal field on the boundary using the divergence theorem .
+///
+/// Basic idea:
+///     TODO: Describe how it works
+///
+class P1NormalFieldWithBlendingEstimator
+{
+ public:
+   P1NormalFieldWithBlendingEstimator( const std::shared_ptr< PrimitiveStorage >& storage,
+                                       size_t                                     minLevel,
+                                       size_t                                     maxLevel )
+   : divT_x( storage, minLevel, maxLevel )
+   , divT_y( storage, minLevel, maxLevel )
+   , divT_z( storage, minLevel, maxLevel )
+   , hasGlobalCells_( storage->hasGlobalCells() )
+   {}
+
+   void apply( const P1StokesFunction< real_t >& dst,
+               const P1StokesFunction< real_t >& tmp,
+               const uint_t                            level ) const
+   {
+      // initialize the pressure to zero
+      tmp.p.interpolate( 1, level, All );
+
+      dst.u.interpolate(0, level, All);
+      dst.v.interpolate(0, level, All);
+      dst.w.interpolate(0, level, All);
+
+      // Assemble int 1 * d\phi_j/dx_i dx
+      divT_x.apply( tmp.p, dst.u, level, All, Replace );
+      divT_y.apply( tmp.p, dst.v, level, All, Replace );
+
+      if( hasGlobalCells_ )
+         divT_z.apply( tmp.p, dst.w, level, All, Replace );
+
+      // Normalization procedure:
+      // Step 1: calculate
+      //    tmp.u = dst.u^2
+      //    tmp.v = dst.v^2
+      tmp.u.multElementwise({ dst.u, dst.u }, level, All);
+      tmp.v.multElementwise({ dst.v, dst.v }, level, All);
+
+      if( hasGlobalCells_ )
+         tmp.w.multElementwise({ dst.w, dst.w }, level, All);
+
+      // tmp.u = dst.u^2 + dst.v^2
+      if ( hasGlobalCells_ )
+         tmp.u.assign({1, 1, 1}, {tmp.u, tmp.v, tmp.w}, level, All);
+      else
+         tmp.u.assign({1, 1}, {tmp.u, tmp.v}, level, All);
+
+      // Step 2: functional for calculating
+      //    dst.u = dst.u / sqrt( dst.u^2 + dst.v^2 )
+      //    ...
+      const auto normalizationFunctional = [](auto, auto val) {
+        auto norm = sqrt(val[1]);
+        WALBERLA_LOG_INFO("norm" << norm );
+        return norm > 1e-4 ? -val[0]/norm: 0;
+      };
+
+      dst.u.interpolateExtended(normalizationFunctional, {dst.u, tmp.u}, level, All);
+      dst.v.interpolateExtended(normalizationFunctional, {dst.v, tmp.u}, level, All);
+
+      if( hasGlobalCells_ )
+         dst.w.interpolateExtended(normalizationFunctional, {dst.w, tmp.u }, level, All);
+
+      // dst.interpolate(0, level, Inner);
+   }
+
+ private:
+   /*
+   P1ToP2ElementwiseBlendingDivTxOperator divT_x;
+   P1ToP2ElementwiseBlendingDivTyOperator divT_y;
+   */
+
+   /*
+   P1ToP2ElementwiseDivTxOperator divT_x;
+   P1ToP2ElementwiseDivTyOperator divT_y;
+   P1ToP2ElementwiseDivTzOperator divT_z;
+    */
+
+   /*
+   P1ToP2ConstantDivTxOperator divT_x;
+   P1ToP2ConstantDivTyOperator divT_y;
+   P1ToP2ConstantDivTzOperator divT_z;
+    */
+
+   P1DivTxOperator divT_x;
+   P1DivTyOperator divT_y;
+   P1DivTzOperator divT_z;
+
+   bool hasGlobalCells_;
+};
+
+
+class P2NormalFieldWithBlendingEstimator
+{
+ public:
+   P2NormalFieldWithBlendingEstimator( const std::shared_ptr< PrimitiveStorage >& storage,
+                                       size_t                                     minLevel,
+                                       size_t                                     maxLevel )
+       : divT_x( storage, minLevel, maxLevel )
+       , divT_y( storage, minLevel, maxLevel )
+       , divT_z( storage, minLevel, maxLevel )
+       , hasGlobalCells_( storage->hasGlobalCells() )
+   {}
+
+   void apply( const P2P1TaylorHoodFunction< real_t >& dst,
+                 const P2P1TaylorHoodFunction< real_t >& tmp,
+                 const uint_t                            level ) const
+   {
+      // initialize the pressure to zero
+      tmp.p.interpolate( 1e13, level, All );
+
+      dst.u.interpolate(0, level, All);
+      dst.v.interpolate(0, level, All);
+      dst.w.interpolate(0, level, All);
+
+      // Assemble int 1 * d\phi_j/dx_i dx
+      //divT_x.getVertexToVertexOpr().apply( tmp.p, dst.u.getVertexDoFFunction(), level, All, Replace );
+      //divT_y.getVertexToVertexOpr().apply( tmp.p, dst.v.getVertexDoFFunction(), level, All, Replace );
+      divT_x.apply( tmp.p, dst.u, level, All, Replace );
+      divT_y.apply( tmp.p, dst.v, level, All, Replace );
+
+      if( hasGlobalCells_ )
+         divT_z.apply( tmp.p, dst.w, level, All, Replace );
+      //divT_z.getVertexToVertexOpr().apply( tmp.p, dst.w.getVertexDoFFunction(), level, All, Replace );
+
+      dst.u.communicate<Vertex, Edge>(level);
+      dst.u.communicate<Edge, Face>(level);
+      dst.u.communicate<Face, Cell>(level);
+
+      // Normalization procedure:
+      // Step 1: calculate
+      //    tmp.u = dst.u^2
+      //    tmp.v = dst.v^2
+      tmp.u.multElementwise({ dst.u, dst.u }, level, All);
+      tmp.v.multElementwise({ dst.v, dst.v }, level, All);
+
+      if( hasGlobalCells_ )
+         tmp.w.multElementwise({ dst.w, dst.w }, level, All);
+
+      // tmp.u = dst.u^2 + dst.v^2
+      if ( hasGlobalCells_ )
+         tmp.u.assign({1, 1, 1}, {tmp.u, tmp.v, tmp.w}, level, All);
+      else
+         tmp.u.assign({1, 1}, {tmp.u, tmp.v}, level, All);
+
+      // Step 2: functional for calculating
+      //    dst.u = dst.u / sqrt( dst.u^2 + dst.v^2 )
+      //    ...
+      const auto normalizationFunctional = [](auto, auto val) {
+         auto norm = sqrt(val[1]);
+        return norm > 1e-12 ? -val[0]/norm: 0;
+      };
+
+      dst.u.interpolate(normalizationFunctional, {dst.u, tmp.u}, level, All);
+      dst.v.interpolate(normalizationFunctional, {dst.v, tmp.u}, level, All);
+
+      if( hasGlobalCells_ )
+         dst.w.interpolate(normalizationFunctional, {dst.w, tmp.u }, level, All);
+
+      communication::syncP2FunctionBetweenPrimitives( dst.u, level );
+      communication::syncP2FunctionBetweenPrimitives( dst.v, level );
+      communication::syncP2FunctionBetweenPrimitives( dst.w, level );
+
+      // dst.interpolate(0, level, Inner);
+   }
+
+ private:
+   P1ToP2ElementwiseBlendingDivTxOperator divT_x;
+   P1ToP2ElementwiseBlendingDivTyOperator divT_y;
+   P1ToP2ElementwiseBlendingDivTzOperator divT_z;
+
+   /*
+   P1ToP2ElementwiseDivTxOperator divT_x;
+   P1ToP2ElementwiseDivTyOperator divT_y;
+   P1ToP2ElementwiseDivTzOperator divT_z;
+    */
+
+   /*
+   P1ToP2ConstantDivTxOperator divT_x;
+   P1ToP2ConstantDivTyOperator divT_y;
+   P1ToP2ConstantDivTzOperator divT_z;
+    */
+
+   bool hasGlobalCells_;
+};
+
+} // namespace hyteg
\ No newline at end of file

tests/hyteg/CMakeLists.txt

@@ -272,6 +272,8 @@ waLBerla_compile_test(FILES edgedofspace/EdgeDoFFunction3DTest.cpp DEPENDS hyteg
 waLBerla_execute_test(NAME EdgeDoFFunction3DTest)
 
 ## Operators ##
+waLBerla_compile_test(FILES operators/P2NormalFieldWithBlendingEstimator.cpp DEPENDS hyteg core)
+waLBerla_execute_test(NAME P2NormalFieldWithBlendingEstimator)
 
 waLBerla_compile_test(FILES operators/EdgeDoFToVertexDoFOperatorTest.cpp DEPENDS hyteg core)
 waLBerla_execute_test(NAME EdgeDoFToVertexDoFOperatorTest)

tests/hyteg/operators/P2NormalFieldWithBlendingEstimator.cpp

+/*
+ * Copyright (c) 2017-2019 Dominik Bartuschat, Dominik Thoennes, Marcus Mohr, Nils Kohl.
+ *
+ * This file is part of HyTeG
+ * (see https://i10git.cs.fau.de/hyteg/hyteg).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+#include <cmath>
+#include <hyteg/elementwiseoperators/P2NormalFieldWithBlendingEstimator.hpp>
+#include <hyteg/geometry/AnnulusMap.hpp>
+#include <hyteg/geometry/IcosahedralShellMap.hpp>
+#include <hyteg/primitivestorage/SetupPrimitiveStorage.hpp>
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/mpi/MPIManager.h"
+
+#include "hyteg/composites/P1StokesFunction.hpp"
+#include "hyteg/dataexport/VTKOutput.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/CGSolver.hpp"
+#include "hyteg/solvers/GaussSeidelSmoother.hpp"
+#include "hyteg/solvers/GeometricMultigridSolver.hpp"
+
+using walberla::real_t;
+using walberla::uint_c;
+using walberla::uint_t;
+
+using namespace hyteg;
+
+/*
+int main( int argc, char** argv )
+{
+   walberla::Environment env( argc, argv );
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   const uint_t level = 2;
+
+   const uint_t gridSize = 1;
+   MeshInfo meshInfo = MeshInfo::meshCuboid( Point3D({-1, -1, -1}), Point3D({ 1, 1, 1}), gridSize, gridSize, gridSize );
+
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, true );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   P2P1TaylorHoodFunction< real_t > u1( "u", storage, level, level );
+   P2P1TaylorHoodFunction< real_t > u2( "u", storage, level, level );
+   P2P1TaylorHoodFunction< real_t > tmp( "tmp", storage, level, level );
+
+   P1ToP2ElementwiseDivTxOperator divT_xe(storage, level, level);
+   P1ToP2ElementwiseDivTyOperator divT_ye(storage, level, level);
+   P1ToP2ElementwiseDivTzOperator divT_ze(storage, level, level);
+
+   P1ToP2ConstantDivTxOperator divT_x(storage, level, level);
+   P1ToP2ConstantDivTyOperator divT_y(storage, level, level);
+   P1ToP2ConstantDivTzOperator divT_z(storage, level, level);
+
+   u1.p.interpolate(1, level, All);
+   u2.p.interpolate(1, level, All);
+
+   u1.p.interpolate([](auto v) { return v[0] + v[1] + v[2]; }, level, All);
+   u2.p.interpolate([](auto v) { return v[0] + v[1] + v[2]; }, level, All);
+
+   divT_xe.apply(u1.p, u1.u, level, All, Replace);
+   divT_ye.apply(u1.p, u1.v, level, All, Replace);
+   divT_ze.apply(u1.p, u1.w, level, All, Replace);
+
+   divT_x.apply(u2.p, u2.u, level, All, Replace);
+   divT_y.apply(u2.p, u2.v, level, All, Replace);
+   divT_z.apply(u2.p, u2.w, level, All, Replace);
+
+   tmp.u.assign({1, -1}, {u1.u, u2.u}, level, All);
+   WALBERLA_LOG_INFO("u error inner " << tmp.u.sumGlobal(level, Inner, true));
+   WALBERLA_LOG_INFO("u error all " << tmp.u.sumGlobal(level, Inner, false));
+
+   tmp.v.assign({1, -1}, {u1.v, u2.v}, level, All);
+   WALBERLA_LOG_INFO("v error inner " << tmp.v.sumGlobal(level, Inner, true));
+   WALBERLA_LOG_INFO("v error all " << tmp.v.sumGlobal(level, Inner, false));
+
+   tmp.w.assign({1, -1}, {u1.w, u2.w}, level, All);
+   WALBERLA_LOG_INFO("w error inner " << tmp.w.sumGlobal(level, Inner, true));
+   WALBERLA_LOG_INFO("w error all " << tmp.w.sumGlobal(level, Inner, false));
+}
+ */
+
+/*
+int main( int argc, char** argv )
+{
+   walberla::Environment env( argc, argv );
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   const uint_t level = 2;
+
+   const uint_t gridSize = 1;
+   // MeshInfo meshInfo = MeshInfo::meshRectangle( Point2D( { -1, -1 } ), Point2D( { 1., 1. } ), MeshInfo::CRISS, gridSize, gridSize );
+   MeshInfo meshInfo = MeshInfo::meshCuboid( Point3D({-1, -1, -1}), Point3D({ 1, 1, 1}), gridSize, gridSize, gridSize );
+
+   // MeshInfo meshInfo = MeshInfo::meshAnnulus(0.5, 1.0, MeshInfo::CRISS, 5, 10);
+
+   // MeshInfo meshInfo = MeshInfo::meshSphericalShell(5, 10, 0.5, 1.0);
+
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, true );
+   //AnnulusMap::setMap( setupStorage );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   // P2P1TaylorHoodFunction< real_t > u( "u", storage, level, level );
+   // P2P1TaylorHoodFunction< real_t > tmp( "tmp", storage, level, level );
+   P1StokesFunction< real_t > u( "u", storage, level, level );
+   P1StokesFunction< real_t > tmp( "tmp", storage, level, level );
+
+   P1NormalFieldWithBlendingEstimator normalEstimator ( storage, level, level );
+
+   tmp.interpolate(0, level, All);
+   u.interpolate(0, level, All);
+
+   normalEstimator.apply( u, tmp, level );
+
+   hyteg::VTKOutput vtkOutput( ".", "TestNormal", storage );
+   vtkOutput.add( u );
+   vtkOutput.write( level );
+}
+ */
+
+
+int main( int argc, char** argv )
+{
+   walberla::Environment env( argc, argv );
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   const uint_t level = 2;
+
+   const uint_t gridSize = 1;
+   // MeshInfo meshInfo = MeshInfo::meshRectangle( Point2D( { -1, -1 } ), Point2D( { 1., 1. } ), MeshInfo::CRISS, gridSize, gridSize );
+   MeshInfo meshInfo = MeshInfo::meshCuboid( Point3D({-1, -1, -1}), Point3D({ 1, 1, 1}), gridSize, gridSize, gridSize );
+
+   // MeshInfo meshInfo = MeshInfo::meshAnnulus(0.5, 1.0, MeshInfo::CRISS, 5, 10);
+
+   // MeshInfo meshInfo = MeshInfo::meshSphericalShell(5, 10, 0.5, 1.0);
+
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   //setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, true );
+   // AnnulusMap::setMap( setupStorage );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   P2P1TaylorHoodFunction< real_t > u( "u", storage, level, level );
+   P2P1TaylorHoodFunction< real_t > tmp( "tmp", storage, level, level );
+
+   P2NormalFieldWithBlendingEstimator normalEstimator ( storage, level, level );
+
+   tmp.interpolate(0, level, All);
+   u.interpolate(0, level, All);
+
+   normalEstimator.apply( u, tmp, level );
+
+   hyteg::VTKOutput vtkOutput( ".", "NormalField", storage );
+   vtkOutput.add( u );
+   vtkOutput.write( level );
+
+   // hyteg::VTKOutput vtkOutput( ".", "TestNormalTwo", storage );
+   // vtkOutput.add( u.u.getVertexDoFFunction() );
+   // vtkOutput.add( u.v.getVertexDoFFunction() );
+   // vtkOutput.add( u.w.getVertexDoFFunction() );
+   // vtkOutput.write( level );
+}
+
+/*
+int main( int argc, char** argv )
+{
+   walberla::Environment env( argc, argv );
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   const uint_t level = 2;
+
+   const uint_t gridSize = 1;
+   // MeshInfo meshInfo = MeshInfo::meshRectangle( Point2D( { -1, -1 } ), Point2D( { 1., 1. } ), MeshInfo::CRISS, gridSize, gridSize );
+   MeshInfo meshInfo = MeshInfo::meshCuboid( Point3D({-1, -1, -1}), Point3D({ 1, 1, 1}), gridSize, gridSize, gridSize );
+
+   // MeshInfo meshInfo = MeshInfo::meshAnnulus(0.5, 1.0, MeshInfo::CRISS, 5, 10);
+
+   // MeshInfo meshInfo = MeshInfo::meshSphericalShell(5, 10, 0.5, 1.0);
+
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, true );
+   // AnnulusMap::setMap( setupStorage );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   P2P1TaylorHoodFunction< real_t > u1( "u", storage, level, level );
+   u1.interpolate(0, level, All);
+   u1.u.interpolate([](auto p){ return p[0]*p[1]; }, level, All); // -> y
+   u1.v.interpolate([](auto p){ return p[0]*p[2]; }, level, All); // -> 0
+   u1.w.interpolate([](auto p){ return p[1]*p[2]; }, level, All); // -> y
+
+   P2P1TaylorHoodFunction< real_t > u2( "u", storage, level, level );
+   u2.interpolate(0, level, All);
+   u2.p.interpolate([](auto p){ return p[1];}, level, All);
+
+   P1ToP2ConstantDivTxOperator divT_x(storage, level, level);
+   P1ToP2ConstantDivTyOperator divT_y(storage, level, level);
+   P1ToP2ConstantDivTzOperator divT_z(storage, level, level);
+
+   divT_x.apply(u2.p, u2.u, level, All);
+   divT_y.apply(u2.p, u2.v, level, All);
+   divT_z.apply(u2.p, u2.w, level, All);
+
+   auto dot = u2.u.dotGlobal(u1.u, level) + u2.v.dotGlobal(u1.v, level) + u2.w.dotGlobal(u1.w, level);
+   WALBERLA_LOG_INFO( u2.u.getMaxValue(level, All) );
+   WALBERLA_LOG_INFO( u1.u.getMaxValue(level, All) );
+   //u2.u.multElementwise({u2.u, u1.u}, level );
+   // WALBERLA_LOG_INFO( u2.u.sumGlobal(level, All) );
+   WALBERLA_LOG_INFO( u2.u.dotGlobal(u1.u, level) );
+   WALBERLA_LOG_INFO( u2.v.dotGlobal(u1.v, level) );
+   WALBERLA_LOG_INFO( u2.w.dotGlobal(u1.w, level) );
+
+   hyteg::VTKOutput vtkOutput( ".", "DoSomething", storage );
+   vtkOutput.add( u2.u );
+   vtkOutput.write( level );
+
+   WALBERLA_LOG_INFO("dot = " << dot);
+
+}
+*/
+
+/*
+int main( int argc, char** argv )
+{
+   walberla::Environment env( argc, argv );
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   const uint_t level = 2;
+
+   const uint_t gridSize = 1;
+   // MeshInfo meshInfo = MeshInfo::meshRectangle( Point2D( { -1, -1 } ), Point2D( { 1., 1. } ), MeshInfo::CRISS, gridSize, gridSize );
+   MeshInfo meshInfo = MeshInfo::meshCuboid( Point3D({-1, -1, -1}), Point3D({ 1, 1, 1}), gridSize, gridSize, gridSize );
+
+   // MeshInfo meshInfo = MeshInfo::meshAnnulus(0.5, 1.0, MeshInfo::CRISS, 5, 10);
+
+   // MeshInfo meshInfo = MeshInfo::meshSphericalShell(5, 10, 0.5, 1.0);
+
+   SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
+   setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true );
+   std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage );
+
+   P2Function< real_t > u1( "u", storage, level, level );
+   //u1.interpolate(1, level, All);
+   u1.interpolate([](auto p){ return p[0];}, level, All);
+   u1.communicate<Vertex, Edge>(level);
+   u1.communicate<Edge, Face>(level);
+   u1.communicate<Face, Cell>(level);
+   hyteg::VTKOutput vtkOutput( ".", "DoInterpolate", storage );
+   vtkOutput.add( u1 );
+   vtkOutput.write( level );
+}
+ */