diff --git a/apps/showcases/PorousMediaGPU/CMakeLists.txt b/apps/showcases/PorousMediaGPU/CMakeLists.txt
index a5cf4003da1439199eb85abfa17b0521426c702c..0d69de5d8ea5d04d949ff031d7490a0edeb103ab 100644
--- a/apps/showcases/PorousMediaGPU/CMakeLists.txt
+++ b/apps/showcases/PorousMediaGPU/CMakeLists.txt
@@ -18,11 +18,11 @@ waLBerla_generate_target_from_python(NAME PorousMediaGPUCodeGen
 
 if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
     waLBerla_add_executable ( NAME PorousMediaGPU
-            FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp
+            FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp TPMS.cpp
             DEPENDS blockforest boundary core field gpu lbm_generated postprocessing stencil timeloop vtk PorousMediaGPUCodeGen )
 else()
     waLBerla_add_executable ( NAME PorousMediaGPU
-            FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp
+            FILES PorousMediaGPU.cpp GridGeneration.h Setup.h Types.h InitializerFunctions.cpp TPMS.cpp
             DEPENDS blockforest boundary core field lbm_generated postprocessing stencil timeloop vtk PorousMediaGPUCodeGen )
 
 endif(WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
diff --git a/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp b/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
index d82111dc24e18f8ed53ff6218e5de108eb537255..0bff94fe9f298e734ec32bf5bbae1fe605cc78ce 100644
--- a/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
+++ b/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
@@ -55,7 +55,7 @@
 #   include "gpu/FieldCopy.h"
 #   include "gpu/HostFieldAllocator.h"
 #   include "gpu/ParallelStreams.h"
-#   include "gpu/communication/UniformGPUScheme.h"
+#   include "gpu/communication/NonUniformGPUScheme.h"
 #endif
 
 #include "lbm_generated/communication/UniformGeneratedPdfPackInfo.h"
@@ -66,7 +66,8 @@
 #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
 #   include "lbm_generated/gpu/AddToStorage.h"
 #   include "lbm_generated/gpu/GPUPdfField.h"
-#   include "lbm_generated/gpu/UniformGeneratedGPUPdfPackInfo.h"
+#   include "lbm_generated/gpu/NonuniformGeneratedGPUPdfPackInfo.h"
+#   include "lbm_generated/gpu/BasicRecursiveTimeStepGPU.h"
 #endif
 
 #include "postprocessing/FieldToSurfaceMesh.h"
@@ -81,6 +82,7 @@
 #include "InitializerFunctions.h"
 #include "GridGeneration.h"
 #include "Setup.h"
+#include "TPMS.h"
 #include "Types.h"
 
 #include "PorousMediaGPUStaticDefines.h"
@@ -93,9 +95,8 @@ using SweepCollection_T = lbm::PorousMediaGPUSweepCollection;
 
 #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
 using GPUPdfField_T = lbm_generated::GPUPdfField< StorageSpecification_T >;
-using gpu::communication::UniformGPUScheme;
-
-using lbm_generated::UniformGeneratedGPUPdfPackInfo;
+using gpu::communication::NonUniformGPUScheme;
+using lbm_generated::NonuniformGeneratedGPUPdfPackInfo;
 #else
 using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
 using blockforest::communication::UniformBufferedScheme;
@@ -117,11 +118,13 @@ int main(int argc, char** argv){
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
    auto parameters        = config->getOneBlock("Parameters");
+   auto tpmsParameters    = config->getOneBlock("TPMS");
    auto domainParameters  = config->getOneBlock("DomainSetup");
    auto loggingParameters = config->getOneBlock("Logging");
-   auto boundariesConfig   = config->getBlock("Boundaries");
+   auto boundariesConfig  = config->getBlock("Boundaries");
 
-   Setup setup(parameters, domainParameters, loggingParameters, infoMap);
+   Setup setup(parameters, tpmsParameters, domainParameters, loggingParameters, infoMap);
+   TPMS tpms(setup);
    WALBERLA_LOG_INFO_ON_ROOT( "Simulation run data:"
                               "\n- simulation parameters:"   << std::setprecision(16) <<
                               "\n   + collision model:     " << setup.collisionModel <<
@@ -176,25 +179,33 @@ int main(int argc, char** argv){
    for (auto& block : *blocks){
       sweepCollection.initialise(&block);
    }
+   sweepCollection.initialiseBlockPointer();
 
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///                                      LB SWEEPS AND BOUNDARY HANDLING                                       ///
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
    // Boundaries
-   if (boundariesConfig){
-      WALBERLA_LOG_INFO_ON_ROOT("Setting boundary conditions")
-      geometry::initBoundaryHandling< FlagField_T >(*blocks, ids.flagField, boundariesConfig);
-      init_TPMS(blocks, ids.flagField, setup.obstacleUID, setup.waveNumber, setup.scaling);
-      geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID);
+   WALBERLA_LOG_INFO_ON_ROOT("Setting boundary conditions")
+   for (auto& block : *blocks){
+      auto flagField = block.getData< FlagField_T >(ids.flagField);
+      setup.inflowFlag   = flagField->registerFlag(setup.inflowUID);
+      setup.outflowFlag  = flagField->registerFlag(setup.outflowUID);
+      setup.obstacleFlag = flagField->registerFlag(setup.obstacleUID);
+      setup.fluidFlag    = flagField->registerFlag(setup.fluidUID);
    }
+
+   geometry::initBoundaryHandling< FlagField_T >(*blocks, ids.flagField, boundariesConfig);
+   tpms.setupBoundary(blocks, ids.flagField);
+   // init_TPMS(blocks, ids.flagField, setup.obstacleUID, setup.waveNumber, setup.scaling);
+   // geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID);
    geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID);
    // BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, fluidFlagUID, latticeVelocity);
    BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfFieldGPU, setup.fluidUID, setup.latticeVelocity, ids.pdfField);
 
    if( loggingParameters.getParameter< bool >("writeSurfaceMeshOfTPMS", false) ){
       WALBERLA_LOG_INFO_ON_ROOT("Writing TPMS as surface file")
-      auto mesh = postprocessing::flagFieldToSurfaceMesh< FlagField_T >(blocks, ids.flagField, 4, true);
+      auto mesh = postprocessing::flagFieldToSurfaceMesh< FlagField_T >(blocks, ids.flagField, setup.obstacleFlag, true);
       WALBERLA_EXCLUSIVE_WORLD_SECTION(0){
          geometry::writeMesh("TPMS.obj", *mesh);
       }
@@ -204,20 +215,26 @@ int main(int argc, char** argv){
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
    const bool gpuEnabledMPI = parameters.getParameter< bool >("gpuEnabledMPI", false);
-   UniformGPUScheme< Stencil_T > communication(blocks, gpuEnabledMPI);
-   auto packInfo = std::make_shared<lbm_generated::UniformGeneratedGPUPdfPackInfo< GPUPdfField_T >>(ids.pdfFieldGPU);
-   communication.addPackInfo(packInfo);
+   auto nonUniformCommunication = std::make_shared< NonUniformGPUScheme< CommunicationStencil_T > >(blocks, gpuEnabledMPI);
+   auto nonUniformPackInfo      = lbm_generated::setupNonuniformGPUPdfCommunication< GPUPdfField_T >(blocks, ids.pdfFieldGPU);
+   nonUniformCommunication->addPackInfo(nonUniformPackInfo);
 
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///                                          TIME STEP DEFINITIONS                                             ///
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    const int streamLowPriority = 0;
    auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority);
+   SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
+
+   lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T >
+      LBMMeshRefinement(blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
 
-   SweepTimeloop timeLoop(blocks->getBlockStorage(), setup.timesteps);
-   timeLoop.add() << BeforeFunction(communication.getCommunicateFunctor(), "communication")
+   LBMMeshRefinement.addRefinementToTimeLoop(timeloop, uint_c(0));
+
+   /*SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
+   timeloop.add() << BeforeFunction(communication.getCommunicateFunctor(), "communication")
                   << Sweep(boundaryCollection.getSweep(BoundaryCollection_T::ALL, defaultStream), "Boundary Conditions");
-   timeLoop.add() << Sweep(sweepCollection.streamCollide(SweepCollection_T::ALL, defaultStream), "LBM StreamCollide");
+   timeloop.add() << Sweep(sweepCollection.streamCollide(SweepCollection_T::ALL, defaultStream), "LBM StreamCollide");*/
 
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    ///                                             VTK Output                                                     ///
@@ -291,17 +308,27 @@ int main(int argc, char** argv){
          auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag");
          vtkOutput->addCellDataWriter(flagWriter);
       }
-      timeLoop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+      timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
    }
 
-      // log remaining time
-   
+   // log remaining time
    if (uint_c(setup.remainingTimeLoggerFrequency) > 0){
-      timeLoop.addFuncAfterTimeStep(
-         timing::RemainingTimeLogger(timeLoop.getNrOfTimeSteps(), setup.remainingTimeLoggerFrequency),
+      timeloop.addFuncAfterTimeStep(
+         timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), setup.remainingTimeLoggerFrequency),
          "remaining time logger");
    }
 
+   auto CheckerParameters      = config->getOneBlock("StabilityChecker");
+   const uint_t checkFrequency = CheckerParameters.getParameter< uint_t >("checkFrequency", uint_t(0));
+   if (checkFrequency > 0)
+   {
+      auto checkFunction = [](VelocityField_T::value_type value) {  return value < math::abs(VelocityField_T::value_type(10)); };
+      timeloop.addFuncAfterTimeStep(
+         makeSharedFunctor(field::makeStabilityChecker< VelocityField_T, FlagField_T >(
+            config, blocks, ids.velocityField, ids.flagField, setup.fluidUID, checkFunction)),
+         "Stability check");
+   }
+
    WALBERLA_MPI_BARRIER()
 #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
    WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
@@ -323,7 +350,7 @@ int main(int argc, char** argv){
 #endif
 
    simTimer.start();
-   timeLoop.run(timeloopTiming);
+   timeloop.run(timeloopTiming);
 #if defined(WALBERLA_BUILD_WITH_GPU_SUPPORT)
    WALBERLA_GPU_CHECK(gpuDeviceSynchronize())
    WALBERLA_GPU_CHECK(gpuPeekAtLastError())
diff --git a/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm b/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
index 2f5c3d2270ee5833374d6895e4b28be4e6bc112a..c11fb59e415655d0ea8a8aafeeddfc983c91ad28 100644
--- a/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
+++ b/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
@@ -7,10 +7,7 @@ Parameters
     latticeVelocity 0.05;
     initialiseWithInletVelocity true;
 
-    waveNumber 6.2831853072;
-    scaling 10.0;
-
-    timesteps 10;
+    timesteps 10001;
 
     processMemoryLimit 512; // MiB
     innerOuterSplit <1, 1, 1>;
@@ -19,6 +16,16 @@ Parameters
     gpuEnabledMPI false;
 }
 
+TPMS
+{
+    waveNumber 6.2831853072;
+    scaling 10.0;
+    strut 1.0;
+    start 2.0;
+    end 6.0;
+    type gyroid;
+}
+
 //! [domainSetup]
 DomainSetup
 {
@@ -48,12 +55,12 @@ StabilityChecker
 
 VTKWriter
 {
-    vtkWriteFrequency 5;
+    vtkWriteFrequency 5000;
     vtkGhostLayers 0;
     velocity true;
     writeVelocityAsMagnitude true;
     density false;
-    flag false;
+    flag true;
     writeOnlySlice true;
     sliceThickness 1;
     writeXZSlice false;
@@ -69,7 +76,7 @@ Logging
     writeSetupForestAndReturn false; // When only one process is used the decomposition is writen and the program terminates
     readSetupFromFile false;
     remainingTimeLoggerFrequency 20; // in seconds
-    writeSurfaceMeshOfTPMS false;
+    writeSurfaceMeshOfTPMS true;
 }
 
 Evaluation
diff --git a/apps/showcases/PorousMediaGPU/PorousMediaGPU.py b/apps/showcases/PorousMediaGPU/PorousMediaGPU.py
index 6e1333dda7dd47254eda58de1a401971b8bd2219..17b08e1b6d8e36248a15b2faa16479ece900bb3c 100644
--- a/apps/showcases/PorousMediaGPU/PorousMediaGPU.py
+++ b/apps/showcases/PorousMediaGPU/PorousMediaGPU.py
@@ -38,7 +38,7 @@ with CodeGeneration() as ctx:
     q = stencil.Q
     dim = stencil.D
 
-    streaming_pattern = 'pull'
+    streaming_pattern = 'aa'
 
     pdfs = fields(f"pdfs({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
     velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
@@ -85,7 +85,7 @@ with CodeGeneration() as ctx:
 
     generate_lbm_package(ctx, name="PorousMediaGPU", collision_rule=collision_rule,
                          lbm_config=lbm_config, lbm_optimisation=lbm_opt,
-                         nonuniform=False, boundaries=[no_slip, ubb, outflow],
+                         nonuniform=True, boundaries=[no_slip, ubb, outflow],
                          macroscopic_fields=macroscopic_fields, gpu_indexing_params=sweep_params,
                          target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
                          max_threads=max_threads)
diff --git a/apps/showcases/PorousMediaGPU/Setup.h b/apps/showcases/PorousMediaGPU/Setup.h
index 889b8bd27a6a05888edbd68677edd670fd86f396..04ed737ae84b9388e99ecdd5d1f50ca145fc67fc 100644
--- a/apps/showcases/PorousMediaGPU/Setup.h
+++ b/apps/showcases/PorousMediaGPU/Setup.h
@@ -32,8 +32,8 @@ namespace walberla
 struct Setup
 {
 
-   Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & domainParameters,
-         const Config::BlockHandle & logging, std::map<std::string, std::string>& infoMap)
+   Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & tpms, 
+         const Config::BlockHandle & domainParameters, const Config::BlockHandle & logging, std::map<std::string, std::string>& infoMap)
    {
       blockForestFilestem = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest");
       numProcesses        = domainParameters.getParameter< uint_t >("numberProcesses");
@@ -53,8 +53,12 @@ struct Setup
       referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
       latticeVelocity   = parameters.getParameter< real_t >("latticeVelocity");
 
-      waveNumber = parameters.getParameter< real_t >("waveNumber");
-      scaling    = parameters.getParameter< real_t >("scaling");
+      waveNumber = tpms.getParameter< real_t >("waveNumber");
+      scaling    = tpms.getParameter< real_t >("scaling");
+      strut      = tpms.getParameter< real_t >("strut");
+      tpmsBeginn = tpms.getParameter< real_t >("start");
+      tpmsEnd    = tpms.getParameter< real_t >("end");
+      tpmsType   = tpms.getParameter< std::string >("type");
 
       const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
       machNumber = latticeVelocity / speedOfSound;
@@ -109,6 +113,10 @@ struct Setup
 
    real_t waveNumber;
    real_t scaling;
+   real_t strut;
+   real_t tpmsBeginn;
+   real_t tpmsEnd;
+   std::string tpmsType;
 
    real_t machNumber;
    real_t reynoldsNumber;
@@ -138,6 +146,11 @@ struct Setup
    FlagUID inflowUID;
    FlagUID outflowUID;
 
+   uint8_t inflowFlag;
+   uint8_t outflowFlag;
+   uint8_t obstacleFlag;
+   uint8_t fluidFlag;
+
    bool writeSetupForestAndReturn;
    real_t remainingTimeLoggerFrequency;
 };
diff --git a/apps/showcases/PorousMediaGPU/TPMS.cpp b/apps/showcases/PorousMediaGPU/TPMS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..82e7ddc3064d8404bd6dc697556487b66a0fb9b7
--- /dev/null
+++ b/apps/showcases/PorousMediaGPU/TPMS.cpp
@@ -0,0 +1,80 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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.
+//
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TPMS.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "TPMS.h"
+
+namespace walberla
+{
+
+real_t TPMS::getDistance(const real_t x, const real_t y, const real_t z, const real_t scaling) const {
+   const real_t kx = setup_.waveNumber * scaling;
+   const real_t ky = setup_.waveNumber * scaling;
+   const real_t kz = setup_.waveNumber * scaling;
+
+   if(tpmsType == 0){
+      return cos(kx * x) * sin(ky * y) + cos(ky * y) * sin(kz * z) + cos(kz * z) * sin(kx * x);
+   }
+   else{
+      return real_c(0.0);
+   }
+}
+
+bool TPMS::contains(const Vector3< real_t >& point) const
+{
+   const real_t dist1 = getDistance(point[0], point[1], point[2], real_c(1.0));
+   const real_t dist2 = getDistance(point[0], point[1], point[2], setup_.scaling);
+   return (dist1 * dist2) >= setup_.strut;
+}
+
+bool TPMS::contains(const AABB& aabb) const
+{
+   Vector3< real_t > p[8];
+   p[0].set(aabb.xMin(), aabb.yMin(), aabb.zMin());
+   p[1].set(aabb.xMax(), aabb.yMin(), aabb.zMin());
+   p[2].set(aabb.xMin(), aabb.yMax(), aabb.zMin());
+   p[3].set(aabb.xMax(), aabb.yMax(), aabb.zMin());
+   p[4].set(aabb.xMin(), aabb.yMin(), aabb.zMax());
+   p[5].set(aabb.xMax(), aabb.yMin(), aabb.zMax());
+   p[6].set(aabb.xMin(), aabb.yMax(), aabb.zMax());
+   p[7].set(aabb.xMax(), aabb.yMax(), aabb.zMax());
+   return contains(p[0]) && contains(p[1]) && contains(p[2]) && contains(p[3]) && contains(p[4]) && contains(p[5]) &&
+          contains(p[6]) && contains(p[7]);
+}
+
+void TPMS::setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID){
+   const real_t start = setup_.tpmsBeginn / setup_.dxC;
+   const real_t end   = setup_.tpmsEnd / setup_.dxC;
+   for (auto& block : *sbfs){
+      const auto flagField = block.template getData< FlagField_T >(flagFieldID);
+      const auto boundaryFlag = flagField->getOrRegisterFlag(setup_.obstacleUID);
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell;
+      sbfs->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+      if (globalCell[0] > start && globalCell[0] <= end){
+         Vector3<real_t> p;
+         sbfs->getCellCenter(p, globalCell, 0);
+         if(contains(p)){
+            addFlag(flagField->get(x, y, z), boundaryFlag);
+         }
+      }
+      )
+   }
+}
+
+} // namespace walberla
\ No newline at end of file
diff --git a/apps/showcases/PorousMediaGPU/TPMS.h b/apps/showcases/PorousMediaGPU/TPMS.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c470d6f343e27109879a4c5ba5aac5d3dbec582
--- /dev/null
+++ b/apps/showcases/PorousMediaGPU/TPMS.h
@@ -0,0 +1,103 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla 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.
+//
+//  waLBerla 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 waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file TPMS.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "blockforest/SetupBlock.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/StructuredBlockForest.h"
+
+#include "domain_decomposition/IBlock.h"
+
+#include "field/FlagUID.h"
+
+#include "core/DataTypes.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+#include "core/cell/Cell.h"
+
+#include "stencil/D3Q7.h"
+#include "stencil/D3Q27.h"
+
+#include "Types.h"
+#include "Setup.h"
+
+namespace walberla
+{
+
+class TPMS
+{
+ public:
+
+   TPMS(const Setup& setup) : setup_(setup){
+      if (setup.tpmsType.compare("gyroid") == 0){
+         tpmsType = 0;
+      }
+      else{
+         WALBERLA_ABORT("Unsoprted TPMS Type: " << setup.tpmsType);
+      }
+   }
+
+   real_t getDistance(const real_t x, const real_t y, const real_t z, const real_t scaling) const;
+   bool operator()(const Vector3< real_t >& point) const { return contains(point); }
+
+   bool contains(const Vector3< real_t >& point) const;
+   bool contains(const AABB& aabb) const;
+
+   Setup getSetup(){ return setup_; }
+   void setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+
+ private:
+   uint_t tpmsType;
+   Setup setup_;
+}; // class TPMS
+
+class TPMSRefinementSelection
+{
+ public:
+   TPMSRefinementSelection(const TPMS& tpms, const uint_t level): tpms_(tpms), level_(level){
+      auto setup = tpms_.getSetup();
+      const real_t start = setup.tpmsBeginn;
+      const real_t end   = setup.tpmsEnd;
+
+      boundingBox_ = AABB(start, real_c(0.0), real_c(0.0),
+                          end  , setup.domainSize[1], setup.domainSize[2]);
+   }
+
+   void operator()(SetupBlockForest& forest)
+   {
+      if(level_ == 0)
+         return;
+      for (auto block = forest.begin(); block != forest.end(); ++block){
+         const AABB& aabb = block->getAABB();
+
+         if (block->getLevel() < level_ && aabb.intersects(boundingBox_) )
+            block->setMarker(true);
+      }
+   }
+
+ private:
+   TPMS tpms_;
+   uint_t level_;
+   AABB boundingBox_;
+
+}; // class TPMSRefinementSelection
+
+
+} // namespace walberla
\ No newline at end of file
diff --git a/python/lbmpy_walberla/templates/LBMSweepCollection.tmpl.h b/python/lbmpy_walberla/templates/LBMSweepCollection.tmpl.h
index a85b9a93bbb676824b7e82ce3a075cbde1979925..c053e7120ebdacdcc3b8a2e31a51134f235969ac 100644
--- a/python/lbmpy_walberla/templates/LBMSweepCollection.tmpl.h
+++ b/python/lbmpy_walberla/templates/LBMSweepCollection.tmpl.h
@@ -86,8 +86,7 @@ class {{class_name}}
       {{kernel_list|generate_constructor(parameter_registration=parameter_scaling) |indent(6)}}
       validInnerOuterSplit_ = true;
 
-      for (auto& iBlock : *blocks)
-      {
+      for (auto& iBlock : *blocks){
          if (int_c(blocks->getNumberOfXCells(iBlock)) <= outerWidth_[0] * 2 ||
              int_c(blocks->getNumberOfYCells(iBlock)) <= outerWidth_[1] * 2 ||
              int_c(blocks->getNumberOfZCells(iBlock)) <= outerWidth_[2] * 2)
@@ -95,8 +94,7 @@ class {{class_name}}
       }
    }
 
-   void initialiseBlockPointer()
-   {
+   void initialiseBlockPointer(){
       {%if block_stream_collide -%}
       blockWise_ = true;
 
@@ -115,8 +113,7 @@ class {{class_name}}
       {%if target is equalto 'gpu' -%} glPropagationBlock_.resize(blocks_->getNumberOfLevels()); {% endif %}
       {%if target is equalto 'gpu' -%} glPropagationGrid_.resize(blocks_->getNumberOfLevels()); {% endif %}
 
-      for( auto it = blocks_->begin(); it != blocks_->end(); ++it )
-      {
+      for( auto it = blocks_->begin(); it != blocks_->end(); ++it ){
          auto* local = dynamic_cast< Block* >(it.get());
          {%- for field in block_stream_collide['all_fields'] %}
          auto {{field.name}} = local->getData< {{field | field_type(is_gpu=is_gpu)}} >({{field.name}}ID);
@@ -137,8 +134,7 @@ class {{class_name}}
          break;
       }
 
-      for( auto it = blocks_->begin(); it != blocks_->end(); ++it )
-      {
+      for( auto it = blocks_->begin(); it != blocks_->end(); ++it ){
          auto* local = dynamic_cast< Block* >(it.get());
          const uint_t level = local->getLevel();
          {%- for field in block_stream_collide['all_fields'] %}
@@ -209,17 +205,18 @@ class {{class_name}}
    {%if block_stream_collide -%}
    ~{{class_name}}() {
       {%if target is equalto 'gpu' -%}
-      for (uint_t level = 0; level < blocks_->getNumberOfLevels(); level++)
-      {
-         {%- for field in block_stream_collide['all_fields'] %}
-         if(!{{field.name}}Pointers[level].empty()){
-            WALBERLA_GPU_CHECK(gpuFree({{field.name}}PointersGPU[level]))
-         }
-         {%- endfor %}
+      if(blockWise_){
+         for (uint_t level = 0; level < blocks_->getNumberOfLevels(); level++){
+            {%- for field in block_stream_collide['all_fields'] %}
+            if(!{{field.name}}Pointers[level].empty()){
+               WALBERLA_GPU_CHECK(gpuFree({{field.name}}PointersGPU[level]))
+            }
+            {%- endfor %}
 
-         for (auto it = glPropagationPDFs[level].begin(); it != glPropagationPDFs[level].end(); it++){
-            if(it->second.empty()){ continue;}
-            WALBERLA_GPU_CHECK(gpuFree(glPropagationPDFsGPU[level][it->first]))
+            for (auto it = glPropagationPDFs[level].begin(); it != glPropagationPDFs[level].end(); it++){
+               if(it->second.empty()){ continue;}
+               WALBERLA_GPU_CHECK(gpuFree(glPropagationPDFsGPU[level][it->first]))
+            }
          }
       }
       {%- endif %}