From aefc230c12c1568e876bd9344c1f88e8f59ad87a Mon Sep 17 00:00:00 2001
From: Soeren Peters <peters@irmb.tu-bs.de>
Date: Wed, 30 May 2018 14:00:28 +0200
Subject: [PATCH] - qs on multiple level works

---
 .../grid/BoundaryConditions/Side.cpp          | 118 ++++++++++--------
 .../grid/BoundaryConditions/Side.h            |  19 +--
 .../grid/GridBuilder/LevelGridBuilder.cpp     |  62 ++++-----
 .../grid/GridBuilder/LevelGridBuilder.h       |   2 +-
 .../grid/GridBuilder/MultipleGridBuilder.cpp  |  14 ++-
 .../Calculation/UpdateGrid27.cpp              |  14 +--
 targets/apps/HULC/main.cpp                    |  16 +--
 7 files changed, 132 insertions(+), 113 deletions(-)

diff --git a/src/GridGenerator/grid/BoundaryConditions/Side.cpp b/src/GridGenerator/grid/BoundaryConditions/Side.cpp
index 479b20388..3bcca926c 100644
--- a/src/GridGenerator/grid/BoundaryConditions/Side.cpp
+++ b/src/GridGenerator/grid/BoundaryConditions/Side.cpp
@@ -59,30 +59,28 @@ uint Side::getIndex(SPtr<Grid> grid, std::string coord, real constant, real v1,
 }
 
 
-void Geometry::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void Geometry::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     auto geometryBoundaryCondition = std::dynamic_pointer_cast<GeometryBoundaryCondition>(boundaryCondition);
 
     std::vector<real> qNode(27);
     bool qFound = false;
 
