diff --git a/CITATION.cff b/CITATION.cff
new file mode 100644
index 0000000000000000000000000000000000000000..50d4989d5c269521392644515d716fa93b3cf6e3
--- /dev/null
+++ b/CITATION.cff
@@ -0,0 +1,40 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+type: software
+authors:
+  - family-names: Kutscher
+    given-names: Konstantin
+    orcid: https://orcid.org/0000-0002-1099-1608
+  - family-names: Schönherr
+    given-names: Martin
+    orcid: https://orcid.org/0000-0002-4774-1776
+  - family-names: Geier
+    given-names: Martin
+    orcid: https://orcid.org/0000-0002-8367-9412
+  - family-names: Krafczyk
+    given-names: Manfred
+    orcid: https://orcid.org/0000-0002-8509-0871
+  - family-names: Alihussein
+    given-names: Hussein
+    orcid: https://orcid.org/0000-0003-3656-7028
+  - family-names: Linxweiler
+    given-names: Jan
+    orcid: https://orcid.org/0000-0002-2755-5087
+  - family-names: Peters
+    given-names: Sören
+    orcid: https://orcid.org/0000-0001-5236-3776
+  - family-names: Wellmann
+    given-names: Anna
+    orcid: https://orcid.org/0000-0002-8825-2995
+  - family-names: Safari
+    given-names: Hesameddin
+    orcid: https://orcid.org/0000-0002-2755-5087
+  - family-names: Marcus
+    given-names: Sven
+    orcid: https://orcid.org/0000-0003-3689-2162
+title: "VirtualFluids"
+version: 0.1.0
+license: GPL-3.0-or-later
+repository-code: "https://git.rz.tu-bs.de/irmb/VirtualFluids"
+date-released: "XXXXXXX"
+
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b3af407acd66ec3223f55de7753df879786ce561..c6498bf19bb021f3ae19d69c4131aa56476149be 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,7 +9,7 @@
 cmake_minimum_required(VERSION 3.15..3.20 FATAL_ERROR)
 
 project(VirtualFluids
-        VERSION 1.0.0
+        VERSION 0.1.0
         DESCRIPTION "CFD code based on the Lattice Boltzmann Method"
         HOMEPAGE_URL "https://www.tu-braunschweig.de/irmb/forschung/virtualfluids"
         LANGUAGES CXX)
diff --git a/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp b/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp
index ed2bd3a9c5be2c7847799ad0a207f8900230eba7..90fcc7ec2e0b426eb54ac91b472ad5af631e497b 100644
--- a/apps/gpu/LBM/DrivenCavity/DrivenCavity.cpp
+++ b/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;
 
