Merge branch 'develop' of https://git.rz.tu-bs.de/irmb/VirtualFluids_dev into develop

apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp

@@ -164,8 +164,8 @@ int main()
         gridBuilder->setNoSlipBoundaryCondition(SideType::MX);
         gridBuilder->setNoSlipBoundaryCondition(SideType::PY);
         gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
-        gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vxLB, vyLB, 0.0);
         gridBuilder->setNoSlipBoundaryCondition(SideType::MZ);
+        gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vxLB, vyLB, 0.0);
 
         BoundaryConditionFactory bcFactory;
 

apps/gpu/LBM/DrivenCavityUniform/CMakeLists.txt

-PROJECT(DrivenCavity LANGUAGES CUDA CXX)
+PROJECT(DrivenCavityUniform LANGUAGES CUDA CXX)
 
 #LIST(APPEND CS_COMPILER_FLAGS_CXX "-DOMPI_SKIP_MPICXX" )
 
@@ -6,4 +6,5 @@ vf_add_library(BUILDTYPE binary PRIVATE_LINK basics VirtualFluids_GPU GridGenera
 
 set_source_files_properties(DrivenCavity.cpp PROPERTIES LANGUAGE CUDA)
 
-set_target_properties(DrivenCavity PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
+set_target_properties(DrivenCavityUniform PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
+

regression-tests/driven_cavity_test.sh

@@ -7,7 +7,7 @@
 # build VirtualFluids accordingly to our specific test scenario.
 # in this case adding -DUSER_APPS="apps/gpu/LBM/DrivenCavity to the cmake command is not necessary, because the DrivenCavity is added to VirtualFluids by default.
 mkdir -p build
-cmake -B build --preset=release_make_gpu -DCMAKE_CUDA_ARCHITECTURES=75 #-DUSER_APPS="apps/gpu/LBM/DrivenCavity"
+cmake -B build --preset=make_gpu -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=75 #-DUSER_APPS="apps/gpu/LBM/DrivenCavity"
 cmake --build build --parallel 8
 
 # execute VirtualFluids

regression-tests/driven_cavity_uniform_test.sh

@@ -7,18 +7,19 @@
 # build VirtualFluids accordingly to our specific test scenario.
 # in this case adding -DUSER_APPS="apps/gpu/LBM/DrivenCavity to the cmake command is not necessary, because the DrivenCavity is added to VirtualFluids by default.
 mkdir -p build
-cmake -B build --preset=release_make_gpu -DCMAKE_CUDA_ARCHITECTURES=75 -DUSER_APPS="apps/gpu/LBM/DrivenCavityUniform"
+cmake -B build --preset=make_gpu -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=75 -DUSER_APPS="apps/gpu/LBM/DrivenCavityUniform"
 cmake --build build --parallel 8
 
 # execute VirtualFluids
-./build/bin/DrivenCavity
+./build/bin/DrivenCavityUniform
+
 
 # set the path to the produced data
 PATH_TO_DIR=output/DrivenCavity_uniform
 
 # set the path to the reference data.
 # `regression-tests/reference_data` is fix `regression_tests/gpu/DrivenCavity_uniform_2022_12_16` must match the structure in https://github.com/irmb/test_data:
-PATH_TO_REFERENCE_DIR=regression-tests/reference_data/regression_tests/gpu/DrivenCavity_uniform_2022_12_16
+PATH_TO_REFERENCE_DIR=regression-tests/reference_data/regression_tests/gpu/DrivenCavity_uniform
 
 # execute fieldcompare (A more comprehensive manual can be found here https://gitlab.com/dglaeser/fieldcompare)
 fieldcompare dir $PATH_TO_DIR $PATH_TO_REFERENCE_DIR --include-files "*.vtu"

src/basics/basics/writer/WbWriter.h

@@ -88,7 +88,12 @@ public:
     {
         throw UbException(UB_EXARGS, "not implemented for " + (std::string) typeid(*this).name());
     }
-
+    virtual std::string writeLinesWithLineData(const std::string & /*filename*/, std::vector<UbTupleFloat3> & /*nodes*/,
+                                               std::vector<UbTupleInt2> & /*lines*/, std::vector<std::string> & /*datanames*/,
+                                               std::vector<std::vector<float>> & /*celldata*/)
+    {
+        throw UbException(UB_EXARGS, "not implemented for " + (std::string) typeid(*this).name());
+    }
     //////////////////////////////////////////////////////////////////////////
     // triangles
     // cell numbering:

src/basics/basics/writer/WbWriterVtkXmlBinary.h

@@ -93,6 +93,9 @@ public:
     // nodedata);
     // FIXME: hides function in base class
 
+    std::string writeLinesWithLineData(const std::string &filename, std::vector<UbTupleFloat3> &nodes, std::vector<UbTupleInt2> &lines,
+                                       std::vector<std::string> &datanames, std::vector<std::vector<float>> &celldata) override;
+
     //////////////////////////////////////////////////////////////////////////
     // triangles
     //                    2

src/gpu/GridGenerator/grid/BoundaryConditions/Side.cpp

@@ -37,12 +37,15 @@
 #include "grid/NodeValues.h"
 
 #include "utilities/math/Math.h"
+#include <array>
+#include <cstddef>
+#include <vector>
 
 using namespace gg;
 
-std::vector<real> Side::getNormal()
+std::array<real, 3> Side::getNormal() const
 {
-    std::vector<real> normal;
+    std::array<real, 3> normal;
     if(this->getCoordinate()==X_INDEX)
         normal = {(real)this->getDirection(), 0.0, 0.0};
     if(this->getCoordinate()==Y_INDEX)
@@ -61,19 +64,17 @@ void Side::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition
         {
             const uint index = getIndex(grid, coord, constant, v1, v2);
 
-            if ((index != INVALID_INDEX) && (   grid->getFieldEntry(index) == vf::gpu::FLUID
+            if(index == INVALID_INDEX)
+                continue;
+
+            if (   grid->getFieldEntry(index) == vf::gpu::FLUID
                                             ||  grid->getFieldEntry(index) == vf::gpu::FLUID_CFC
                                             ||  grid->getFieldEntry(index) == vf::gpu::FLUID_CFF
                                             ||  grid->getFieldEntry(index) == vf::gpu::FLUID_FCC
                                             ||  grid->getFieldEntry(index) == vf::gpu::FLUID_FCF
                                             ||  grid->getFieldEntry(index) == vf::gpu::FLUID_FCF
-
-                                            //! Enforce overlap of BCs on edge nodes
-                                            ||  grid->getFieldEntry(index)  == vf::gpu::BC_PRESSURE
-                                            ||  grid->getFieldEntry(index)  == vf::gpu::BC_VELOCITY
-                                            ||  grid->getFieldEntry(index)  == vf::gpu::BC_NOSLIP
-                                            ||  grid->getFieldEntry(index)  == vf::gpu::BC_SLIP
-                                            ||  grid->getFieldEntry(index)  == vf::gpu::BC_STRESS ))
+                                            // Overlap of BCs on edge nodes
+                                            || grid->nodeHasBC(index) )
             {
                 grid->setFieldEntry(index, boundaryCondition->getType());
                 boundaryCondition->indices.push_back(index);
@@ -86,6 +87,10 @@ void Side::addIndices(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition
             }
         }
     }
+
+    const auto currentBCSide = this->whoAmI();
+    if(currentBCSide != SideType::GEOMETRY)
+        grid->addBCalreadySet(currentBCSide);
 }
 
 void Side::setPressureNeighborIndices(SPtr<BoundaryCondition> boundaryCondition, SPtr<Grid> grid, const uint index)
@@ -138,55 +143,111 @@ void Side::setStressSamplingIndices(SPtr<BoundaryCondition> boundaryCondition, S
 
 void Side::setQs(SPtr<Grid> grid, SPtr<BoundaryCondition> boundaryCondition, uint index,  real q)
 {
-
     std::vector<real> qNode(grid->getEndDirection() + 1);
 
-    for (int dir = 0; dir <= grid->getEndDirection(); dir++)
-    {
-        real x,y,z;
-        grid->transIndexToCoords( index, x, y, z );
+    for (int dir = 0; dir <= grid->getEndDirection(); dir++) {
+        real x, y, z;
+        grid->transIndexToCoords(index, x, y, z);
 
-        real coords[3] = {x,y,z};
+        std::array<real, 3> coords = { x, y, z };
+        std::array<real, 3> neighborCoords = getNeighborCoordinates(grid.get(), coords, (size_t)dir);
 
-        real neighborX = x + grid->getDirection()[dir * DIMENSION + 0] * grid->getDelta();
-        real neighborY = y + grid->getDirection()[dir * DIMENSION + 1] * grid->getDelta();
-        real neighborZ = z + grid->getDirection()[dir * DIMENSION + 2] * grid->getDelta();
+        correctNeighborForPeriodicBoundaries(grid.get(), coords, neighborCoords);
 
-        // correct neighbor coordinates in case of periodic boundaries
-        if( grid->getPeriodicityX() && grid->getFieldEntry( grid->transCoordToIndex( neighborX, y, z ) ) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY )
-        {
-            if( neighborX > x ) neighborX = grid->getFirstFluidNode( coords, 0, grid->getStartX() );
-            else                neighborX = grid->getLastFluidNode ( coords, 0, grid->getEndX() );
-        }
+        const uint neighborIndex = grid->transCoordToIndex(neighborCoords[0], neighborCoords[1], neighborCoords[2]);
 
-        if( grid->getPeriodicityY() && grid->getFieldEntry( grid->transCoordToIndex( x, neighborY, z ) ) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY )
-        {
-            if( neighborY > y ) neighborY = grid->getFirstFluidNode( coords, 1, grid->getStartY() );
-            else                neighborY = grid->getLastFluidNode ( coords, 1, grid->getEndY() );
+        //! Only setting q's that partially point in the Side-normal direction
+        const bool alignedWithNormal = this->isAlignedWithMyNormal(grid.get(), dir);
+        if (grid->isStopperForBC(neighborIndex) && alignedWithNormal) {
+            qNode[dir] = 0.5;
+        } else {
+            qNode[dir] = -1.0;
         }
 
-        if( grid->getPeriodicityZ() && grid->getFieldEntry( grid->transCoordToIndex( x, y, neighborZ ) ) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY )
-        {
-            if( neighborZ > z ) neighborZ = grid->getFirstFluidNode( coords, 2, grid->getStartZ() );
-            else                neighborZ = grid->getLastFluidNode ( coords, 2, grid->getEndZ() );
+        // reset diagonals in case they were set by another bc
+        resetDiagonalsInCaseOfOtherBC(grid.get(), qNode, dir, coords);
+    }
+
+    boundaryCondition->qs.push_back(qNode);
+}
+
+std::array<real, 3> Side::getNeighborCoordinates(Grid *grid, const std::array<real, 3> &coordinates, size_t direction) const
+{
+    return { coordinates[0] + grid->getDirection()[direction * DIMENSION + 0] * grid->getDelta(),
+             coordinates[1] + grid->getDirection()[direction * DIMENSION + 1] * grid->getDelta(),
+             coordinates[2] + grid->getDirection()[direction * DIMENSION + 2] * grid->getDelta() };
+}
+
+bool Side::neighborNormalToSideIsAStopper(Grid *grid, const std::array<real, 3> &coordinates, SideType side) const
+{
+    const auto neighborCoords = getNeighborCoordinates(grid, coordinates, sideToD3Q27.at(side));
+    const auto neighborIndex = grid->transCoordToIndex(neighborCoords[0], neighborCoords[1], neighborCoords[2]);
+    return grid->isStopperForBC(neighborIndex);
+}
+
+void Side::resetDiagonalsInCaseOfOtherBC(Grid *grid, std::vector<real> &qNode, int dir,
+                                         const std::array<real, 3> &coordinates) const
+{
+    // When to reset a diagonal q to -1:
+    // - it is normal to another boundary condition which was already set
+    // - and it actually is influenced by the other bc:
+    //   We check if its neighbor in the regular direction to the other bc is a stopper. If it is a stopper, it is influenced by the other bc.
+
+    if (qNode[dir] == 0.5 && grid->getBCAlreadySet().size() > 0) {
+        for (int i = 0; i < (int)grid->getBCAlreadySet().size(); i++) {
+            SideType otherDir = grid->getBCAlreadySet()[i];
+
+            // only reset normals for nodes on edges and corners, not on faces
+            if (!neighborNormalToSideIsAStopper(grid, coordinates, otherDir))
+                continue;
+
+            const auto otherNormal = normals.at(otherDir);
+            if (isAlignedWithNormal(grid, dir, otherNormal)) {
+                qNode[dir] = -1.0;
+            }
         }
+    }
+}
 
-        //! Only seting q's that partially point in the Side-normal direction
-        bool alignedWithNormal = (this->getNormal()[0]*grid->getDirection()[dir * DIMENSION + 0]+
-                                  this->getNormal()[1]*grid->getDirection()[dir * DIMENSION + 1]+
-                                  this->getNormal()[2]*grid->getDirection()[dir * DIMENSION + 2] ) > 0;
+bool Side::isAlignedWithMyNormal(const Grid *grid, int dir) const
+{
+    std::array<real, 3> normal = this->getNormal();
+    return isAlignedWithNormal(grid, dir, normal);
+}
 
-        uint neighborIndex = grid->transCoordToIndex( neighborX, neighborY, neighborZ );
-        if((grid->getFieldEntry(neighborIndex) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY ||
-            grid->getFieldEntry(neighborIndex) == vf::gpu::STOPPER_OUT_OF_GRID          ||
-            grid->getFieldEntry(neighborIndex) == vf::gpu::STOPPER_SOLID)               &&
-            alignedWithNormal )
-            qNode[dir] = q;
+bool Side::isAlignedWithNormal(const Grid *grid, int dir, const std::array<real, 3> &normal) const
+{
+    return (normal[0] * grid->getDirection()[dir * DIMENSION + 0] +
+            normal[1] * grid->getDirection()[dir * DIMENSION + 1] +
+            normal[2] * grid->getDirection()[dir * DIMENSION + 2]) > 0;
+}
+
+void Side::correctNeighborForPeriodicBoundaries(const Grid *grid, std::array<real, 3>& coords, std::array<real, 3>& neighborCoords) const
+{
+    // correct neighbor coordinates in case of periodic boundaries
+    if (grid->getPeriodicityX() &&
+        grid->getFieldEntry(grid->transCoordToIndex(neighborCoords[0], coords[1], coords[2])) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY) {
+        if (neighborCoords[0] > coords[0])
+            neighborCoords[0] = grid->getFirstFluidNode(coords.data(), 0, grid->getStartX());
         else
-            qNode[dir] = -1.0;
+            neighborCoords[0] = grid->getLastFluidNode(coords.data(), 0, grid->getEndX());
     }
 
-    boundaryCondition->qs.push_back(qNode);
+    if (grid->getPeriodicityY() &&
+        grid->getFieldEntry(grid->transCoordToIndex(coords[0], neighborCoords[1], coords[2])) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY) {
+        if (neighborCoords[1] > coords[1])
+            neighborCoords[1] = grid->getFirstFluidNode(coords.data(), 1, grid->getStartY());
+        else
+            neighborCoords[1] = grid->getLastFluidNode(coords.data(), 1, grid->getEndY());
+    }
+
+    if (grid->getPeriodicityZ() &&
+        grid->getFieldEntry(grid->transCoordToIndex(coords[0], coords[1], neighborCoords[2])) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY) {
+        if (neighborCoords[2] > coords[2])
+            neighborCoords[2] = grid->getFirstFluidNode(coords.data(), 2, grid->getStartZ());
+        else
+            neighborCoords[2] = grid->getLastFluidNode(coords.data(), 2, grid->getEndZ());
+    }
 }
 
 uint Side::getIndex(SPtr<Grid> grid, std::string coord, real constant, real v1, real v2)

src/gpu/GridGenerator/grid/BoundaryConditions/Side.h

@@ -33,10 +33,14 @@
 #ifndef SIDE_H
 #define SIDE_H
 
+#include <cstddef>
 #include <string>
 #include <vector>
+#include <map>
+#include <array>
 
 #include "gpu/GridGenerator/global.h"
+#include "lbm/constants/D3Q27.h"
 
 #define X_INDEX 0
 #define Y_INDEX 1
@@ -59,8 +63,6 @@ enum class SideType
     MX, PX, MY, PY, MZ, PZ, GEOMETRY
 };
 
-
-
 class Side
 {
 public:
@@ -72,7 +74,7 @@ public:
 
     virtual SideType whoAmI() const = 0;
 
-    std::vector<real> getNormal();
+    std::array<real, 3> getNormal() const;
 
 protected:
     void addIndices(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, std::string coord, real constant,
@@ -84,8 +86,29 @@ protected:
 
     void setQs(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint index, real q);
 
+    virtual void correctNeighborForPeriodicBoundaries(const Grid *grid, std::array<real, 3>& coords, std::array<real, 3>& neighbors) const;
+
+    virtual bool isAlignedWithMyNormal(const Grid *grid, int dir) const;
+    bool isAlignedWithNormal(const Grid *grid, int dir, const std::array<real, 3>& normal) const;
+
 private:
     static uint getIndex(SPtr<Grid> grid, std::string coord, real constant, real v1, real v2);
+    void resetDiagonalsInCaseOfOtherBC(Grid *grid, std::vector<real>& qNode, int dir, const std::array<real, 3> &coordinates) const;
+    std::array<real, 3> getNeighborCoordinates(Grid *grid, const std::array<real, 3> &coordinates,
+                                               size_t direction) const;
+    bool neighborNormalToSideIsAStopper(Grid *grid, const std::array<real, 3> &coordinates, SideType side) const;
+
+protected:
+    const std::map<SideType, const std::array<real, 3>> normals = {
+        { SideType::MX, { NEGATIVE_DIR, 0.0, 0.0 } }, { SideType::PX, { POSITIVE_DIR, 0.0, 0.0 } },
+        { SideType::MY, { 0.0, NEGATIVE_DIR, 0.0 } }, { SideType::PY, { 0.0, POSITIVE_DIR, 0.0 } },
+        { SideType::MZ, { 0.0, 0.0, NEGATIVE_DIR } }, { SideType::PZ, { 0.0, 0.0, POSITIVE_DIR } }
+    };
+    const std::map<SideType, size_t> sideToD3Q27 = {
+        { SideType::MX, vf::lbm::dir::DIR_M00 }, { SideType::PX, vf::lbm::dir::DIR_P00 },
+        { SideType::MY, vf::lbm::dir::DIR_0M0 }, { SideType::PY, vf::lbm::dir::DIR_0P0 },
+        { SideType::MZ, vf::lbm::dir::DIR_00M }, { SideType::PZ, vf::lbm::dir::DIR_00P }
+    };
 };
 
 class Geometry : public Side

src/gpu/GridGenerator/grid/Grid.h

@@ -47,6 +47,7 @@ struct Triangle;
 class GridInterface;
 class Object;
 class BoundingBox;
+enum class SideType;
 
 class GRIDGENERATOR_EXPORT Grid
 {
@@ -84,6 +85,8 @@ public:
     virtual void getGridInterfaceIndices(uint* iCellCfc, uint* iCellCff, uint* iCellFcc, uint* iCellFcf) const = 0;
     virtual bool isSparseIndexInFluidNodeIndicesBorder(uint &sparseIndex) const = 0;
 
+    virtual bool isStopperForBC(uint index) const = 0;
+
     virtual int *getNeighborsX() const = 0;
     virtual int *getNeighborsY() const = 0;
     virtual int *getNeighborsZ() const = 0;
@@ -133,9 +136,9 @@ public:
     virtual void setPeriodicityY(bool periodicity) = 0;
     virtual void setPeriodicityZ(bool periodicity) = 0;
 
-    virtual bool getPeriodicityX() = 0;
-    virtual bool getPeriodicityY() = 0;
-    virtual bool getPeriodicityZ() = 0;
+    virtual bool getPeriodicityX() const = 0;
+    virtual bool getPeriodicityY() const = 0;
+    virtual bool getPeriodicityZ() const = 0;
 
     virtual void setEnableFixRefinementIntoTheWall(bool enableFixRefinementIntoTheWall) = 0;
 
@@ -170,6 +173,11 @@ public:
 
     virtual void repairCommunicationIndices(int direction) = 0;
 
+    virtual bool nodeHasBC(uint index) const = 0;
+
+    virtual std::vector<SideType> getBCAlreadySet() = 0;
+    virtual void addBCalreadySet(SideType side) = 0;
+
     // needed for CUDA Streams 
     virtual void findFluidNodeIndices(bool onlyBulk) = 0;
     virtual uint getNumberOfFluidNodes() const = 0;
@@ -192,7 +200,6 @@ public:
     virtual void getFluidNodeIndicesMacroVars(uint *fluidNodeIndicesMacroVars) const = 0;
     virtual void getFluidNodeIndicesApplyBodyForce(uint *fluidNodeIndicesApplyBodyForce) const = 0;
     virtual void getFluidNodeIndicesAllFeatures(uint *fluidNodeIndicesAllFeatures) const = 0;
-
 };
 
 #endif

src/gpu/GridGenerator/grid/GridImp.cpp

@@ -720,6 +720,12 @@ void GridImp::setNonStopperOutOfGridCellTo(uint index, char type)
     }
 }
 
+bool GridImp::nodeHasBC(uint index) const
+{
+    return (getFieldEntry(index) == vf::gpu::BC_PRESSURE || getFieldEntry(index) == vf::gpu::BC_VELOCITY ||
+            getFieldEntry(index) == vf::gpu::BC_NOSLIP   || getFieldEntry(index) == vf::gpu::BC_SLIP     ||
+            getFieldEntry(index) == vf::gpu::BC_STRESS);
+}
 
 void GridImp::setPeriodicity(bool periodicityX, bool periodicityY, bool periodicityZ)
 {
@@ -743,17 +749,17 @@ void GridImp::setPeriodicityZ(bool periodicity)
     this->periodicityZ = periodicity;
 }
 
-bool GridImp::getPeriodicityX()
+bool GridImp::getPeriodicityX() const
 {
     return this->periodicityX;
 }
 
-bool GridImp::getPeriodicityY()
+bool GridImp::getPeriodicityY() const
 {
     return this->periodicityY;
 }
 
-bool GridImp::getPeriodicityZ()
+bool GridImp::getPeriodicityZ() const
 {
     return this->periodicityZ;
 }
@@ -2208,8 +2214,14 @@ void GridImp::getFluidNodeIndicesAllFeatures(uint *_fluidNodeIndicesAllFeatures)
 }
 
 
+std::vector<SideType> GridImp::getBCAlreadySet() {
+    return this->bcAlreadySet;
+}
 
-
+void GridImp::addBCalreadySet(SideType side)
+{
+    this->bcAlreadySet.push_back(side);
+}
 
 
 void GridImp::print() const
@@ -2219,3 +2231,10 @@ void GridImp::print() const
     if(this->gridInterface)
         this->gridInterface->print();
 }
+
+bool GridImp::isStopperForBC(uint index) const
+{
+    return (this->getFieldEntry(index) == vf::gpu::STOPPER_OUT_OF_GRID_BOUNDARY ||
+            this->getFieldEntry(index) == vf::gpu::STOPPER_OUT_OF_GRID ||
+            this->getFieldEntry(index) == vf::gpu::STOPPER_SOLID);
+}

src/gpu/GridGenerator/grid/GridImp.h

@@ -34,6 +34,7 @@
 #define GRID_IMP_H
 
 #include <array>
+#include <vector>
 
 #include "Core/LbmOrGks.h"
 
@@ -52,6 +53,7 @@ class Object;
 class BoundingBox;
 class TriangularMeshDiscretizationStrategy;
 
+
 #ifdef __GNUC__
     #ifndef __clang__
         #pragma push
@@ -76,7 +78,7 @@ protected:
 
 public:
     static SPtr<GridImp> makeShared(Object* object, real startX, real startY, real startZ, real endX, real endY, real endZ, real delta, std::string d3Qxx, uint level);
-    virtual ~GridImp() = default;
+    ~GridImp() override = default;
 
 private:
     void initalNumberOfNodesAndSize();
@@ -92,6 +94,7 @@ private:
     bool nodeInPreviousCellIs(int index, char type) const;
     bool nodeInCellIs(Cell& cell, char type) const override;
 
+
     uint getXIndex(real x) const;
     uint getYIndex(real y) const;
     uint getZIndex(real z) const;
@@ -135,6 +138,8 @@ private:
 
     bool enableFixRefinementIntoTheWall;
 
+    std::vector<SideType> bcAlreadySet;
+
 protected:
     Field field;
     int *neighborIndexX, *neighborIndexY, *neighborIndexZ, *neighborIndexNegative;
@@ -149,9 +154,9 @@ public:
     void setPeriodicityY(bool periodicity) override;
     void setPeriodicityZ(bool periodicity) override;
 
-    bool getPeriodicityX() override;
-    bool getPeriodicityY() override;
-    bool getPeriodicityZ() override;
+    bool getPeriodicityX() const override;
+    bool getPeriodicityY() const override;
+    bool getPeriodicityZ() const override;
 
     void setEnableFixRefinementIntoTheWall(bool enableFixRefinementIntoTheWall) override;
 
@@ -185,6 +190,9 @@ public:
 
     void setNumberOfLayers(uint numberOfLayers) override;
 
+    std::vector<SideType> getBCAlreadySet() override;
+    void addBCalreadySet(SideType side) override;
+
 public:
     Distribution distribution;
 
@@ -219,6 +227,7 @@ public:
     bool nodeInNextCellIs(int index, char type) const;
     bool hasAllNeighbors(uint index) const;
     bool hasNeighborOfType(uint index, char type) const;
+    bool nodeHasBC(uint index) const override;
     bool cellContainsOnly(Cell &cell, char type) const;
     bool cellContainsOnly(Cell &cell, char typeA, char typeB) const;
 
@@ -259,6 +268,8 @@ public:
     static void getGridInterface(uint *gridInterfaceList, const uint *oldGridInterfaceList, uint size);
 
     bool isSparseIndexInFluidNodeIndicesBorder(uint &sparseIndex) const override;
+    
+    bool isStopperForBC(uint index) const override;
 
     int *getNeighborsX() const override;
     int* getNeighborsY() const override;
@@ -276,7 +287,7 @@ public:
     void print() const;
 
 public:
-    virtual void findSparseIndices(SPtr<Grid> fineGrid) override;
+    void findSparseIndices(SPtr<Grid> fineGrid) override;
 
     void findForGridInterfaceNewIndices(SPtr<GridImp> fineGrid);
     void updateSparseIndices();

src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp

@@ -933,7 +933,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             //preprocessing
             real* QQ = para->getParH(i)->pressureBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->pressureBC.numberOfBCnodes;
-            QforBoundaryConditions Q;
+            QforBoundaryConditions &Q = para->getParH(i)->pressureBC;
             getPointersToBoundaryConditions(Q, QQ, sizeQ);
 
             builder->getPressureQs(Q.q27, i);
@@ -980,7 +980,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             //preprocessing
             real* QQ = para->getParH(i)->slipBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->slipBC.numberOfBCnodes;
-            QforBoundaryConditions Q;
+            QforBoundaryConditions &Q = para->getParH(i)->slipBC;
             getPointersToBoundaryConditions(Q, QQ, sizeQ);
 
             builder->getSlipQs(Q.q27, i);
@@ -1000,7 +1000,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             //preprocessing
             real* QQ = para->getParH(i)->stressBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->stressBC.numberOfBCnodes;
-            QforBoundaryConditions Q;
+            QforBoundaryConditions &Q = para->getParH(i)->stressBC;
             getPointersToBoundaryConditions(Q, QQ, sizeQ);
 
             builder->getStressQs(Q.q27, i);
@@ -1020,7 +1020,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             //preprocessing
             real* QQ = para->getParH(i)->velocityBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->velocityBC.numberOfBCnodes;
-            QforBoundaryConditions Q;
+            QforBoundaryConditions &Q = para->getParH(i)->velocityBC;
             getPointersToBoundaryConditions(Q, QQ, sizeQ);
             builder->getVelocityQs(Q.q27, i);
 
@@ -1120,7 +1120,7 @@ void GridGenerator::allocArrays_BoundaryQs()
             //preprocessing
             real* QQ = para->getParH(i)->geometryBC.q27[0];
             unsigned int sizeQ = para->getParH(i)->geometryBC.numberOfBCnodes;
-            QforBoundaryConditions Q;
+            QforBoundaryConditions &Q = para->getParH(i)->geometryBC;
             getPointersToBoundaryConditions(Q, QQ, sizeQ);
             //////////////////////////////////////////////////////////////////
 
@@ -1362,31 +1362,31 @@ std::string GridGenerator::checkNeighbor(int level, real x, real y, real z, int
 }
 
 void GridGenerator::getPointersToBoundaryConditions(QforBoundaryConditions& boundaryConditionStruct, real* subgridDistances, const unsigned int numberOfBCnodes){
-    boundaryConditionStruct.q27[DIR_P00] =    &subgridDistances[DIR_P00   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_M00] =    &subgridDistances[DIR_M00   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0P0] =    &subgridDistances[DIR_0P0   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0M0] =    &subgridDistances[DIR_0M0   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_00P] =    &subgridDistances[DIR_00P   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_00M] =    &subgridDistances[DIR_00M   * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PP0] =   &subgridDistances[DIR_PP0  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MM0] =   &subgridDistances[DIR_MM0  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PM0] =   &subgridDistances[DIR_PM0  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MP0] =   &subgridDistances[DIR_MP0  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_P0P] =   &subgridDistances[DIR_P0P  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_M0M] =   &subgridDistances[DIR_M0M  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_P0M] =   &subgridDistances[DIR_P0M  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_M0P] =   &subgridDistances[DIR_M0P  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0PP] =   &subgridDistances[DIR_0PP  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0MM] =   &subgridDistances[DIR_0MM  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0PM] =   &subgridDistances[DIR_0PM  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_0MP] =   &subgridDistances[DIR_0MP  * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_000] = &subgridDistances[DIR_000* numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PPP] =  &subgridDistances[DIR_PPP * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MMP] =  &subgridDistances[DIR_MMP * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PMP] =  &subgridDistances[DIR_PMP * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MPP] =  &subgridDistances[DIR_MPP * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PPM] =  &subgridDistances[DIR_PPM * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MMM] =  &subgridDistances[DIR_MMM * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_PMM] =  &subgridDistances[DIR_PMM * numberOfBCnodes];
-    boundaryConditionStruct.q27[DIR_MPM] =  &subgridDistances[DIR_MPM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_P00] = &subgridDistances[DIR_P00 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_M00] = &subgridDistances[DIR_M00 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0P0] = &subgridDistances[DIR_0P0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0M0] = &subgridDistances[DIR_0M0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_00P] = &subgridDistances[DIR_00P * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_00M] = &subgridDistances[DIR_00M * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PP0] = &subgridDistances[DIR_PP0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MM0] = &subgridDistances[DIR_MM0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PM0] = &subgridDistances[DIR_PM0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MP0] = &subgridDistances[DIR_MP0 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_P0P] = &subgridDistances[DIR_P0P * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_M0M] = &subgridDistances[DIR_M0M * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_P0M] = &subgridDistances[DIR_P0M * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_M0P] = &subgridDistances[DIR_M0P * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0PP] = &subgridDistances[DIR_0PP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0MM] = &subgridDistances[DIR_0MM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0PM] = &subgridDistances[DIR_0PM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_0MP] = &subgridDistances[DIR_0MP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_000] = &subgridDistances[DIR_000 * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PPP] = &subgridDistances[DIR_PPP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MMP] = &subgridDistances[DIR_MMP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PMP] = &subgridDistances[DIR_PMP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MPP] = &subgridDistances[DIR_MPP * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PPM] = &subgridDistances[DIR_PPM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MMM] = &subgridDistances[DIR_MMM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_PMM] = &subgridDistances[DIR_PMM * numberOfBCnodes];
+    boundaryConditionStruct.q27[DIR_MPM] = &subgridDistances[DIR_MPM * numberOfBCnodes];
 }

src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriter.hpp

+#ifndef QVTKWRITER_HPP
+#define QVTKWRITER_HPP
+
+#include <array>
+#include <vector>
+
+#include "basics/Core/StringUtilities/StringUtil.h"
+#include "basics/utilities/UbSystem.h"
+#include "basics/writer/WbWriterVtkXmlBinary.h"
+#include "lbm/constants/D3Q27.h"
+#include "logger/Logger.h"
+
+#include "gpu/GridGenerator/grid/NodeValues.h"
+#include "gpu/VirtualFluids_GPU/Communication/Communicator.h"
+#include "gpu/VirtualFluids_GPU/LBM/LB.h"
+#include "gpu/VirtualFluids_GPU/Parameter/Parameter.h"
+#include "gpu/VirtualFluids_GPU/Utilities/FindNeighbors.h"
+
+namespace QDebugVtkWriter
+{
+
+using namespace vf::lbm::dir;
+
+namespace
+{
+inline void modifyLineLengthsForQs(const std::array<double, 3> &coords, std::array<double, 3> &neighborCoords, real q)
+{
+    if (q == 1.0 || q <= 0.0)
+        return;
+
+    const auto dx = neighborCoords[0] - coords[0];
+    const auto dy = neighborCoords[1] - coords[1];
+    const auto dz = neighborCoords[2] - coords[2];
+
+    neighborCoords[0] = coords[0] + q * dx;
+    neighborCoords[1] = coords[1] + q * dy;
+    neighborCoords[2] = coords[2] + q * dz;
+}
+
+inline void writeQLines(LBMSimulationParameter *parH, QforBoundaryConditions &boundaryQ, const std::string &filepath,
+                        WbWriter *writer)
+{
+    VF_LOG_INFO("Write qs in for boundary condition to {}.", filepath);
+
+    const auto numberOfNodes = boundaryQ.numberOfBCnodes;
+    std::vector<UbTupleFloat3> nodes;
+    nodes.reserve(numberOfNodes * 8 * 2);
+    std::vector<UbTupleInt2> lines;
+    lines.reserve(numberOfNodes * 8);
+
+    std::vector<std::string> dataNames = { "nodeIndex", "q" };
+    std::vector<std::vector<float>> lineData(2);
+
+    for (size_t i = 0; i < numberOfNodes; i++) {
+        const auto nodeIndex = boundaryQ.k[i];
+        const std::array<double, 3> coords = { parH->coordinateX[nodeIndex], parH->coordinateY[nodeIndex],
+                                               parH->coordinateZ[nodeIndex] };
+
+        for (size_t direction = 1; direction < ENDDIR; direction++) {
+
+            const auto q = boundaryQ.q27[direction][i];
+            if (q <= (real)0.0) {
+                continue;
+            }
+
+            const auto positionNeighbor = getNeighborIndex(parH, (uint)nodeIndex, (int)direction);
+
+            std::array<double, 3> neighborCoords = { parH->coordinateX[positionNeighbor],
+                                                     parH->coordinateY[positionNeighbor],
+                                                     parH->coordinateZ[positionNeighbor] };
+
+            modifyLineLengthsForQs(coords, neighborCoords, q);
+
+            nodes.emplace_back(float(coords[0]), float(coords[1]), coords[2]);
+            nodes.emplace_back(float(neighborCoords[0]), float(neighborCoords[1]), float(neighborCoords[2]));
+
+            lines.emplace_back((int)nodes.size() - 2, (int)nodes.size() - 1);
+            lineData[0].push_back(nodeIndex);
+            lineData[1].push_back(q);
+        }
+    }
+
+    writer->writeLinesWithLineData(filepath, nodes, lines, dataNames, lineData);
+}
+} // namespace
+
+inline void writeQLinesDebug(Parameter *para, QforBoundaryConditions &boundaryQ, uint level, std::string fileName)
+{
+    const auto filePath = para->getFName() + "_" + fileName + ".vtk";
+    auto writer = WbWriterVtkXmlBinary::getInstance();
+    writeQLines(para->getParH(level).get(), boundaryQ, filePath, writer);
+}
+
+} // namespace QDebugVtkWriter
+
+#endif

src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriterTest.cpp

+#include "gmock/gmock.h"
+#include <cmath>
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+#include "QDebugVtkWriter.hpp"
+#include <tuple>
+
+MATCHER(DoubleNear5, "") {
+    return abs(std::get<0>(arg) - std::get<1>(arg)) < 0.00001;
+}
+
+using namespace QDebugVtkWriter;
+
+double calcVectorLength(const std::array<double, 3> coords, const std::array<double, 3> neighborCoords)
+{
+    return std::sqrt(std::pow((neighborCoords[0] - coords[0]), 2) + std::pow((neighborCoords[1] - coords[1]), 2) +
+                     std::pow((neighborCoords[2] - coords[2]), 2));
+}
+
+TEST(QDebugVtkWriterTest, modifyLineLengthsForQsSameCoords3)
+{
+    const std::array<double, 3> coords = { 0, 0, 0 };
+    std::array<double, 3> neighborCoords = { 1, 1, 1 };
+    const real q = 0.3;
+    const real initialLength = calcVectorLength(coords, neighborCoords);
+
+    modifyLineLengthsForQs(coords, neighborCoords, q);
+
+    std::array<double, 3> expectedNeighborCoords = { 0.3, 0.3, 0.3 };
+    EXPECT_THAT(neighborCoords,testing::Pointwise(DoubleNear5(), expectedNeighborCoords));
+    EXPECT_THAT(calcVectorLength(coords, neighborCoords), testing::DoubleNear(q*initialLength, 0.00001));
+}
+
+TEST(QDebugVtkWriterTest, modifyLineLengthDifferentCoords)
+{
+    const std::array<double, 3> coords = { 0, 0, 0 };
+    std::array<double, 3> neighborCoords = { 1, 2, 3 };
+    const real q = 0.3;
+    const real initialLength = calcVectorLength(coords, neighborCoords);
+
+    modifyLineLengthsForQs(coords, neighborCoords, q);
+
+    std::array<double, 3> expectedNeighborCoords = { 0.3, 0.6, 0.9 };
+    EXPECT_THAT(neighborCoords,testing::Pointwise(DoubleNear5(), expectedNeighborCoords));
+    EXPECT_THAT(calcVectorLength(coords, neighborCoords), testing::DoubleNear(q*initialLength, 0.00001));
+}
+
+TEST(QDebugVtkWriterTest, modifyLineLengthNegativeCoord)
+{
+    const std::array<double, 3> coords = { 0, 0, 0 };
+    std::array<double, 3> neighborCoords = { 1, 2, -3 };
+    const real q = 0.3;
+    const real initialLength = calcVectorLength(coords, neighborCoords);
+
+    modifyLineLengthsForQs(coords, neighborCoords, q);
+
+    std::array<double, 3> expectedNeighborCoords = { 0.3, 0.6, -0.9 };
+    EXPECT_THAT(neighborCoords,testing::Pointwise(DoubleNear5(), expectedNeighborCoords));
+    EXPECT_THAT(calcVectorLength(coords, neighborCoords), testing::DoubleNear(q*initialLength, 0.00001));
+}
\ No newline at end of file

src/gpu/VirtualFluids_GPU/Output/QDebugWriter.hpp

@@ -13,8 +13,6 @@
 #include <basics/writer/WbWriterVtkXmlBinary.h>
 #include "Core/StringUtilities/StringUtil.h"
 
-//using namespace std;
-
 namespace QDebugWriter
 {
     void writeQValues(QforBoundaryConditions &Q, int* k, int kq, const std::string &name)