-    for (int i = 0; i < grid->getSize(); i++)
+    for (int i = 0; i < grid[level]->getSize(); i++)
     {
-        if (grid->getFieldEntry(i) != BC_GEOMETRY)
+        if (grid[level]->getFieldEntry(i) != BC_GEOMETRY)
             continue;
 
-        for (int dir = 0; dir < grid->getEndDirection(); dir++)
+        for (int dir = 0; dir < grid[level]->getEndDirection(); dir++)
         {
-            const int qIndex = dir * grid->getSize() + i;
-            const real q = grid->getDistribution()[qIndex];
+            const int qIndex = dir * grid[level]->getSize() + i;
+            const real q = grid[level]->getDistribution()[qIndex];
 
             qNode[dir] = q;
             if (q > 0)
                 qFound = true;
             else
                 qNode[dir] = -1.0;
-
-
         }
 
         if (qFound)
@@ -97,87 +95,101 @@ void Geometry::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondi
 
 
 
-void MX::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void MX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
-    real startInner = grid->getStartY();
-    real endInner = grid->getEndY();
+    real startInner = grid[level]->getStartY();
+    real endInner = grid[level]->getEndY();
 
-    real startOuter = grid->getStartZ();
-    real endOuter = grid->getEndZ();
+    real startOuter = grid[level]->getStartZ();
+    real endOuter = grid[level]->getEndZ();
 
-    real coords[3] = { grid->getStartX(), grid->getStartY() + (grid->getEndY() - grid->getStartY()) / 2.0, grid->getStartZ() + (grid->getEndZ() - grid->getStartZ()) / 2.0 };
-    real startCoord = grid->getFirstFluidNode(coords, 0, grid->getStartX());
+    real coords[3] = { grid[level]->getStartX(), grid[level]->getStartY() + (grid[level]->getEndY() - grid[level]->getStartY()) / 2.0, grid[level]->getStartZ() + (grid[level]->getEndZ() - grid[level]->getStartZ()) / 2.0 };
+    real startCoord = grid[level]->getFirstFluidNode(coords, level, grid[level]->getStartX());
 
-    Side::addIndices(grid, boundaryCondition, "x", startCoord, startInner,
+    if(!isBoundaryOnFineGrid(level, grid, startCoord))
+        return;
+
+    Side::addIndices(grid[level], boundaryCondition, "x", startCoord, startInner,
         endInner, startOuter, endOuter);
 }
 
 
 
-void PX::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+bool MX::isBoundaryOnFineGrid(uint level, std::vector<SPtr<Grid>> grid, real startCoord)
 {
-    real startInner = grid->getStartY();
-    real endInner = grid->getEndY();
+    if (level > 0) {
+        real coords[3] = { grid[level - 1]->getStartX(), grid[level - 1]->getStartY() + (grid[level - 1]->getEndY() - grid[level - 1]->getStartY()) / 2.0, grid[level - 1]->getStartZ() + (grid[level - 1]->getEndZ() - grid[level - 1]->getStartZ()) / 2.0 };
+        real startCoordCoarser = grid[level - 1]->getFirstFluidNode(coords, level, grid[level - 1]->getStartX());
+        if (startCoord > startCoordCoarser)
+            return false;
+    }
+}
 
-    real startOuter = grid->getStartZ();
-    real endOuter = grid->getEndZ();
 
-    real coords[3] = { grid->getEndX(), grid->getStartY() + (grid->getEndY() - grid->getStartY()) / 2.0, grid->getStartZ() + (grid->getEndZ() - grid->getStartZ()) / 2.0 };
-    real startCoord = grid->getLastFluidNode(coords, 0, grid->getEndX());
+void PX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+{
+    //real startInner = grid->getStartY();
+    //real endInner = grid->getEndY();
 
-    Side::addIndices(grid, boundaryCondition, "x", startCoord, startInner,
-        endInner, startOuter, endOuter);
+    //real startOuter = grid->getStartZ();
+    //real endOuter = grid->getEndZ();
+
+    //real coords[3] = { grid->getEndX(), grid->getStartY() + (grid->getEndY() - grid->getStartY()) / 2.0, grid->getStartZ() + (grid->getEndZ() - grid->getStartZ()) / 2.0 };
+    //real startCoord = grid->getLastFluidNode(coords, 0, grid->getEndX());
+
+    //Side::addIndices(grid, boundaryCondition, "x", startCoord, startInner,
+    //    endInner, startOuter, endOuter);
 }
 
-void MY::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void MY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
-    real startInner = grid->getStartX();
-    real endInner = grid->getEndX();
+    //real startInner = grid->getStartX();
+    //real endInner = grid->getEndX();
 
-    real startOuter = grid->getStartZ();
-    real endOuter = grid->getEndZ();
+    //real startOuter = grid->getStartZ();
+    //real endOuter = grid->getEndZ();
 
 
-    Side::addIndices(grid, boundaryCondition, "y", grid->getStartY() + grid->getDelta(), startInner,
-        endInner, startOuter, endOuter);
+    //Side::addIndices(grid, boundaryCondition, "y", grid->getStartY() + grid->getDelta(), startInner,
+    //    endInner, startOuter, endOuter);
 }
 
 
-void PY::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void PY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
-    real startInner = grid->getStartX();
-    real endInner = grid->getEndX();
+    //real startInner = grid->getStartX();
+    //real endInner = grid->getEndX();
 
-    real startOuter = grid->getStartZ();
-    real endOuter = grid->getEndZ();
+    //real startOuter = grid->getStartZ();
+    //real endOuter = grid->getEndZ();
 
 
-    Side::addIndices(grid, boundaryCondition, "y", grid->getEndY() - grid->getDelta(), startInner,
-        endInner, startOuter, endOuter);
+    //Side::addIndices(grid, boundaryCondition, "y", grid->getEndY() - grid->getDelta(), startInner,
+    //    endInner, startOuter, endOuter);
 }
 
 
-void MZ::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void MZ::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
-    real startInner = grid->getStartX();
-    real endInner = grid->getEndX();
+    //real startInner = grid->getStartX();
+    //real endInner = grid->getEndX();
 
-    real startOuter = grid->getStartY();
-    real endOuter = grid->getEndY();
+    //real startOuter = grid->getStartY();
+    //real endOuter = grid->getEndY();
 
-    Side::addIndices(grid, boundaryCondition, "z", grid->getStartZ() + grid->getDelta(), startInner,
-        endInner, startOuter, endOuter);
+    //Side::addIndices(grid, boundaryCondition, "z", grid->getStartZ() + grid->getDelta(), startInner,
+    //    endInner, startOuter, endOuter);
 }
 