diff --git a/apps/gpu/LBM/DrivenCavityUniform/CMakeLists.txt b/apps/gpu/LBM/DrivenCavityUniform/CMakeLists.txt
index 8384e1bc6fcfa3fd2514434b620b266e96b3626a..40b4f08d7500c56efae7378df6398d065e4ecbfb 100644
--- a/apps/gpu/LBM/DrivenCavityUniform/CMakeLists.txt
+++ b/apps/gpu/LBM/DrivenCavityUniform/CMakeLists.txt
@@ -1,4 +1,4 @@
-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)
+
diff --git a/metadata.xml b/metadata.xml
deleted file mode 100644
index 7cbae3ae7e1d5d7d48af2f0e5577253a89f953f5..0000000000000000000000000000000000000000
--- a/metadata.xml
+++ /dev/null
@@ -1,204 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<resource xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://datacite.org/schema/kernel-4" xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4.3/metadata.xsd">
-	<identifier identifierType="DOI">PLACEHOLDER</identifier>
-	<titles>
-		<title xml:lang="en">VirtualFluids</title>
-	</titles>
-	<language>en</language>
-	<creators>
-		<creator>
-			<creatorName nameType="Personal">Krafczyk, Manfred</creatorName>
-			<givenName>Manfred</givenName>
-			<familyName>Krafczyk</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID">0000-0002-8509-0871</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="de">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</creator>
-		<creator>
-			<creatorName nameType="Organizational">Institut für rechnergestützte Modellierung im Bauingenieurwesen</creatorName>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-		</creator>
-	</creators>
-	<publisher xml:lang="de">Institut für rechnergestützte Modellierung im Bauingenieurwesen</publisher>
-	<publicationYear>2021</publicationYear>
-	<resourceType resourceTypeGeneral="Software">Computational Fluid Dynamics Solver</resourceType>
-	<subjects>
-		<subject subjectScheme="DDC" schemeURI="https://www.oclc.org/en/dewey.html">532 Fluid Mechanics, liquid mechanics</subject>
-	</subjects>
-	<contributors>
-		<contributor contributorType="Researcher">
-			<contributorName>Ahrenholz, Benjamin</contributorName>
-			<givenName>Benjamin</givenName>
-			<familyName>Ahrenholz</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Alihussein, Hussein</contributorName>
-			<givenName>Hussein</givenName>
-			<familyName>Alihussein</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0003-3656-7028</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Bindick, Sebastian</contributorName>
-			<givenName>Sebastian</givenName>
-			<familyName>Bindick</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Brendel, Aileen</contributorName>
-			<givenName>Aileen</givenName>
-			<familyName>Brendel</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Geier, Martin</contributorName>
-			<givenName>Martin</givenName>
-			<familyName>Geier</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0002-8367-9412</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Geller, Sebastian</contributorName>
-			<givenName>Sebastian</givenName>
-			<familyName>Geller</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Goraki Fard, Ehsan</contributorName>
-			<givenName>Ehsan</givenName>
-			<familyName>Goraki Fard</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Hegewald, Jan</contributorName>
-			<givenName>Jan</givenName>
-			<familyName>Hegewald</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Janßen, Christian</contributorName>
-			<givenName>Christian</givenName>
-			<familyName>Janßen</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Kutscher, Konstantin</contributorName>
-			<givenName>Konstantin</givenName>
-			<familyName>Kutscher</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0002-1099-1608</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Lenz, Stephan</contributorName>
-			<givenName>Stephan</givenName>
-			<familyName>Lenz</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Linxweiler, Jan</contributorName>
-			<givenName>Jan</givenName>
-			<familyName>Linxweiler</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0002-2755-5087</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Lux, Lennard</contributorName>
-			<givenName>Lennard</givenName>
-			<familyName>Lux</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Marcus, Sven</contributorName>
-			<givenName>Sven</givenName>
-			<familyName>Marcus</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0003-3689-2162</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Universitätsbibliothek Braunschweig</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Peters, Sören</contributorName>
-			<givenName>Sören</givenName>
-			<familyName>Peters</familyName>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Safari, Hesameddin</contributorName>
-			<givenName>Hesameddin</givenName>
-			<familyName>Safari</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Schönherr, Martin</contributorName>
-			<givenName>Martin</givenName>
-			<familyName>Schönherr</familyName>
-			<nameIdentifier nameIdentifierScheme="ORCID" schemeURI="http://orcid.org/">0000-0002-4774-1776</nameIdentifier>
-			<affiliation xml:lang="de">TU Braunschweig</affiliation>
-			<affiliation xml:lang="en">Institut für rechnergestützte Modellierung im Bauingenieurwesen</affiliation>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Stiebler, Maik</contributorName>
-			<givenName>Maik</givenName>
-			<familyName>Stiebler</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Textor, Sören</contributorName>
-			<givenName>Sören</givenName>
-			<familyName>Textor</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Tölke, Jonas</contributorName>
-			<givenName>Jonas</givenName>
-			<familyName>Tölke</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Uphoff, Sonja</contributorName>
-			<givenName>Sonja</givenName>
-			<familyName>Uphoff</familyName>
-		</contributor>
-
-		<contributor contributorType="Researcher">
-			<contributorName>Wellmann, Anna</contributorName>
-			<givenName>Anna</givenName>
-			<familyName>Wellmann</familyName>
-		</contributor>
-	</contributors>
-	<dates>
-		<date dateType="Created">2000</date>
-	</dates>
-	<formats>
-		<format>text/x-c</format>
-		<format>text/x-h</format>
-		<format>text/x-script.python</format>
-	</formats>
-	<relatedIdentifiers>
-		<relatedIdentifier relatedIdentifierType="URL" relationType="Requires" resourceTypeGeneral="Software">https://www.open-mpi.org/software/ompi/v4.1/</relatedIdentifier>
-		<relatedIdentifier relatedIdentifierType="URL" relationType="IsCompiledBy" resourceTypeGeneral="Software">https://cmake.org</relatedIdentifier>
-		<relatedIdentifier relatedIdentifierType="URL" relationType="IsCompiledBy" resourceTypeGeneral="Software">https://gcc.gnu.org</relatedIdentifier>
-		<relatedIdentifier relatedIdentifierType="URL" relationType="IsCompiledBy" resourceTypeGeneral="Software">https://clang.llvm.org</relatedIdentifier>
-		<relatedIdentifier relatedIdentifierType="URL" relationType="IsCompiledBy" resourceTypeGeneral="Software">https://visualstudio.microsoft.com/vs/features/cplusplus/</relatedIdentifier>
-	</relatedIdentifiers>
-	<rightsList>
-		<rights xml:lang="en" schemeURI="https://spdx.org/licenses/" rightsIdentifierScheme="SPDX" rightsIdentifier="GPL-3.0-only" rightsURI="https://www.gnu.org/licenses/gpl-3.0-standalone.html">GNU General Public License Version 3</rights>
-	</rightsList>
-	<descriptions>
-		<description descriptionType="Abstract">
-			VirtualFluids (VF) is a research code developed at the Institute for Computational Modeling in Civil Engineering (iRMB). The code is a Computational Fluid Dynamics (CFD) solver based on the Lattice Boltzmann Method (LBM) for turbulent, thermal, multiphase and multicomponent flow problems as well as for multi-field problems such as Fluid-Structure-interaction including distributed pre- and postprocessing capabilities for simulations with more than 100 billion degrees of freedom.
-		</description>
-	</descriptions>
-</resource>
diff --git a/regression-tests/driven_cavity_test.sh b/regression-tests/driven_cavity_test.sh
index a00e550f0d00181d6e6ab00fc97a619ad8301c38..7f799facb4459ddafcd8b210a5477954af1444cb 100755
--- a/regression-tests/driven_cavity_test.sh
+++ b/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
diff --git a/regression-tests/driven_cavity_uniform_test.sh b/regression-tests/driven_cavity_uniform_test.sh
index 9337eb7f9e607a88e884728b599e863f5b746a63..95e2bab635d3a6a73fb514a1f67902083c98e5d3 100755
--- a/regression-tests/driven_cavity_uniform_test.sh
+++ b/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"
diff --git a/setup.cfg b/setup.cfg
index f3bc36f29efe64b81916efa04b9d7831a2abeb09..5894f9dec06953c3eeb909af96db9cb19d202d65 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -5,7 +5,7 @@ long_description = file: README.md
 long_description_content_type = text/markdown
 platforms = any
 url = https://git.rz.tu-bs.de/irmb/virtualfluids
-version = 0.0.1
+version = 0.1.0
 
 [options]
 python_requires = >=3.6
diff --git a/src/basics/basics/writer/WbWriter.h b/src/basics/basics/writer/WbWriter.h
index 26d43464c03311a2cbc14cd4fc9fe717d4b01531..55dceb7cb4a64dc90f0677796cab52135b726f56 100644
--- a/src/basics/basics/writer/WbWriter.h
+++ b/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:
diff --git a/src/basics/basics/writer/WbWriterVtkXmlBinary.cpp b/src/basics/basics/writer/WbWriterVtkXmlBinary.cpp
index 6731fa56026ca284ad671cb6ce59000a609bbb8c..55c3541983ea4248512508146792832a34a1c563 100644
--- a/src/basics/basics/writer/WbWriterVtkXmlBinary.cpp
+++ b/src/basics/basics/writer/WbWriterVtkXmlBinary.cpp
@@ -34,6 +34,8 @@
 #include <basics/writer/WbWriterVtkXmlASCII.h>
 #include <basics/writer/WbWriterVtkXmlBinary.h>
 #include <cstring>