-void PZ::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition)
+void PZ::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
-    real startInner = grid->getStartX();
-    real endInner = grid->getEndX();
+    //real startInner = grid->getStartX();
+    //real endInner = grid->getEndX();
 
-    real startOuter = grid->getStartY();
-    real endOuter = grid->getEndY();
+    //real startOuter = grid->getStartY();
+    //real endOuter = grid->getEndY();
 
 
-    Side::addIndices(grid, boundaryCondition, "z", grid->getEndZ() - grid->getDelta(), startInner,
-        endInner, startOuter, endOuter);
+    //Side::addIndices(grid, boundaryCondition, "z", grid->getEndZ() - grid->getDelta(), startInner,
+    //    endInner, startOuter, endOuter);
 }
diff --git a/src/GridGenerator/grid/BoundaryConditions/Side.h b/src/GridGenerator/grid/BoundaryConditions/Side.h
index d4a44e798..175348923 100644
--- a/src/GridGenerator/grid/BoundaryConditions/Side.h
+++ b/src/GridGenerator/grid/BoundaryConditions/Side.h
@@ -2,6 +2,7 @@
 #define SIDE_H
 
 #include <string>
+#include <vector>
 
 #include <VirtualFluidsDefinitions.h>
 
@@ -30,7 +31,7 @@ enum class SideType
 class Side
 {
 public:
-    virtual void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) = 0;
+    virtual void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) = 0;
 
     virtual int getCoordinate() const = 0;
     virtual int getDirection() const = 0;
@@ -48,7 +49,7 @@ private:
 class Geometry : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -64,7 +65,9 @@ public:
 class MX : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
+
+    bool isBoundaryOnFineGrid(uint level, std::vector<SPtr<Grid>> grid, real startCoord);
 
     int getCoordinate() const override
     {
@@ -80,7 +83,7 @@ public:
 class PX : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -97,7 +100,7 @@ public:
 class MY : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -113,7 +116,7 @@ public:
 class PY : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -130,7 +133,7 @@ public:
 class MZ : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -146,7 +149,7 @@ public:
 class PZ : public Side
 {
 public:
-    void addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition) override;
+    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
index 65d3da23e..b82352e72 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
@@ -66,9 +66,9 @@ void LevelGridBuilder::setVelocityBoundaryCondition(SideType sideType, real vx,
             auto side = SideFactory::make(sideType);
 
             velocityBoundaryCondition->side = side;
-            velocityBoundaryCondition->side->addIndices(grids[level], velocityBoundaryCondition);
+            velocityBoundaryCondition->side->addIndices(grids, level, velocityBoundaryCondition);
 
-            boundaryConditions[level].velocityBoundaryConditions.push_back(velocityBoundaryCondition);
+            boundaryConditions[level]->velocityBoundaryConditions.push_back(velocityBoundaryCondition);
         }
     }
 
@@ -80,10 +80,10 @@ void LevelGridBuilder::setVelocityGeometryBoundaryCondition(real vx, real vy, re
 
     for (int level = 0; level < getNumberOfGridLevels(); level++)
     {
-        boundaryConditions[level].geometryBoundaryCondition->vx = vx;
-        boundaryConditions[level].geometryBoundaryCondition->vy = vy;
-        boundaryConditions[level].geometryBoundaryCondition->vz = vz;
-        boundaryConditions[level].geometryBoundaryCondition->side->addIndices(grids[level], boundaryConditions[level].geometryBoundaryCondition);
+        boundaryConditions[level]->geometryBoundaryCondition->vx = vx;
+        boundaryConditions[level]->geometryBoundaryCondition->vy = vy;
+        boundaryConditions[level]->geometryBoundaryCondition->vz = vz;
+        boundaryConditions[level]->geometryBoundaryCondition->side->addIndices(grids, level, boundaryConditions[level]->geometryBoundaryCondition);
     }
 }
 
@@ -95,9 +95,9 @@ void LevelGridBuilder::setPressureBoundaryCondition(SideType sideType, real rho)
 
         auto side = SideFactory::make(sideType);
         pressureBoundaryCondition->side = side;
-        pressureBoundaryCondition->side->addIndices(grids[level], pressureBoundaryCondition);
+        pressureBoundaryCondition->side->addIndices(grids, level, pressureBoundaryCondition);
 
-        boundaryConditions[level].pressureBoundaryConditions.push_back(pressureBoundaryCondition);
+        boundaryConditions[level]->pressureBoundaryConditions.push_back(pressureBoundaryCondition);
     }
 }
 
@@ -115,9 +115,9 @@ void LevelGridBuilder::setNoSlipBoundaryCondition(SideType sideType)
         auto side = SideFactory::make(sideType);
 
         noSlipBoundaryCondition->side = side;
-        noSlipBoundaryCondition->side->addIndices(grids[level], noSlipBoundaryCondition);
+        noSlipBoundaryCondition->side->addIndices(grids, level, noSlipBoundaryCondition);
 
-        boundaryConditions[level].noSlipBoundaryConditions.push_back(noSlipBoundaryCondition);
+        boundaryConditions[level]->noSlipBoundaryConditions.push_back(noSlipBoundaryCondition);
     }
 }
 
@@ -275,7 +275,7 @@ void LevelGridBuilder::setOutflowValues(real* RhoBC, int* kN, int channelSide, i
 uint LevelGridBuilder::getVelocitySize(int level) const
 {
     uint size = 0;
-    for (auto boundaryCondition : boundaryConditions[level].velocityBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->velocityBoundaryConditions)
     {
         size += uint(boundaryCondition->indices.size());
     }
@@ -285,7 +285,7 @@ uint LevelGridBuilder::getVelocitySize(int level) const
 void LevelGridBuilder::getVelocityValues(real* vx, real* vy, real* vz, int* indices, int level) const
 {
     int allIndicesCounter = 0;
-    for (auto boundaryCondition : boundaryConditions[level].velocityBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->velocityBoundaryConditions)
     {
         for(int i = 0; i < boundaryCondition->indices.size(); i++)
         {
@@ -302,7 +302,7 @@ void LevelGridBuilder::getVelocityValues(real* vx, real* vy, real* vz, int* indi
 uint LevelGridBuilder::getPressureSize(int level) const
 {
     uint size = 0;
-    for (auto boundaryCondition : boundaryConditions[level].pressureBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->pressureBoundaryConditions)
     {
         size += uint(boundaryCondition->indices.size());
     }
@@ -312,7 +312,7 @@ uint LevelGridBuilder::getPressureSize(int level) const
 void LevelGridBuilder::getPressureValues(real* rho, int* indices, int* neighborIndices, int level) const
 {
     int allIndicesCounter = 0;
-    for (auto boundaryCondition : boundaryConditions[level].pressureBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->pressureBoundaryConditions)
     {
         for (int i = 0; i < boundaryCondition->indices.size(); i++)
         {
@@ -330,7 +330,7 @@ void LevelGridBuilder::getPressureValues(real* rho, int* indices, int* neighborI
 void LevelGridBuilder::getVelocityQs(real* qs[27], int level) const
 {
     int allIndicesCounter = 0;
-    for (auto boundaryCondition : boundaryConditions[level].velocityBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->velocityBoundaryConditions)
     {
         for (int i = 0; i < boundaryCondition->indices.size(); i++)
         {
@@ -350,7 +350,7 @@ void LevelGridBuilder::getVelocityQs(real* qs[27], int level) const
 void LevelGridBuilder::getPressureQs(real* qs[27], int level) const
 {
     int allIndicesCounter = 0;
-    for (auto boundaryCondition : boundaryConditions[level].pressureBoundaryConditions)
+    for (auto boundaryCondition : boundaryConditions[level]->pressureBoundaryConditions)
     {
         for (int i = 0; i < boundaryCondition->indices.size(); i++)
         {
@@ -369,17 +369,17 @@ void LevelGridBuilder::getPressureQs(real* qs[27], int level) const
 
 uint LevelGridBuilder::getGeometrySize(int level) const
 {
-    if (boundaryConditions[level].geometryBoundaryCondition)
-        return  boundaryConditions[level].geometryBoundaryCondition->indices.size();
+    if (boundaryConditions[level]->geometryBoundaryCondition)
+        return  boundaryConditions[level]->geometryBoundaryCondition->indices.size();
     
     return 0;
 }
 
 void LevelGridBuilder::getGeometryIndices(int* indices, int level) const
 {
-    for (uint i = 0; i <  boundaryConditions[level].geometryBoundaryCondition->indices.size(); i++)
+    for (uint i = 0; i <  boundaryConditions[level]->geometryBoundaryCondition->indices.size(); i++)
     {
-        indices[i] = grids[level]->getSparseIndex(boundaryConditions[level].geometryBoundaryCondition->indices[i]) + 1;
+        indices[i] = grids[level]->getSparseIndex(boundaryConditions[level]->geometryBoundaryCondition->indices[i]) + 1;
     }
 }
 
@@ -391,22 +391,22 @@ bool LevelGridBuilder::hasGeometryValues() const
 
 void LevelGridBuilder::getGeometryValues(real* vx, real* vy, real* vz, int level) const
 {
-    for (uint i = 0; i < boundaryConditions[level].geometryBoundaryCondition->indices.size(); i++)
+    for (uint i = 0; i < boundaryConditions[level]->geometryBoundaryCondition->indices.size(); i++)
     {
-        vx[i] = boundaryConditions[level].geometryBoundaryCondition->vx;
-        vy[i] = boundaryConditions[level].geometryBoundaryCondition->vy;
-        vz[i] = boundaryConditions[level].geometryBoundaryCondition->vz;
+        vx[i] = boundaryConditions[level]->geometryBoundaryCondition->vx;
+        vy[i] = boundaryConditions[level]->geometryBoundaryCondition->vy;
+        vz[i] = boundaryConditions[level]->geometryBoundaryCondition->vz;
     }
 }
 
 
 void LevelGridBuilder::getGeometryQs(real* qs[27], int level) const
 {
-    for (int i = 0; i < boundaryConditions[level].geometryBoundaryCondition->indices.size(); i++)
+    for (int i = 0; i < boundaryConditions[level]->geometryBoundaryCondition->indices.size(); i++)
     {
         for (int dir = 0; dir <= grids[level]->getEndDirection(); dir++)
         {
-            qs[dir][i] = boundaryConditions[level].geometryBoundaryCondition->qs[i][dir];
+            qs[dir][i] = boundaryConditions[level]->geometryBoundaryCondition->qs[i][dir];
         }
     }
 }
@@ -626,15 +626,15 @@ void LevelGridBuilder::writeArrows(std::string fileName) const
     std::vector<UbTupleInt2> cells;
 
     int actualNodeNumber = 0;
-    for (int index = 0; index < boundaryConditions[getNumberOfGridLevels() - 1].geometryBoundaryCondition->indices.size(); index++)
+    for (int index = 0; index < boundaryConditions[getNumberOfGridLevels() - 1]->geometryBoundaryCondition->indices.size(); index++)
     {
-        Vertex startNode = getVertex(boundaryConditions[getNumberOfGridLevels() - 1].geometryBoundaryCondition->indices[index]);
+        Vertex startNode = getVertex(boundaryConditions[getNumberOfGridLevels() - 1]->geometryBoundaryCondition->indices[index]);
         for (int qi = 0; qi <= 26; qi++)
         {
-            real qval = boundaryConditions[getNumberOfGridLevels() - 1].geometryBoundaryCondition->qs[index][qi];
+            real qval = boundaryConditions[getNumberOfGridLevels() - 1]->geometryBoundaryCondition->qs[index][qi];
             if (qval > 0.0f)
             {
-                Vertex dir((real)grids[0]->getDirection()[qi * DIMENSION + 0], (real)grids[0]->getDirection()[qi* DIMENSION + 1], (real)grids[0]->getDirection()[qi * DIMENSION + 2]);
+                Vertex dir((real)grids[0]->getDirection()[qi * DIMENSION + 0], (real)grids[0]->getDirection()[qi * DIMENSION + 1], (real)grids[0]->getDirection()[qi * DIMENSION + 2]);
                 Vertex nodeOnGeometry(startNode + (dir * qval));
 
                 nodes.push_back(makeUbTuple(float(startNode.x), float(startNode.y), float(startNode.z)));
@@ -652,7 +652,7 @@ void LevelGridBuilder::writeArrows(std::string fileName) const
 Vertex LevelGridBuilder::getVertex(int matrixIndex) const
 {
     real x, y, z;
-    this->grids[0]->transIndexToCoords(matrixIndex, x, y, z);
+    this->grids[getNumberOfGridLevels() - 1]->transIndexToCoords(matrixIndex, x, y, z);
     return Vertex(x,y,z);
 }
 
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
index 80b4bd5b3..8a289ece6 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
@@ -108,7 +108,7 @@ protected:
     bool geometryHasValues = false;
 
     std::vector<std::shared_ptr<Grid> > grids;
-    std::vector<BoundaryConditions> boundaryConditions;
+    std::vector<SPtr<BoundaryConditions> > boundaryConditions;
 
 
     void checkLevel(int level);
diff --git a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
index 83f62061a..b2df06da1 100644
--- a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
+++ b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
@@ -26,7 +26,7 @@ SPtr<MultipleGridBuilder> MultipleGridBuilder::makeShared(SPtr<GridFactory> grid
 
 void MultipleGridBuilder::addCoarseGrid(real startX, real startY, real startZ, real endX, real endY, real endZ, real delta)
 {
-    boundaryConditions.push_back(BoundaryConditions());
+    boundaryConditions.push_back(SPtr<BoundaryConditions>(new BoundaryConditions));
 
     const auto grid = this->makeGrid(new Cuboid(startX, startY, startZ, endX, endY, endZ), startX, startY, startZ, endX, endY, endZ, delta);
     addGridToList(grid);
@@ -38,8 +38,8 @@ void MultipleGridBuilder::addGeometry(Object* solidObject)
 
     for(auto bcs : boundaryConditions)
     {
-        bcs.geometryBoundaryCondition = GeometryBoundaryCondition::make();
-        bcs.geometryBoundaryCondition->side = SideFactory::make(SideType::GEOMETRY);
+        bcs->geometryBoundaryCondition = GeometryBoundaryCondition::make();
+        bcs->geometryBoundaryCondition->side = SideFactory::make(SideType::GEOMETRY);
     }
 }
 
@@ -57,7 +57,7 @@ void MultipleGridBuilder::addGeometry(Object* solidObject, uint level)
 
 void MultipleGridBuilder::addGrid(Object* gridShape)
 {
-    boundaryConditions.push_back(BoundaryConditions());
+    boundaryConditions.push_back(SPtr<BoundaryConditions>(new BoundaryConditions));
 
     if (!coarseGridExists())
         return emitNoCoarseGridExistsWarning();
@@ -69,7 +69,7 @@ void MultipleGridBuilder::addGrid(Object* gridShape)
 
 void MultipleGridBuilder::addGrid(Object* gridShape, uint levelFine)
 {
-    boundaryConditions.push_back(BoundaryConditions());
+    boundaryConditions.push_back(SPtr<BoundaryConditions>(new BoundaryConditions));
 
     if (!coarseGridExists())
         return emitNoCoarseGridExistsWarning();
@@ -272,14 +272,16 @@ void MultipleGridBuilder::buildGrids()
         grids[grids.size() - 1]->findQs(solidObject);
     }
 
+
     for (size_t i = 0; i < grids.size() - 1; i++)
         grids[i]->findGridInterface(grids[i + 1]);
 
+
     for (size_t i = 0; i < grids.size() - 1; i++)
         grids[i]->findSparseIndices(grids[i + 1]);
 
-    grids[grids.size() - 1]->findSparseIndices(nullptr);
 
+    grids[grids.size() - 1]->findSparseIndices(nullptr);
 
 
     //for(auto velocityBC : velocityBoundaryConditions)
diff --git a/src/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index e49b71a2f..3f1b24d19 100644
--- a/src/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -742,13 +742,13 @@ void updateGrid27(Parameter* para, Communicator* comm, PorousMedia** pm, int lev
 		  //          para->getParD(level)->size_Mat_SP,     para->getParD(level)->evenOrOdd);
 		  //getLastCudaError("QVelDev27 execution failed");
 
-		     // QVelDevComp27(para->getParD(level)->numberofthreads, para->getParD(level)->nx,           para->getParD(level)->ny,
-							//para->getParD(level)->QGeom.Vx,        para->getParD(level)->QGeom.Vy,     para->getParD(level)->QGeom.Vz,
-							//para->getParD(level)->d0SP.f[0],       para->getParD(level)->QGeom.k,      para->getParD(level)->QGeom.q27[0], 
-							//para->getParD(level)->QGeom.kQ,        para->getParD(level)->QGeom.kQ,     para->getParD(level)->omega,
-							//para->getParD(level)->neighborX_SP,    para->getParD(level)->neighborY_SP, para->getParD(level)->neighborZ_SP,
-							//para->getParD(level)->size_Mat_SP,     para->getParD(level)->evenOrOdd);
-		     // getLastCudaError("QVelDevComp27 execution failed");
+		      QVelDevComp27(para->getParD(level)->numberofthreads, para->getParD(level)->nx,           para->getParD(level)->ny,
+							para->getParD(level)->QGeom.Vx,        para->getParD(level)->QGeom.Vy,     para->getParD(level)->QGeom.Vz,
+							para->getParD(level)->d0SP.f[0],       para->getParD(level)->QGeom.k,      para->getParD(level)->QGeom.q27[0], 
+							para->getParD(level)->QGeom.kQ,        para->getParD(level)->QGeom.kQ,     para->getParD(level)->omega,
+							para->getParD(level)->neighborX_SP,    para->getParD(level)->neighborY_SP, para->getParD(level)->neighborZ_SP,
+							para->getParD(level)->size_Mat_SP,     para->getParD(level)->evenOrOdd);
+		      getLastCudaError("QVelDevComp27 execution failed");
 		 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
          if (para->getParD(level)->evenOrOdd==true)  para->getParD(level)->evenOrOdd=false;
          else                                        para->getParD(level)->evenOrOdd=true;
diff --git a/targets/apps/HULC/main.cpp b/targets/apps/HULC/main.cpp
index 6055d6690..e7833c4b9 100644
--- a/targets/apps/HULC/main.cpp
+++ b/targets/apps/HULC/main.cpp
@@ -276,21 +276,25 @@ void multipleLevel(const std::string& configPath)
     //gridBuilder->addCoarseGrid(-16, -14, -14, 59, 28, 29, 1.0);
     //TriangularMesh* triangularMesh = TriangularMesh::make("D:/GRIDGENERATION/STL/input/local_input/bruecke.stl");
 
-    gridBuilder->addCoarseGrid(-10, -10, -10, 9, 9, 9, 0.25);
+    gridBuilder->addCoarseGrid(-14, -14, -14, 14, 14, 14, 1.0);
     TriangularMesh* triangularMesh = TriangularMesh::make("D:/GRIDGENERATION/STL/cubeBinaer1x1.stl");
+    gridBuilder->addGrid(new Cuboid(-4, -4, -4, 5, 5, 5), 1);
 
 
     gridBuilder->addGeometry(triangularMesh);
 
-    gridBuilder->addGrid(new Cuboid(-2,-2,-2,2,2,2));
+    gridBuilder->setPeriodicBoundaryCondition(true, true, true);
 
 
     gridBuilder->buildGrids(); // buildGrids() has to be called before setting the BCs!!!!
 
-    gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
-    gridBuilder->setVelocityBoundaryCondition(SideType::MX, 0.001, 0.0, 0.0);
+    gridBuilder->writeGridsToVtk("D:/GRIDGENERATION/");
+    //gridBuilder->writeArrows("D:/arrows");
+
+
+    //gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
+    //gridBuilder->setVelocityBoundaryCondition(SideType::MX, 0.001, 0.0, 0.0);
 
-    gridBuilder->setPeriodicBoundaryCondition(false, true, true);
 
     //gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
 
@@ -327,8 +331,6 @@ void multipleLevel(const std::string& configPath)
 
     //SimulationFileWriter::write("D:/GRIDGENERATION/files/", gridBuilder, FILEFORMAT::ASCII);
 
-    //gridBuilder->writeGridsToVtk("D:/GRIDGENERATION/");
-    //gridBuilder->writeArrows("D:/arrows");
 
 
     SPtr<Parameter> para = Parameter::make();
-- 
GitLab