+#include <fstream>
+#include <string>
 
 using namespace std;
 
@@ -154,12 +156,13 @@ string WbWriterVtkXmlBinary::writeParallelFile(const string &filename, vector<st
 
     return vtkfilename;
 }
+
 /*===============================================================================*/
-string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFloat3> &nodes,
-                                        vector<UbTupleInt2> &lines)
+
+// helper functions
+
+ofstream createFileStream(std::string vtkfilename)
 {
-    string vtkfilename = filename + getFileExtension();
-    UBLOG(logDEBUG1, "WbWriterVtkXmlBinary::writeLines to " << vtkfilename << " - start");
 
     ofstream out(vtkfilename.c_str(), ios::out | ios::binary);
     if (!out) {
@@ -172,89 +175,199 @@ string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFl
         if (!out)
             throw UbException(UB_EXARGS, "couldn't open file " + vtkfilename);
     }
+    return out;
+}
 
-    int nofNodes = (int)nodes.size();
-    int nofCells = (int)lines.size();
-
-    int bytesPerByteVal      = 4; //==sizeof(int)
-    int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 2 /*nodes per line */ * nofCells * sizeof(int);
-    int bytesCellOffsets     = 1 /*offset per line */ * nofCells * sizeof(int);
-    int bytesCellTypes       = 1 /*type of line */ * nofCells * sizeof(unsigned char);
-
-    int offset = 0;
-    // VTK FILE
+void writeVtkHeader(ofstream &out, int numberOfNodes, int numberOfCells)
+{
     out << "<?xml version=\"1.0\"?>\n";
     out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" >"
         << "\n";
     out << "   <UnstructuredGrid>"
         << "\n";
-    out << "      <Piece NumberOfPoints=\"" << nofNodes << "\" NumberOfCells=\"" << nofCells << "\">\n";
+    out << "      <Piece NumberOfPoints=\"" << numberOfNodes << "\" NumberOfCells=\"" << numberOfCells << "\">\n";
+}
 
-    // POINTS SECTION
+int writePointHeader(ofstream &out, int offset, int bytesPerByteVal, int bytesPoints)
+{
     out << "         <Points>\n";
     out << "            <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"appended\" offset=\"" << offset
         << "\"  />\n";
     out << "         </Points>\n";
     offset += (bytesPerByteVal + bytesPoints);
+    return offset;
+}
 
-    // CELLS SECTION
+int writeCellHeader(ofstream &out, int offset, int bytesPerByteVal, int bytesCellConnectivity, int bytesCellOffsets,
+                    int bytesCellTypes)
+{
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
     out << "            <DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\"" << offset << "\" />\n ";
     offset += (bytesPerByteVal + bytesCellTypes);
     out << "         </Cells>\n";
+    return offset;
+}
 
+int writeDataHeader(ofstream &out, vector<string> &datanames, int offset, int bytesPerByteVal, int bytesScalarData)
+{
+    out << "         <CellData>\n";
+    for (size_t s = 0; s < datanames.size(); ++s) {
+        out << "            <DataArray type=\"Float32\" Name=\"" << datanames[s] << "\" format=\"appended\" offset=\""
+            << offset << "\" /> \n";
+        offset += (bytesPerByteVal + bytesScalarData);
+    }
+    out << "         </CellData>\n";
+    return offset;
+}
+
+void writeAppendDataHeader(ofstream &out)
+{
     out << "      </Piece>\n";
     out << "   </UnstructuredGrid>\n";
-
-    // AppendedData SECTION
     out << "   <AppendedData encoding=\"raw\">\n";
     out << "_";
+}
 
-    // POINTS SECTION
+void writePoints(ofstream &out, int bytesPerByteVal, int bytesPoints, vector<UbTupleFloat3> &nodes)
+{
     out.write((char *)&bytesPoints, bytesPerByteVal);
-    for (int n = 0; n < nofNodes; n++) {
+    for (int n = 0; n < (int)nodes.size(); n++) {
         out.write((char *)&val<1>(nodes[n]), sizeof(float));
         out.write((char *)&val<2>(nodes[n]), sizeof(float));
         out.write((char *)&val<3>(nodes[n]), sizeof(float));
     }
+}
 
-    // CELLS SECTION
-    // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
-    for (int c = 0; c < nofCells; c++) {
-        out.write((char *)&val<1>(lines[c]), sizeof(int));
-        out.write((char *)&val<2>(lines[c]), sizeof(int));
+void writeCellConnectivity(ofstream &out, int bytesPerByteVal, int bytesCellConnectivity, vector<UbTupleInt2> &cells)
+{
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
+    for (int c = 0; c < (int)cells.size(); c++) {
+        out.write((char *)&val<1>(cells[c]), sizeof(int));
+        out.write((char *)&val<2>(cells[c]), sizeof(int));
     }
+}
 
-    // cellOffsets
+void writeCellOffsets(ofstream &out, int bytesPerByteVal, int bytesCellOffsets, int numberOfCells)
+{
     out.write((char *)&bytesCellOffsets, bytesPerByteVal);
     int itmp;
-    for (int c = 1; c <= nofCells; c++) {
+    for (int c = 1; c <= numberOfCells; c++) {
         itmp = 2 * c;
         out.write((char *)&itmp, sizeof(int));
     }
+}
 
-    // cellTypes
+void writeCellTypes(ofstream &out, int bytesPerByteVal, int bytesCellTypes, int numberOfCells)
+{
     out.write((char *)&bytesCellTypes, bytesPerByteVal);
     unsigned char vtkCellType = 3;
-    for (int c = 0; c < nofCells; c++) {
+    for (int c = 0; c < numberOfCells; c++) {
         out.write((char *)&vtkCellType, sizeof(unsigned char));
     }
+}
+
+void writeCellData(ofstream &out, int bytesPerByteVal, int bytesScalarData, vector<string> &datanames,
+                   vector<vector<float>> &celldata)
+{
+    for (size_t s = 0; s < datanames.size(); ++s) {
+        out.write((char *)&bytesScalarData, bytesPerByteVal);
+        for (size_t d = 0; d < celldata[s].size(); ++d) {
+            // loake kopie machen, da in celldata "doubles" sind
+            float tmp = (float)celldata[s][d];
+            out.write((char *)&tmp, sizeof(float));
+        }
+    }
+}
+
+void writeEndOfFile(ofstream &out)
+{
     out << "\n</AppendedData>\n";
     out << "</VTKFile>";
     out << endl;
     out.close();
+}
+
+/*===============================================================================*/
+string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFloat3> &nodes,
+                                        vector<UbTupleInt2> &lines)
+{
+    string vtkfilename = filename + getFileExtension();
+    UBLOG(logDEBUG1, "WbWriterVtkXmlBinary::writeLines to " << vtkfilename << " - start");
+
+    ofstream out = createFileStream(vtkfilename);
+
+    int nofNodes = (int)nodes.size();
+    int nofCells = (int)lines.size();
+
+    int bytesPerByteVal = 4; //==sizeof(int)
+    int bytesPoints = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
+    int bytesCellConnectivity = 2 /*nodes per line */ * nofCells * sizeof(int);
+    int bytesCellOffsets = 1 /*offset per line */ * nofCells * sizeof(int);
+    int bytesCellTypes = 1 /*type of line */ * nofCells * sizeof(unsigned char);
+
+    int offset = 0;
+
+    writeVtkHeader(out, nofNodes, nofCells);
+    offset = writePointHeader(out, offset, bytesPerByteVal, bytesPoints);
+    writeCellHeader(out, offset, bytesPerByteVal, bytesCellConnectivity, bytesCellOffsets, bytesCellTypes);
+    writeAppendDataHeader(out);
+
+    writePoints(out, bytesPerByteVal, bytesPoints, nodes);
+    writeCellConnectivity(out, bytesPerByteVal, bytesCellConnectivity, lines);
+    writeCellOffsets(out, bytesPerByteVal, bytesCellOffsets, nofCells);
+    writeCellTypes(out, bytesPerByteVal, bytesCellTypes, nofCells);
+    writeEndOfFile(out);
     UBLOG(logDEBUG1, "WbWriterVtkXmlBinary::writeLines to " << vtkfilename << " - end");
 
     return vtkfilename;
 }
+
+/*===============================================================================*/
+string WbWriterVtkXmlBinary::writeLinesWithLineData(const string &filename, vector<UbTupleFloat3> &nodes,
+                                                    vector<UbTupleInt2> &lines, vector<string> &datanames,
+                                                    vector<vector<float>> &celldata)
+{
+    string vtkfilename = filename + getFileExtension();
+    UBLOG(logDEBUG1, "WbWriterVtkXmlBinary::writeLinesWithLineData to " << vtkfilename << " - start");
+
+    ofstream out = createFileStream(vtkfilename);
+
+    int nofNodes = (int)nodes.size();
+    int nofCells = (int)lines.size();
+
+    int bytesPerByteVal = 4; //==sizeof(int)
+    int bytesPoints = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
+    int bytesCellConnectivity = 2 /*nodes per line */ * nofCells * sizeof(int);
+    int bytesCellOffsets = 1 /*offset per line */ * nofCells * sizeof(int);
+    int bytesCellTypes = 1 /*type of line */ * nofCells * sizeof(unsigned char);
+    int bytesScalarData = 1 /*scalar        */ * nofCells * sizeof(float);
+
+    int offset = 0;
+
+    writeVtkHeader(out, nofNodes, nofCells);
+    offset = writePointHeader(out, offset, bytesPerByteVal, bytesPoints);
+    offset = writeCellHeader(out, offset, bytesPerByteVal, bytesCellConnectivity, bytesCellOffsets, bytesCellTypes);
+    writeDataHeader(out, datanames, offset, bytesPerByteVal, bytesScalarData);
+    writeAppendDataHeader(out);
+
+    writePoints(out, bytesPerByteVal, bytesPoints, nodes);
+    writeCellConnectivity(out, bytesPerByteVal, bytesCellConnectivity, lines);
+    writeCellOffsets(out, bytesPerByteVal, bytesCellOffsets, nofCells);
+    writeCellTypes(out, bytesPerByteVal, bytesCellTypes, nofCells);
+    writeCellData(out, bytesPerByteVal, bytesScalarData, datanames, celldata);
+    writeEndOfFile(out);
+
+    UBLOG(logDEBUG1, "WbWriterVtkXmlBinary::writeLinesWithLineData to " << vtkfilename << " - end");
+
+    return vtkfilename;
+}
+
 /*===============================================================================*/
 // std::string WbWriterVtkXmlBinary::writeLinesWithNodeData(const string& filename,vector<UbTupleFloat3 >& nodes,
 // vector<UbTupleInt2 >& lines, std::vector< std::string >& datanames, std::vector< std::vector< double > >& nodedata)
@@ -276,7 +389,7 @@ string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFl
 //
 //   int bytesPerByteVal      = 4; //==sizeof(int)
 //   int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-//   int bytesCellConnectivty = 2 /*nodes per line  */ * nofCells * sizeof(int  );
+//   int bytesCellConnectivity = 2 /*nodes per line  */ * nofCells * sizeof(int  );
 //   int bytesCellOffsets     = 1 /*offset per line */ * nofCells * sizeof(int  );
 //   int bytesCellTypes       = 1 /*type of line    */ * nofCells * sizeof(unsigned char);
 //   int bytesScalarData      = 1 /*scalar          */ * nofNodes * sizeof(float);
@@ -296,7 +409,7 @@ string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFl
 //   //CELLS SECTION
 //   out<<"         <Cells>\n";
 //   out<<"            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\""<< offset <<"\"
-//   />\n"; offset += (bytesPerByteVal + bytesCellConnectivty); out<<"            <DataArray type=\"Int32\"
+//   />\n"; offset += (bytesPerByteVal + bytesCellConnectivity); out<<"            <DataArray type=\"Int32\"
 //   Name=\"offsets\" format=\"appended\" offset=\""<< offset <<"\" />\n"; offset += (bytesPerByteVal +
 //   bytesCellOffsets); out<<"            <DataArray type=\"UInt8\" Name=\"types\" format=\"appended\" offset=\""<<
 //   offset <<"\" />\n "; offset += (bytesPerByteVal + bytesCellTypes); out<<"         </Cells>\n";
@@ -328,7 +441,7 @@ string WbWriterVtkXmlBinary::writeLines(const string &filename, vector<UbTupleFl
 //
 //   //CELLS SECTION
 //   //cellConnectivity
-//   out.write( (char*)&bytesCellConnectivty, bytesPerByteVal );
+//   out.write( (char*)&bytesCellConnectivity, bytesPerByteVal );
 //   for(int c=0; c<nofCells; c++)
 //   {
 //      out.write( (char*)&val<1>(lines[c]), sizeof(int) );
@@ -397,7 +510,7 @@ string WbWriterVtkXmlBinary::writeTriangles(const string &filename, vector<UbTup
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3 - coord    */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 3 /*nodes per triangle  */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 3 /*nodes per triangle  */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per triangle */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of triangle    */ * nofCells * sizeof(unsigned char);
 
@@ -421,7 +534,7 @@ string WbWriterVtkXmlBinary::writeTriangles(const string &filename, vector<UbTup
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -446,7 +559,7 @@ string WbWriterVtkXmlBinary::writeTriangles(const string &filename, vector<UbTup
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(triangles[c]), sizeof(int));
         out.write((char *)&val<2>(triangles[c]), sizeof(int));
@@ -502,7 +615,7 @@ string WbWriterVtkXmlBinary::writeTrianglesWithNodeData(const string &filename,
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 3 /*nodes per tri   */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 3 /*nodes per tri   */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per tri  */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of tri     */ * nofCells * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar          */ * nofNodes * sizeof(float);
@@ -527,7 +640,7 @@ string WbWriterVtkXmlBinary::writeTrianglesWithNodeData(const string &filename,
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -561,7 +674,7 @@ string WbWriterVtkXmlBinary::writeTrianglesWithNodeData(const string &filename,
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -625,7 +738,7 @@ string WbWriterVtkXmlBinary::writeQuads(const string &filename, vector<UbTupleFl
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 4 /*nodes per quad  */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 4 /*nodes per quad  */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per quad */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of quad    */ * nofCells * sizeof(unsigned char);
 
@@ -649,7 +762,7 @@ string WbWriterVtkXmlBinary::writeQuads(const string &filename, vector<UbTupleFl
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -674,7 +787,7 @@ string WbWriterVtkXmlBinary::writeQuads(const string &filename, vector<UbTupleFl
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -730,7 +843,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeData(const string &filename, vect
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 4 /*nodes per quad  */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 4 /*nodes per quad  */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per quad */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of quad    */ * nofCells * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar          */ * nofNodes * sizeof(float);
@@ -755,7 +868,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeData(const string &filename, vect
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -789,7 +902,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeData(const string &filename, vect
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -855,7 +968,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithCellData(const string &filename, vect
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 4 /*nodes per quad  */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 4 /*nodes per quad  */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per quad */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of quad    */ * nofCells * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar          */ * nofCells * sizeof(float);
@@ -880,7 +993,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithCellData(const string &filename, vect
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -914,7 +1027,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithCellData(const string &filename, vect
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -984,7 +1097,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeAndCellData(const string &filenam
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 4 /*nodes per quad  */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 4 /*nodes per quad  */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per quad */ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of quad    */ * nofCells * sizeof(unsigned char);
     int bytesScalarDataPoint = 1 /*scalar          */ * nofNodes * sizeof(float);
@@ -1010,7 +1123,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeAndCellData(const string &filenam
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1052,7 +1165,7 @@ string WbWriterVtkXmlBinary::writeQuadsWithNodeAndCellData(const string &filenam
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -1128,7 +1241,7 @@ string WbWriterVtkXmlBinary::writeOctsWithCellData(const string &filename, vecto
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3      */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 8 /*nodes per oct */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 8 /*nodes per oct */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per oct*/ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of oct   */ * nofCells * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar        */ * nofCells * sizeof(float);
@@ -1153,7 +1266,7 @@ string WbWriterVtkXmlBinary::writeOctsWithCellData(const string &filename, vecto
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1187,7 +1300,7 @@ string WbWriterVtkXmlBinary::writeOctsWithCellData(const string &filename, vecto
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -1257,7 +1370,7 @@ string WbWriterVtkXmlBinary::writeOctsWithNodeData(const string &filename, vecto
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3      */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 8 /*nodes per oct */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 8 /*nodes per oct */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per oct*/ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of oct   */ * nofCells * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar        */ * nofNodes * sizeof(double);
@@ -1282,7 +1395,7 @@ string WbWriterVtkXmlBinary::writeOctsWithNodeData(const string &filename, vecto
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1316,7 +1429,7 @@ string WbWriterVtkXmlBinary::writeOctsWithNodeData(const string &filename, vecto
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -1386,7 +1499,7 @@ string WbWriterVtkXmlBinary::writeOcts(const string &filename, vector<UbTupleFlo
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3      */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 8 /*nodes per oct */ * nofCells * sizeof(int);
+    int bytesCellConnectivity = 8 /*nodes per oct */ * nofCells * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per oct*/ * nofCells * sizeof(int);
     int bytesCellTypes       = 1 /*type of oct   */ * nofCells * sizeof(unsigned char);
     // int bytesScalarData      = 1 /*scalar        */ * nofNodes * sizeof(float);
@@ -1411,7 +1524,7 @@ string WbWriterVtkXmlBinary::writeOcts(const string &filename, vector<UbTupleFlo
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1436,7 +1549,7 @@ string WbWriterVtkXmlBinary::writeOcts(const string &filename, vector<UbTupleFlo
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofCells; c++) {
         out.write((char *)&val<1>(cells[c]), sizeof(int));
         out.write((char *)&val<2>(cells[c]), sizeof(int));
@@ -1491,7 +1604,7 @@ std::string WbWriterVtkXmlBinary::writeNodes(const std::string &filename, std::v
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3        */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 1 /*nodes per cell  */ * nofNodes * sizeof(int);
+    int bytesCellConnectivity = 1 /*nodes per cell  */ * nofNodes * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per cell */ * nofNodes * sizeof(int);
     int bytesCellTypes       = 1 /*type of line    */ * nofNodes * sizeof(unsigned char);
 
@@ -1515,7 +1628,7 @@ std::string WbWriterVtkXmlBinary::writeNodes(const std::string &filename, std::v
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1540,7 +1653,7 @@ std::string WbWriterVtkXmlBinary::writeNodes(const std::string &filename, std::v
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofNodes; c++)
         out.write((char *)&c, sizeof(int));
 
@@ -1586,7 +1699,7 @@ std::string WbWriterVtkXmlBinary::writeNodesWithNodeData(const std::string &file
 
     int bytesPerByteVal      = 4; //==sizeof(int)
     int bytesPoints          = 3 /*x1/x2/x3       */ * nofNodes * sizeof(float);
-    int bytesCellConnectivty = 1 /*nodes per cell */ * nofNodes * sizeof(int);
+    int bytesCellConnectivity = 1 /*nodes per cell */ * nofNodes * sizeof(int);
     int bytesCellOffsets     = 1 /*offset per cell*/ * nofNodes * sizeof(int);
     int bytesCellTypes       = 1 /*type of oct    */ * nofNodes * sizeof(unsigned char);
     int bytesScalarData      = 1 /*scalar         */ * nofNodes * sizeof(double);
@@ -1611,7 +1724,7 @@ std::string WbWriterVtkXmlBinary::writeNodesWithNodeData(const std::string &file
     out << "         <Cells>\n";
     out << "            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
-    offset += (bytesPerByteVal + bytesCellConnectivty);
+    offset += (bytesPerByteVal + bytesCellConnectivity);
     out << "            <DataArray type=\"Int32\" Name=\"offsets\" format=\"appended\" offset=\"" << offset
         << "\" />\n";
     offset += (bytesPerByteVal + bytesCellOffsets);
@@ -1645,7 +1758,7 @@ std::string WbWriterVtkXmlBinary::writeNodesWithNodeData(const std::string &file
 
     // CELLS SECTION
     // cellConnectivity
-    out.write((char *)&bytesCellConnectivty, bytesPerByteVal);
+    out.write((char *)&bytesCellConnectivity, bytesPerByteVal);
     for (int c = 0; c < nofNodes; c++)
         out.write((char *)&c, sizeof(int));
 
diff --git a/src/basics/basics/writer/WbWriterVtkXmlBinary.h b/src/basics/basics/writer/WbWriterVtkXmlBinary.h
index 421148d90497e3628ed274439c0b2fd7636b7fd2..0f2c31eda81ad0c1975c9715ac1b7fb37a06339b 100644
--- a/src/basics/basics/writer/WbWriterVtkXmlBinary.h
+++ b/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
diff --git a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.cpp b/src/gpu/GridGenerator/grid/BoundaryConditions/Side.cpp
index 5b191ee4e3fcdc0ec71633111085f70c5dc43479..ba4eea50ffb6bc136528db31207274d626fe9b15 100644
--- a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.cpp
+++ b/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)
 {
-
     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] = 0.5;
+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)
@@ -201,7 +262,7 @@ uint Side::getIndex(SPtr<Grid> grid, std::string coord, real constant, real v1,
 }
 
 
-void Geometry::addIndices(std::vector<SPtr<Grid> > grids, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void Geometry::addIndices(const std::vector<SPtr<Grid>> &grids, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     auto geometryBoundaryCondition = std::dynamic_pointer_cast<GeometryBoundaryCondition>(boundaryCondition);
 
@@ -242,7 +303,7 @@ void Geometry::addIndices(std::vector<SPtr<Grid> > grids, uint level, SPtr<Bound
 
 
 
-void MX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void MX::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartY();
     real endInner = grid[level]->getEndY();
@@ -258,7 +319,7 @@ void MX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCond
 
 }
 
-void PX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void PX::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartY();
     real endInner = grid[level]->getEndY();
@@ -273,7 +334,7 @@ void PX::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCond
     Side::addIndices(grid[level], boundaryCondition, "x", coordinateNormal, startInner, endInner, startOuter, endOuter);
 }
 
-void MY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void MY::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartX();
     real endInner = grid[level]->getEndX();
@@ -289,7 +350,7 @@ void MY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCond
 }
 
 
-void PY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void PY::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartX();
     real endInner = grid[level]->getEndX();
@@ -305,7 +366,7 @@ void PY::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCond
 }
 
 
-void MZ::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void MZ::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartX();
     real endInner = grid[level]->getEndX();
@@ -320,7 +381,7 @@ void MZ::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCond
     Side::addIndices(grid[level], boundaryCondition, "z", coordinateNormal, startInner, endInner, startOuter, endOuter);
 }
 
-void PZ::addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
+void PZ::addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<BoundaryCondition> boundaryCondition)
 {
     real startInner = grid[level]->getStartX();
     real endInner = grid[level]->getEndX();
diff --git a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h b/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h
index 53a763bc562ee978042b28d24856fbcca256c5f9..624b3722a1c909ba26063b49565779b924d34adc 100644
--- a/src/gpu/GridGenerator/grid/BoundaryConditions/Side.h
+++ b/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,20 +63,19 @@ enum class SideType
     MX, PX, MY, PY, MZ, PZ, GEOMETRY
 };
 
-
-
 class Side
 {
 public:
     virtual ~Side() = default;
-    virtual void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) = 0;
+    virtual void addIndices(const std::vector<SPtr<Grid>> &grid, uint level,
+                            SPtr<gg::BoundaryCondition> boundaryCondition) = 0;
 
     virtual int getCoordinate() const = 0;
     virtual int getDirection() const = 0;
 
     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,14 +87,35 @@ protected:
 
     void setQs(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint index);
 
+    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
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -112,7 +136,7 @@ public:
 class MX : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -133,7 +157,7 @@ public:
 class PX : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -155,7 +179,7 @@ public:
 class MY : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -176,7 +200,7 @@ public:
 class PY : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -198,7 +222,7 @@ public:
 class MZ : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
@@ -219,7 +243,7 @@ public:
 class PZ : public Side
 {
 public:
-    void addIndices(std::vector<SPtr<Grid> > grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override;
 
     int getCoordinate() const override
     {
diff --git a/src/gpu/GridGenerator/grid/BoundaryConditions/SideTest.cpp b/src/gpu/GridGenerator/grid/BoundaryConditions/SideTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..36a286a8766db4af7e109eb3f8d47add401779f9
--- /dev/null
+++ b/src/gpu/GridGenerator/grid/BoundaryConditions/SideTest.cpp
@@ -0,0 +1,873 @@
+#include "Side.h"
+#include "PointerDefinitions.h"
+#include "gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h"
+#include "grid/GridImp.h"
+#include "grid/NodeValues.h"
+#include "lbm/constants/D3Q27.h"
+#include "gmock/gmock.h"
+#include <algorithm>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+using namespace vf::gpu;
+using namespace vf::lbm::dir;
+
+class SideTestSpecificSubclass : public Side
+{
+
+public:
+    void setQs(SPtr<Grid> grid, SPtr<gg::BoundaryCondition> boundaryCondition, uint index)
+    {
+        Side::setQs(grid, boundaryCondition, index);
+    };
+    int sideDirection = POSITIVE_DIR;
+    int coordinateDirection = X_INDEX;
+    SideType mySide = SideType::PX;
+
+private:
+    void correctNeighborForPeriodicBoundaries(const Grid *grid, std::array<real, 3>& coords, std::array<real, 3>& neighbors) const override
+    {
+    }
+
+    int getDirection() const override
+    {
+        return sideDirection;
+    }
+
+    void addIndices(const std::vector<SPtr<Grid>> &grid, uint level, SPtr<gg::BoundaryCondition> boundaryCondition) override
+    {
+    }
+
+    int getCoordinate() const override
+    {
+        return coordinateDirection;
+    }
+
+    SideType whoAmI() const override
+    {
+        return mySide;
+    }
+};
+
+class GridDouble : public GridImp
+{
+
+public:
+    int endDirection = -1;
+
+    GridDouble()
+    {
+        this->distribution = DistributionHelper::getDistribution27();
+    }
+
+    void transIndexToCoords(uint index, real &x, real &y, real &z) const override
+    {
+        x = 0;
+        y = 0;
+        z = 0;
+    }
+
+    real getDelta() const override
+    {
+        return 1.0;
+    }
+
+    uint transCoordToIndex(const real &x, const real &y, const real &z) const override
+    {
+        return 0;
+    }
+
+    char getFieldEntry(uint /*matrixIndex*/) const override
+    {
+        return STOPPER_OUT_OF_GRID_BOUNDARY;
+    }
+
+    int getEndDirection() const override
+    {
+        return endDirection;
+    }
+};
+
+class BoundaryConditionSpy : public gg::BoundaryCondition
+{
+public:
+    char getType() const override
+    {
+        return 't';
+    };
+    const std::vector<std::vector<real>> &getQs()
+    {
+        return this->qs;
+    }
+    void resetQVector()
+    {
+        this->qs.clear();
+    }
+};
+
+class SideTestBC : public testing::Test
+{
+protected:
+    SideTestSpecificSubclass side;
+    SPtr<GridDouble> grid = std::make_shared<GridDouble>();
+    SPtr<BoundaryConditionSpy> bc = std::make_shared<BoundaryConditionSpy>();
+    uint index = 0;
+
+    std::vector<real> noBC;
+
+    void SetUp() override
+    {
+        grid->endDirection = 26;
+    }
+};
+
+TEST_F(SideTestBC, setQs2D_whenSettingPX_setAllQsNormalToBC)
+{
+    grid->endDirection = 10;
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(11, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs2D_givenPYhasBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    grid->endDirection = 10;
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(11, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMXhasBeenSet_thenSetPX_setAllQsNormalToPX)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+
+    // no previous BC on this node
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+
+    // node already has BC in MX direction, but this does not change anything
+
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    actualQs = bc->getQs()[0];
+
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenGeometryBCInVector_thenSetPX_throws)
+{
+    // do not add Geometry BC to this vector, as it has an invalid normal
+    grid->addBCalreadySet(SideType::GEOMETRY);
+
+    EXPECT_THROW(side.setQs(grid, bc, index), std::out_of_range);
+}
+
+TEST_F(SideTestBC, setQs3D_whenSettingPX_setAllQsNormalToBC)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYhasBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYhasBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPZhasBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMZhasBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_P00] = 0.5;
+    expectedQs[DIR_PP0] = 0.5;
+    expectedQs[DIR_PM0] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandMZhaveBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_P00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PM0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandPZhaveBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_P00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PM0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandPZhaveBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_P00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PP0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandMZhaveBeenSet_thenSetPX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_P00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PP0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PPP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_whenSettingMX_setAllQsNormalToBC)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_M00] = 0.5;
+    expectedQs[DIR_MP0] = 0.5;
+    expectedQs[DIR_MM0] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYhasBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_M00] = 0.5;
+    expectedQs[DIR_MM0] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYhasBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_M00] = 0.5;
+    expectedQs[DIR_MP0] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPZhasBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_M00] = 0.5;
+    expectedQs[DIR_MP0] = 0.5;
+    expectedQs[DIR_MM0] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMZhasBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_M00] = 0.5;
+    expectedQs[DIR_MP0] = 0.5;
+    expectedQs[DIR_MM0] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandMZhaveBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_M00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MM0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandPZhaveBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_M00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MM0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandPZhaveBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::PZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_M00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MP0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandMZhaveBeenSet_thenSetMX_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = X_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::MZ);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_M00] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MP0] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MPP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_whenSettingMZ_setAllQsNormalToBC)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00M] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_0PM] = 0.5;
+    expectedQs[DIR_0MM] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYhasBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00M] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_0PM] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYhasBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00M] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_0MM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPXhasBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00M] = 0.5;
+    expectedQs[DIR_M0M] = 0.5;
+    expectedQs[DIR_0PM] = 0.5;
+    expectedQs[DIR_0MM] = 0.5;
+    expectedQs[DIR_MPM] = 0.5;
+    expectedQs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMXhasBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00M] = 0.5;
+    expectedQs[DIR_P0M] = 0.5;
+    expectedQs[DIR_0PM] = 0.5;
+    expectedQs[DIR_0MM] = 0.5;
+    expectedQs[DIR_PPM] = 0.5;
+    expectedQs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandPXhaveBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0PM] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandMXhaveBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0PM] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PPM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandPXhaveBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0MM] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandMXhaveBeenSet_thenSetMZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = NEGATIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0M] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0MM] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PMM] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_whenSettingPZ_setAllQsNormalToBC)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00P] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_0PP] = 0.5;
+    expectedQs[DIR_0MP] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYhasBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00P] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_0PP] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYhasBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00P] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_0MP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPXhasBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00P] = 0.5;
+    expectedQs[DIR_M0P] = 0.5;
+    expectedQs[DIR_0PP] = 0.5;
+    expectedQs[DIR_0MP] = 0.5;
+    expectedQs[DIR_MPP] = 0.5;
+    expectedQs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMXhasBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQs(27, -1);
+    expectedQs[DIR_00P] = 0.5;
+    expectedQs[DIR_P0P] = 0.5;
+    expectedQs[DIR_0PP] = 0.5;
+    expectedQs[DIR_0MP] = 0.5;
+    expectedQs[DIR_PPP] = 0.5;
+    expectedQs[DIR_PMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandPXhaveBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0PP] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MPP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenMYandMXhaveBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::MY);
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0PP] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PPP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandPXhaveBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::PX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_M0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0MP] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_MMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
+
+TEST_F(SideTestBC, setQs3D_givenPYandMXhaveBeenSet_thenSetPZ_doNotSetSameQsAgain)
+{
+    side.coordinateDirection = Z_INDEX;
+    side.sideDirection = POSITIVE_DIR;
+    grid->addBCalreadySet(SideType::PY);
+    grid->addBCalreadySet(SideType::MX);
+
+    side.setQs(grid, bc, index);
+    auto actualQs = bc->getQs()[0];
+
+    std::vector<real> expectedQsForTwoPreviousBCs(27, -1);
+    expectedQsForTwoPreviousBCs[DIR_00P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_P0P] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_0MP] = 0.5;
+    expectedQsForTwoPreviousBCs[DIR_PMP] = 0.5;
+    EXPECT_THAT(actualQs, testing::Eq(expectedQsForTwoPreviousBCs));
+}
diff --git a/src/gpu/GridGenerator/grid/Grid.h b/src/gpu/GridGenerator/grid/Grid.h
index 85b19bedd470c9856954a1ca20fb446c3d875da2..ad2ce473fb65fe4414f6da5c4caf0d3e140b7e02 100644
--- a/src/gpu/GridGenerator/grid/Grid.h
+++ b/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
diff --git a/src/gpu/GridGenerator/grid/GridImp.cpp b/src/gpu/GridGenerator/grid/GridImp.cpp
index 05c684410166e329ba63bbe3bdbf0c09e3a881ab..32cf9d07da87149695a5bf548ed357be2b2f71b4 100644
--- a/src/gpu/GridGenerator/grid/GridImp.cpp
+++ b/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;
 }
@@ -2188,8 +2194,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
@@ -2199,3 +2211,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);
+}
diff --git a/src/gpu/GridGenerator/grid/GridImp.h b/src/gpu/GridGenerator/grid/GridImp.h
index 8283bf569e266b84f020334a306d93756b01c394..2cd322ebed78daaf135ad97b881923ca5831bbcd 100644
--- a/src/gpu/GridGenerator/grid/GridImp.h
+++ b/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();
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
index 6f885190d25eeccf7bf4792f46be0f84880f7947..38a7eef7e356e2f2da4c1a819d8375035a37313a 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
+++ b/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];
 }
diff --git a/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d567c695a0e33b7a88c2c8cf3bcb88093ce5b802
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriter.hpp
@@ -0,0 +1,96 @@
+#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, const 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
diff --git a/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriterTest.cpp b/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriterTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9eecb25c663fcfc8fde353b76ccf20cbcb9cf272
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/Output/QDebugVtkWriterTest.cpp
@@ -0,0 +1,60 @@
+#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));
+}
diff --git a/src/gpu/VirtualFluids_GPU/Output/QDebugWriter.hpp b/src/gpu/VirtualFluids_GPU/Output/QDebugWriter.hpp
index d006636572377477aeb3599a8ae843ea2b1e31ff..c1a3658d318eb47e84530bf437afa0bb6ba91743 100644
--- a/src/gpu/VirtualFluids_GPU/Output/QDebugWriter.hpp
+++ b/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)