diff --git a/apps/cpu/FlowAroundCylinder/cylinder.cfg b/apps/cpu/FlowAroundCylinder/cylinder.cfg
index 0a7066ed9bc3351736c511d7aaeecaa04604fe55..97ece40e65d4ffe47a75e5377db49bd0018bbff6 100644
--- a/apps/cpu/FlowAroundCylinder/cylinder.cfg
+++ b/apps/cpu/FlowAroundCylinder/cylinder.cfg
@@ -1,6 +1,6 @@
-pathOut = d:/temp/cylinder_test
+pathOut = d:/temp/cylinder_test_naming
 
-numOfThreads = 4
+numOfThreads = 8
 availMem = 15e9
 refineLevel = 0
 blockNx = 25 41 41
@@ -16,7 +16,7 @@ restartStep = 1000
 cpStart = 1000
 cpStep = 1000
 
-outTime = 10000
-endTime = 100000
+outTime = 10
+endTime = 100
 
 nupsStep = 100 100 10000000
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
index b7e4d4b9a6cc3276728b5d07d89ac3c723ab027c..ea45bb14110a071724f816b3c7840ce0dfbd7327 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
@@ -9,7 +9,7 @@
 #include "D3Q27System.h"
 #include "DataSet3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 #include "UbFileInputASCII.h"
 #include "UbFileOutputASCII.h"
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
index fc993d95808036d1645e7b1f43fff40336f63641..c0d88a853b0a739ee1fc0529273d9dec26fd2232 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -228,14 +228,14 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                         ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
                         (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
                         ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
-                            ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+                            ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000];
                     if (distributionsH2) {
                     distributionsH2->getDistribution(f, ix1, ix2, ix3);
                     (*phaseField2)(ix1, ix2, ix3) =
                         ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
                         (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
                         ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
-                            ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+                            ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000];
                 }
                     else { (*phaseField2)(ix1, ix2, ix3) = 999.0; }
                     
@@ -273,8 +273,8 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     nodes.push_back(UbTupleFloat3(float(worldCoordinates[0]), float(worldCoordinates[1]),
                                                   float(worldCoordinates[2])));
 
-                    phi[REST] = (*phaseField)(ix1, ix2, ix3);
-                    phi2[REST] = (*phaseField2)(ix1, ix2, ix3);
+                    phi[DIR_000] = (*phaseField)(ix1, ix2, ix3);
+                    phi2[DIR_000] = (*phaseField2)(ix1, ix2, ix3);
 
                     if ((ix1 == 0) || (ix2 == 0) || (ix3 == 0)) {
                         dX1_phi = 0.0;
@@ -314,7 +314,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                         dX1_phi  = 0.0 * gradX1_phi(phi);
                         dX2_phi  = 0.0 * gradX2_phi(phi);
                         dX3_phi  = 0.0 * gradX3_phi(phi);
-                        mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
+                        mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi(phi);
 
                         //phi2[E] = (*phaseField2)(ix1 + DX1[E], ix2 + DX2[E], ix3 + DX3[E]);
                         //phi2[N] = (*phaseField2)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
@@ -358,7 +358,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     LBMReal rhoToPhi = (rhoH - rhoL) / (phiH - phiL);
 
                     // rho = phi[ZERO] + (1.0 - phi[ZERO])*1.0/densityRatio;
-                    rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                    rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
                     if (pressure) {
                         vx1 =
@@ -398,7 +398,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     p1 = (((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
                           (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
                            ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
-                          ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST]) +
+                          ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000]) +
                          (vx1 * rhoToPhi * dX1_phi * c1o3 + vx2 * rhoToPhi * dX2_phi * c1o3 +
                           vx3 * rhoToPhi * dX3_phi * c1o3) /
                              2.0;
@@ -426,7 +426,7 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
                                            UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
 
-                    if (UbMath::isNaN(phi[REST]) || UbMath::isInfinity(phi[REST]))
+                    if (UbMath::isNaN(phi[DIR_000]) || UbMath::isInfinity(phi[DIR_000]))
                         UB_THROW(UbException(
                             UB_EXARGS, "phi is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
@@ -436,12 +436,12 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                         UB_THROW( UbException(UB_EXARGS,"p1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
                         ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
 
-                    data[index++].push_back(phi[REST]);
+                    data[index++].push_back(phi[DIR_000]);
                     data[index++].push_back(vx1);
                     data[index++].push_back(vx2);
                     data[index++].push_back(vx3);
                     data[index++].push_back(p1);
-                    data[index++].push_back(phi2[REST]);
+                    data[index++].push_back(phi2[DIR_000]);
                     if (pressure) data[index++].push_back((*pressure)(ix1, ix2, ix3));
                 }
             }
@@ -502,7 +502,7 @@ LBMReal WriteMultiphaseQuantitiesCoProcessor::nabla2_phi(const LBMReal *const &h
     using namespace D3Q27System;
     LBMReal sum = 0.0;
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
-        sum += WEIGTH[k] * (h[k] - h[REST]);
+        sum += WEIGTH[k] * (h[k] - h[DIR_000]);
     }
     return 6.0 * sum;
 }
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
index 739efcddb9ceea5c0951df83833d64ad90bb02c5..21dbee2f9e0f931d1a6400d03784c1612e61e3d2 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
@@ -18,7 +18,7 @@ void OneDistributionFullVectorConnector::init()
 
     int anz = 27;
     switch (sendDir) {
-        case D3Q27System::REST:
+        case D3Q27System::DIR_000:
             UB_THROW(UbException(UB_EXARGS, "ZERO not allowed"));
             break;
         case D3Q27System::E:
diff --git a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
index 3e314fbce45e8d29affa09777a7cbeb1b92ca2ac..7d2e5e00a798fb3fecca426f601254d9f7d0729e 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
@@ -60,7 +60,7 @@ void ThreeDistributionsDoubleGhostLayerFullVectorConnector::init()
    int anz = 3*27+1;
    switch (sendDir)
    {
-   case D3Q27System::REST: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
+   case D3Q27System::DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
    case D3Q27System::E:
    case D3Q27System::W: sender->getData().resize(maxX2*maxX3*anz*2, 0.0);   break;
    case D3Q27System::N:
diff --git a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
index 2e726fc7b88c9ef229e503924eadcf53a9b06dfd..5cd49a163a559eba2ef74c498f54befeeb5e696f 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
@@ -59,7 +59,7 @@ void ThreeDistributionsFullVectorConnector::init()
    int anz = 3*27;
    switch (sendDir)
    {
-   case D3Q27System::REST: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
+   case D3Q27System::DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
    case D3Q27System::E:
    case D3Q27System::W: sender->getData().resize(maxX2*maxX3*anz, 0.0);   break;
    case D3Q27System::N:
diff --git a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
index b831cb840b3de012c1df35efd237e3e0c3f81056..15d18279c315257df777b009534f093f667db064 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
@@ -59,7 +59,7 @@ void TwoDistributionsDoubleGhostLayerFullVectorConnector::init()
    int anz = 2*27+1;
    switch (sendDir)
    {
-   case D3Q27System::REST: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
+   case D3Q27System::DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
    case D3Q27System::E:
    case D3Q27System::W: sender->getData().resize(maxX2*maxX3*anz*2, 0.0);   break;
    case D3Q27System::N:
diff --git a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
index 7fe8bc3643c337323ef25ee35c260597744e6191..008b6630d41fa8b0c9f93d99ded1e044ec455b73 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
@@ -58,7 +58,7 @@ void TwoDistributionsFullVectorConnector::init()
    int anz = 2*27;
    switch (sendDir)
    {
-   case D3Q27System::REST: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
+   case D3Q27System::DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
    case D3Q27System::E:
    case D3Q27System::W: sender->getData().resize(maxX2*maxX3*anz, 0.0);   break;
    case D3Q27System::N:
diff --git a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSoA.cpp b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSoA.cpp
index bd6d46c2bdaeb72244578b4e3f3625cd2dfe7ff1..765461b34b40151ef8f7364d916d2d9a1b885c11 100644
--- a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSoA.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSoA.cpp
@@ -119,7 +119,7 @@ void D3Q27EsoTwist3DSoA::getDistribution(LBMReal *const f, size_t x1, size_t x2,
     f[D3Q27System::BNW] = (*d.BNW)(x1p, x2, x3p);
     f[D3Q27System::BNE] = (*d.BNE)(x1, x2, x3p);
 
-    f[D3Q27System::REST] = (*d.REST)(x1, x2, x3);
+    f[D3Q27System::DIR_000] = (*d.REST)(x1, x2, x3);
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSoA::setDistribution(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -156,7 +156,7 @@ void D3Q27EsoTwist3DSoA::setDistribution(const LBMReal *const f, size_t x1, size
     (*d.BNW)(x1p, x2, x3p)  = f[D3Q27System::INV_BNW];
     (*d.BNE)(x1, x2, x3p)   = f[D3Q27System::INV_BNE];
 
-    (*d.REST)(x1, x2, x3) = f[D3Q27System::REST];
+    (*d.REST)(x1, x2, x3) = f[D3Q27System::DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSoA::getDistributionInv(LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -189,7 +189,7 @@ void D3Q27EsoTwist3DSoA::getDistributionInv(LBMReal *const f, size_t x1, size_t
     f[D3Q27System::INV_BNW] = (*d.BNW)(x1 + 1, x2, x3 + 1);
     f[D3Q27System::INV_BNE] = (*d.BNE)(x1, x2, x3 + 1);
 
-    f[D3Q27System::REST] = (*d.REST)(x1, x2, x3);
+    f[D3Q27System::DIR_000] = (*d.REST)(x1, x2, x3);
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSoA::setDistributionInv(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
diff --git a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
index 5e762c68bab806ee7c892c000869bce8c76431af..ecb53fda11681ada356db6af2ad40996aebef422 100644
--- a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
@@ -84,7 +84,7 @@ void D3Q27EsoTwist3DSplittedVector::getDistribution(LBMReal *const f, size_t x1,
     f[D3Q27System::BNW] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
     f[D3Q27System::BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
 
-    f[D3Q27System::REST] = (*this->zeroDistributions)(x1, x2, x3);
+    f[D3Q27System::DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::setDistribution(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -117,7 +117,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistribution(const LBMReal *const f, size
     (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1)     = f[D3Q27System::INV_BNW];
     (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1)         = f[D3Q27System::INV_BNE];
 
-    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::getDistributionInv(LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -150,7 +150,7 @@ void D3Q27EsoTwist3DSplittedVector::getDistributionInv(LBMReal *const f, size_t
     f[D3Q27System::INV_BNW] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
     f[D3Q27System::INV_BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
 
-    f[D3Q27System::REST] = (*this->zeroDistributions)(x1, x2, x3);
+    f[D3Q27System::DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::setDistributionInv(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -183,7 +183,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInv(const LBMReal *const f, s
     (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1)     = f[D3Q27System::BNW];
     (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1)         = f[D3Q27System::BNE];
 
-    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(const LBMReal *const f, size_t x1, size_t x2, size_t x3,
@@ -242,7 +242,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(const LBMReal *c
     if ((direction & EsoTwistD3Q27System::etTSW) == EsoTwistD3Q27System::etTSW)
         (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[D3Q27System::TSW];
     if ((direction & EsoTwistD3Q27System::REST) == EsoTwistD3Q27System::REST)
-        (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+        (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(LBMReal f, size_t x1, size_t x2, size_t x3,
@@ -327,7 +327,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(LBMReal f, size_
         case D3Q27System::TSW:
             (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f;
             break;
-        case D3Q27System::REST:
+        case D3Q27System::DIR_000:
             (*this->zeroDistributions)(x1, x2, x3) = f;
             break;
         default:
@@ -391,7 +391,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(const LBMReal
     if ((direction & EsoTwistD3Q27System::etTSW) == EsoTwistD3Q27System::etTSW)
         (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f[D3Q27System::TSW];
     if ((direction & EsoTwistD3Q27System::REST) == EsoTwistD3Q27System::REST)
-        (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+        (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(LBMReal f, size_t x1, size_t x2, size_t x3,
@@ -476,7 +476,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(LBMReal f, si
         case D3Q27System::TSW:
             (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f;
             break;
-        case D3Q27System::REST:
+        case D3Q27System::DIR_000:
             (*this->zeroDistributions)(x1, x2, x3) = f;
             break;
         default:
@@ -539,7 +539,7 @@ LBMReal D3Q27EsoTwist3DSplittedVector::getDistributionForDirection(size_t x1, si
             return (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
         case D3Q27System::BNE:
             return (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
-        case D3Q27System::REST:
+        case D3Q27System::DIR_000:
             return (*this->zeroDistributions)(x1, x2, x3);
         default:
             UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
@@ -601,7 +601,7 @@ LBMReal D3Q27EsoTwist3DSplittedVector::getDistributionInvForDirection(size_t x1,
             return (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
         case D3Q27System::TSW:
             return (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
-        case D3Q27System::REST:
+        case D3Q27System::DIR_000:
             return (*this->zeroDistributions)(x1, x2, x3);
         default:
             UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
diff --git a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
index c456be678449744475a0ac6932850dceb0ee6f1c..21d1141a126993876242d55bed562bbf680597a1 100644
--- a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
@@ -49,7 +49,7 @@ const int EsoTwistD3Q27System::etINVDIR[EsoTwistD3Q27System::ENDF + 1] = {
     D3Q27System::INV_TE,  D3Q27System::INV_BW,  D3Q27System::INV_BE,  D3Q27System::INV_TW,  D3Q27System::INV_TN,
     D3Q27System::INV_BS,  D3Q27System::INV_BN,  D3Q27System::INV_TS,  D3Q27System::INV_TNE, D3Q27System::INV_TNW,
     D3Q27System::INV_TSE, D3Q27System::INV_TSW, D3Q27System::INV_BNE, D3Q27System::INV_BNW, D3Q27System::INV_BSE,
-    D3Q27System::INV_BSW, D3Q27System::REST
+    D3Q27System::INV_BSW, D3Q27System::DIR_000
 };
 
 const unsigned long int EsoTwistD3Q27System::etDIR[EsoTwistD3Q27System::ENDF + 1] = {
diff --git a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
index 21752cc48a84b02bc24cb7efe9e3c5912f476dfd..c2ea5a9ddaba852ec5aaf28708b0cc02f5341913 100644
--- a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
@@ -47,7 +47,7 @@ struct EsoTwistD3Q27System {
     const static int STARTDIR = D3Q27System::STARTDIR;
     const static int ENDDIR   = D3Q27System::ENDDIR;
 
-    static const int REST = D3Q27System::REST; /*f0 */
+    static const int REST = D3Q27System::DIR_000; /*f0 */
     static const int E    = D3Q27System::E;    /*f1 */
     static const int W    = D3Q27System::W;    /*f2 */
     static const int N    = D3Q27System::N;    /*f3 */
diff --git a/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp b/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
index 79753c144f5cfff831f1d0415e9434c50b11bcea..3a55c6ad7d4faecfd8d058f8d044c7da1247f962 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
@@ -34,7 +34,7 @@
 #include "Block3D.h"
 
 #include "Block3DConnector.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 
 int Block3D::counter = 0;
@@ -335,13 +335,13 @@ std::string Block3D::toString()
     for (std::size_t i = 0; i < connectors.size(); i++)
         if (connectors[i]) {
             if (connectors[i]->isLocalConnector())
-                ss << "l." << Grid3DSystem::getDirectionString(connectors[i]->getSendDir()) << ", ";
+                ss << "l." << D3Q27System::getDirectionString(connectors[i]->getSendDir()) << ", ";
             if (connectors[i]->isRemoteConnector())
-                ss << "r." << Grid3DSystem::getDirectionString(connectors[i]->getSendDir()) << ", ";
+                ss << "r." << D3Q27System::getDirectionString(connectors[i]->getSendDir()) << ", ";
             if (connectors[i]->isInterpolationConnectorCF())
-                ss << "cf." << Grid3DSystem::getDirectionString(connectors[i]->getSendDir()) << ", ";
+                ss << "cf." << D3Q27System::getDirectionString(connectors[i]->getSendDir()) << ", ";
             if (connectors[i]->isInterpolationConnectorFC())
-                ss << "fc." << Grid3DSystem::getDirectionString(connectors[i]->getSendDir()) << ", ";
+                ss << "fc." << D3Q27System::getDirectionString(connectors[i]->getSendDir()) << ", ";
         }
     return ss.str();
 }
diff --git a/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp b/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
index 41235087e4d48f81c1581e39dd8c0a3a54a7779a..b4d3ba655677893d0139c15ff6d25385f1e42f4d 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
@@ -39,21 +39,21 @@
 #include <geometry3d/CoordinateTransformation3D.h>
 
 #include "Block3DVisitor.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "Grid3DVisitor.h"
 #include "Interactor3D.h"
-#include "LBMSystem.h"
+#include "D3Q27System.h"
 #include <Block3D.h>
 #include <Communicator.h>
 
 using namespace std;
 
-Grid3D::Grid3D() { levelSet.resize(Grid3DSystem::MAXLEVEL + 1); }
+Grid3D::Grid3D() { levelSet.resize(D3Q27System::MAXLEVEL + 1); }
 //////////////////////////////////////////////////////////////////////////
 Grid3D::Grid3D(std::shared_ptr<vf::mpi::Communicator> comm)
 
 {
-    levelSet.resize(Grid3DSystem::MAXLEVEL + 1);
+    levelSet.resize(D3Q27System::MAXLEVEL + 1);
     bundle = comm->getBundleID();
     rank = comm->getProcessID();
 }
@@ -63,7 +63,7 @@ Grid3D::Grid3D(std::shared_ptr<vf::mpi::Communicator> comm, int blockNx1, int bl
 
       blockNx1(blockNx1), blockNx2(blockNx2), blockNx3(blockNx2), nx1(gridNx1), nx2(gridNx2), nx3(gridNx3)
 {
-    levelSet.resize(Grid3DSystem::MAXLEVEL + 1);
+    levelSet.resize(D3Q27System::MAXLEVEL + 1);
     bundle = comm->getBundleID();
     rank  = comm->getProcessID();
     trafo = std::make_shared<CoordinateTransformation3D>(0.0, 0.0, 0.0, (double)blockNx1, (double)blockNx2,
@@ -88,7 +88,7 @@ void Grid3D::accept(Block3DVisitor &blockVisitor)
     int startLevel = blockVisitor.getStartLevel();
     int stopLevel  = blockVisitor.getStopLevel();
 
-    if (startLevel < 0 || stopLevel < 0 || startLevel > Grid3DSystem::MAXLEVEL || stopLevel > Grid3DSystem::MAXLEVEL)
+    if (startLevel < 0 || stopLevel < 0 || startLevel > D3Q27System::MAXLEVEL || stopLevel > D3Q27System::MAXLEVEL)
         throw UbException(UB_EXARGS, "not valid level!");
 
     bool dir = startLevel < stopLevel;
@@ -158,8 +158,8 @@ bool Grid3D::deleteBlock(int ix1, int ix2, int ix3, int level)
 void Grid3D::deleteBlocks()
 {
     std::vector<std::vector<SPtr<Block3D>>> blocksVector(25);
-    int minInitLevel = Grid3DSystem::MINLEVEL;
-    int maxInitLevel = Grid3DSystem::MAXLEVEL;
+    int minInitLevel = D3Q27System::MINLEVEL;
+    int maxInitLevel = D3Q27System::MAXLEVEL;
     for (int level = minInitLevel; level < maxInitLevel; level++) {
         getBlocks(level, blocksVector[level]);
         for (SPtr<Block3D> block : blocksVector[level]) //	blocks of the current level
@@ -265,7 +265,7 @@ void Grid3D::getSubBlocks(int ix1, int ix2, int ix3, int level, int levelDepth,
         return;
     if (level > 0 && !this->getSuperBlock(ix1, ix2, ix3, level))
         return;
-    if (level >= Grid3DSystem::MAXLEVEL)
+    if (level >= D3Q27System::MAXLEVEL)
         throw UbException(UB_EXARGS, "Level bigger then MAXLEVEL");
 
     int x1[] = { ix1 << 1, (ix1 << 1) + 1 };
@@ -300,7 +300,7 @@ bool Grid3D::expandBlock(int ix1, int ix2, int ix3, int level)
     ix3 = block->getX3();
 
     int l = level + 1;
-    if (l > Grid3DSystem::MAXLEVEL)
+    if (l > D3Q27System::MAXLEVEL)
         throw UbException(UB_EXARGS, "level > Grid3D::MAXLEVEL");
 
     int west   = ix1 << 1;
@@ -584,7 +584,7 @@ void Grid3D::checkLevel(int level)
     if (level < 0) {
         throw UbException(UB_EXARGS, "l(" + UbSystem::toString(level) + (string) ")<0");
     }
-    if (level > Grid3DSystem::MAXLEVEL) {
+    if (level > D3Q27System::MAXLEVEL) {
         throw UbException(UB_EXARGS, "l(" + UbSystem::toString(level) + (string) ")>MAXLEVEL");
     }
     if (this->levelSet[level].size() == 0) {
@@ -596,7 +596,7 @@ bool Grid3D::hasLevel(int level) const
 {
     if (level < 0)
         return false;
-    if (level > Grid3DSystem::MAXLEVEL)
+    if (level > D3Q27System::MAXLEVEL)
         return false;
     if (this->levelSet[level].size() == 0)
         return false;
@@ -616,7 +616,7 @@ UbTupleInt3 Grid3D::getBlockNX() const { return makeUbTuple(blockNx1, blockNx2,
 
 SPtr<Block3D> Grid3D::getNeighborBlock(int dir, int ix1, int ix2, int ix3, int level) const
 {
-    return this->getBlock(ix1 + Grid3DSystem::EX1[dir], ix2 + Grid3DSystem::EX2[dir], ix3 + Grid3DSystem::EX3[dir],
+    return this->getBlock(ix1 + D3Q27System::DX1[dir], ix2 + D3Q27System::DX2[dir], ix3 + D3Q27System::DX3[dir],
                           level);
 }
 //////////////////////////////////////////////////////////////////////////
@@ -631,8 +631,8 @@ SPtr<Block3D> Grid3D::getNeighborBlock(int dir, SPtr<Block3D> block) const
 //////////////////////////////////////////////////////////////////////////
 void Grid3D::getAllNeighbors(int ix1, int ix2, int ix3, int level, int levelDepth, std::vector<SPtr<Block3D>> &blocks)
 {
-    for (int dir = Grid3DSystem::STARTDIR; dir <= Grid3DSystem::ENDDIR; dir++)
-    // for (int dir = Grid3DSystem::STARTDIR; dir<=Grid3DSystem::TS; dir++)
+    for (int dir = D3Q27System::STARTDIR; dir <= D3Q27System::ENDDIR; dir++)
+    // for (int dir = D3Q27System::STARTDIR; dir<=D3Q27System::TS; dir++)
     {
         this->getNeighborBlocksForDirection(dir, ix1, ix2, ix3, level, levelDepth, blocks);
     }
@@ -1100,82 +1100,82 @@ void Grid3D::getNeighborBlocksForDirection(int dir, int ix1, int ix2, int ix3, i
                                            std::vector<SPtr<Block3D>> &blocks)
 {
     switch (dir) {
-        case Grid3DSystem::E:
+        case D3Q27System::E:
             this->getNeighborsEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::W:
+        case D3Q27System::W:
             this->getNeighborsWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::N:
+        case D3Q27System::N:
             this->getNeighborsNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::S:
+        case D3Q27System::S:
             this->getNeighborsSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::T:
+        case D3Q27System::T:
             this->getNeighborsTop(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::B:
+        case D3Q27System::B:
             this->getNeighborsBottom(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::NE:
+        case D3Q27System::NE:
             this->getNeighborsNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::SW:
+        case D3Q27System::SW:
             this->getNeighborsSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::SE:
+        case D3Q27System::SE:
             this->getNeighborsSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::NW:
+        case D3Q27System::NW:
             this->getNeighborsNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TE:
+        case D3Q27System::TE:
             this->getNeighborsTopEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BW:
+        case D3Q27System::BW:
             this->getNeighborsBottomWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BE:
+        case D3Q27System::BE:
             this->getNeighborsBottomEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TW:
+        case D3Q27System::TW:
             this->getNeighborsTopWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TN:
+        case D3Q27System::TN:
             this->getNeighborsTopNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BS:
+        case D3Q27System::BS:
             this->getNeighborsBottomSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BN:
+        case D3Q27System::BN:
             this->getNeighborsBottomNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TS:
+        case D3Q27System::TS:
             this->getNeighborsTopSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TNE:
+        case D3Q27System::TNE:
             this->getNeighborsTopNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TNW:
+        case D3Q27System::TNW:
             this->getNeighborsTopNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TSE:
+        case D3Q27System::TSE:
             this->getNeighborsTopSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TSW:
+        case D3Q27System::TSW:
             this->getNeighborsTopSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BNE:
+        case D3Q27System::BNE:
             this->getNeighborsBottomNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BNW:
+        case D3Q27System::BNW:
             this->getNeighborsBottomNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BSE:
+        case D3Q27System::BSE:
             this->getNeighborsBottomSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BSW:
+        case D3Q27System::BSW:
             this->getNeighborsBottomSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
         default:
@@ -1263,85 +1263,85 @@ void Grid3D::getNeighborBlocksForDirectionWithDirZero(int dir, int ix1, int ix2,
                                                       std::vector<SPtr<Block3D>> &blocks)
 {
     switch (dir) {
-        case Grid3DSystem::E:
+        case D3Q27System::E:
             this->getNeighborsEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::W:
+        case D3Q27System::W:
             this->getNeighborsWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::N:
+        case D3Q27System::N:
             this->getNeighborsNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::S:
+        case D3Q27System::S:
             this->getNeighborsSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::T:
+        case D3Q27System::T:
             this->getNeighborsTop(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::B:
+        case D3Q27System::B:
             this->getNeighborsBottom(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::NE:
+        case D3Q27System::NE:
             this->getNeighborsNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::SW:
+        case D3Q27System::SW:
             this->getNeighborsSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::SE:
+        case D3Q27System::SE:
             this->getNeighborsSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::NW:
+        case D3Q27System::NW:
             this->getNeighborsNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TE:
+        case D3Q27System::TE:
             this->getNeighborsTopEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BW:
+        case D3Q27System::BW:
             this->getNeighborsBottomWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BE:
+        case D3Q27System::BE:
             this->getNeighborsBottomEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TW:
+        case D3Q27System::TW:
             this->getNeighborsTopWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TN:
+        case D3Q27System::TN:
             this->getNeighborsTopNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BS:
+        case D3Q27System::BS:
             this->getNeighborsBottomSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BN:
+        case D3Q27System::BN:
             this->getNeighborsBottomNorth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TS:
+        case D3Q27System::TS:
             this->getNeighborsTopSouth(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TNE:
+        case D3Q27System::TNE:
             this->getNeighborsTopNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TNW:
+        case D3Q27System::TNW:
             this->getNeighborsTopNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TSE:
+        case D3Q27System::TSE:
             this->getNeighborsTopSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::TSW:
+        case D3Q27System::TSW:
             this->getNeighborsTopSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BNE:
+        case D3Q27System::BNE:
             this->getNeighborsBottomNorthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BNW:
+        case D3Q27System::BNW:
             this->getNeighborsBottomNorthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BSE:
+        case D3Q27System::BSE:
             this->getNeighborsBottomSouthEast(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::BSW:
+        case D3Q27System::BSW:
             this->getNeighborsBottomSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
-        case Grid3DSystem::REST:
+        case D3Q27System::DIR_000:
             this->getNeighborsZero(ix1, ix2, ix3, level, levelDepth, blocks);
             break;
         default:
@@ -1980,7 +1980,7 @@ void Grid3D::getBlocks(int level, int rank, bool active, std::vector<SPtr<Block3
 //////////////////////////////////////////////////////////////////////////
 int Grid3D::getFinestInitializedLevel()
 {
-    for (int i = Grid3DSystem::MAXLEVEL; i >= 0; i--)
+    for (int i = D3Q27System::MAXLEVEL; i >= 0; i--)
         if (this->levelSet[i].size() > 0)
             return (i);
     return (-1);
@@ -1988,7 +1988,7 @@ int Grid3D::getFinestInitializedLevel()
 //////////////////////////////////////////////////////////////////////////
 int Grid3D::getCoarsestInitializedLevel()
 {
-    for (int i = 0; i <= Grid3DSystem::MAXLEVEL; i++)
+    for (int i = 0; i <= D3Q27System::MAXLEVEL; i++)
         if (this->levelSet[i].size() > 0)
             return (i);
     return (-1);
@@ -2343,7 +2343,7 @@ void Grid3D::updateDistributedBlocks(std::shared_ptr<vf::mpi::Communicator> comm
             levelSet[l].clear();
         }
         this->levelSet.clear();
-        levelSet.resize(Grid3DSystem::MAXLEVEL + 1);
+        levelSet.resize(D3Q27System::MAXLEVEL + 1);
 
         int rsize = (int)blocks.size();
         for (int i = 0; i < rsize; i += 5) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
index 20851b019a3a0abd2c8865c7c40530e73bcf6245..4d4dff77128f3e0a254aee1fc937820bf35aa253 100644
--- a/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
@@ -91,7 +91,7 @@ void BGKLBMKernel::calculate(int step)
                     //////////////////////////////////////////////////////////////////////////
                     // read distribution
                     ////////////////////////////////////////////////////////////////////////////
-                    f[REST] = (*this->zeroDistributions)(x1, x2, x3);
+                    f[DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
 
                     f[E]   = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
                     f[N]   = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
@@ -122,7 +122,7 @@ void BGKLBMKernel::calculate(int step)
                     f[BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
                     //////////////////////////////////////////////////////////////////////////
 
-                    drho = f[REST] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + f[NE] + f[SW] + f[SE] + f[NW] + f[TE] +
+                    drho = f[DIR_000] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + f[NE] + f[SW] + f[SE] + f[NW] + f[TE] +
                            f[BW] + f[BE] + f[TW] + f[TN] + f[BS] + f[BN] + f[TS] + f[TNE] + f[TSW] + f[TSE] + f[TNW] +
                            f[BNE] + f[BSW] + f[BSE] + f[BNW];
 
@@ -137,7 +137,7 @@ void BGKLBMKernel::calculate(int step)
 
                     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-                    feq[REST] = c8o27 * (drho - cu_sq);
+                    feq[DIR_000] = c8o27 * (drho - cu_sq);
                     feq[E]    = c2o27 * (drho + 3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq);
                     feq[W]    = c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
                     feq[N]    = c2o27 * (drho + 3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq);
@@ -174,7 +174,7 @@ void BGKLBMKernel::calculate(int step)
                                          c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
 
                     // Relaxation
-                    f[REST] += (feq[REST] - f[REST]) * collFactor;
+                    f[DIR_000] += (feq[DIR_000] - f[DIR_000]) * collFactor;
                     f[E] += (feq[E] - f[E]) * collFactor;
                     f[W] += (feq[W] - f[W]) * collFactor;
                     f[N] += (feq[N] - f[N]) * collFactor;
@@ -214,7 +214,7 @@ void BGKLBMKernel::calculate(int step)
                         forcingX2 = muForcingX2.Eval();
                         forcingX3 = muForcingX3.Eval();
 
-                        f[REST] += 0.0;
+                        f[DIR_000] += 0.0;
                         f[E] += 3.0 * c2o27 * (forcingX1);
                         f[W] += 3.0 * c2o27 * (-forcingX1);
                         f[N] += 3.0 * c2o27 * (forcingX2);
@@ -244,7 +244,7 @@ void BGKLBMKernel::calculate(int step)
                     }
                     //////////////////////////////////////////////////////////////////////////
 #ifdef PROOF_CORRECTNESS
-                    LBMReal rho_post = f[REST] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + f[NE] + f[SW] + f[SE] +
+                    LBMReal rho_post = f[DIR_000] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + f[NE] + f[SW] + f[SE] +
                                        f[NW] + f[TE] + f[BW] + f[BE] + f[TW] + f[TN] + f[BS] + f[BN] + f[TS] + f[TNE] +
                                        f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW];
                     LBMReal dif = drho - rho_post;
@@ -291,7 +291,7 @@ void BGKLBMKernel::calculate(int step)
                     (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p)  = f[D3Q27System::INV_BNW];
                     (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p)   = f[D3Q27System::INV_BNE];
 
-                    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+                    (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
                     //////////////////////////////////////////////////////////////////////////
                 }
             }
diff --git a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.cpp
index 036b35379ec218585a43a67d03f2a03deb79d6e5..21b23b61d145530b5b36a73d9056039dea25d75e 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetInterpolationProcessor.cpp
@@ -497,7 +497,7 @@ void CompressibleOffsetInterpolationProcessor::calcInterpolatedNodeCF(LBMReal* f
    f[BSW]  = f_TNE  + xs*x_TNE  + ys*y_TNE  + zs*z_TNE  + xs*ys*xy_TNE  + xs*zs*xz_TNE  + ys*zs*yz_TNE  + feq[BSW];
    f[BSE]  = f_TNW  + xs*x_TNW  + ys*y_TNW  + zs*z_TNW  + xs*ys*xy_TNW  + xs*zs*xz_TNW  + ys*zs*yz_TNW  + feq[BSE];
    f[BNW]  = f_TSE  + xs*x_TSE  + ys*y_TSE  + zs*z_TSE  + xs*ys*xy_TSE  + xs*zs*xz_TSE  + ys*zs*yz_TSE  + feq[BNW];
-   f[REST] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[REST];
+   f[DIR_000] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 //Position SWB -0.25, -0.25, -0.25
@@ -691,7 +691,7 @@ void CompressibleOffsetInterpolationProcessor::calcInterpolatedNodeFC(LBMReal* f
    f[BNW]  = f_TSE  + feq[BNW];
    f[BSE]  = f_TNW  + feq[BSE];
    f[BSW]  = f_TNE  + feq[BSW];
-   f[REST] = f_ZERO + feq[REST];
+   f[DIR_000] = f_ZERO + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void CompressibleOffsetInterpolationProcessor::calcInterpolatedVelocity(LBMReal x, LBMReal y, LBMReal z, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3)
diff --git a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
index 4dec637580458cfa77d151b810df04a853116de8..936cf18ab89d7153b0621df61ac0023d324a6f43 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
@@ -803,7 +803,7 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeCF(LBM
    f[BS]   = mfbaa;
    f[BN]   = mfbca;
    f[TS]   = mfbac;
-   f[REST] = mfbbb;
+   f[DIR_000] = mfbbb;
    f[TNE]  = mfccc;
    f[TSE]  = mfcac;
    f[BNE]  = mfcca;
@@ -1251,7 +1251,7 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeFC(LBM
    f[BS]   = mfbaa;
    f[BN]   = mfbca;
    f[TS]   = mfbac;
-   f[REST] = mfbbb;
+   f[DIR_000] = mfbbb;
    f[TNE]  = mfccc;
    f[TSE]  = mfcac;
    f[BNE]  = mfcca;
diff --git a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
index 7a19f156e4447acd9d4451ce4c1a1de7bf5c990d..721dc5dad88bdec9a53ab6dc643e504767239398 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetSquarePressureInterpolationProcessor.cpp
@@ -798,7 +798,7 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    f[BS]   = mfbaa;
    f[BN]   = mfbca;
    f[TS]   = mfbac;
-   f[REST] = mfbbb;
+   f[DIR_000] = mfbbb;
    f[TNE]  = mfccc;
    f[TSE]  = mfcac;
    f[BNE]  = mfcca;
@@ -1251,7 +1251,7 @@ void CompressibleOffsetSquarePressureInterpolationProcessor::calcInterpolatedNod
    f[BS]   = mfbaa;
    f[BN]   = mfbca;
    f[TS]   = mfbac;
-   f[REST] = mfbbb;
+   f[DIR_000] = mfbbb;
    f[TNE]  = mfccc;
    f[TSE]  = mfcac;
    f[BNE]  = mfcca;
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
index b82616cbcd06de1dbdfb8d991a18266ee354e483..62adaa06b250f524266e23a88959154cd122bc42 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
@@ -18,14 +18,23 @@ const int DX1[] = { 0,  1, -1,  0,  0,  0,  0,  1, -1,  1, -1,  1, -1,  1, -1,
 const int DX2[] = { 0,  0,  0,  1, -1,  0,  0,  1, -1, -1,  1,  0,  0,  0,  0,  1, -1,  1, -1,   1,  1, -1, -1,  1,  1, -1, -1 };
 const int DX3[] = { 0,  0,  0,  0,  0,  1, -1,  0,  0,  0,  0,  1, -1, -1,  1,  1, -1, -1,  1,   1,  1,  1,  1, -1, -1, -1, -1 };
 
-const double WEIGTH[] = { c2o27,  c2o27,  c2o27,  c2o27,  c2o27,  c2o27,  c1o54,  c1o54,  c1o54,
-                          c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,
-                          c1o216, c1o216, c1o216, c1o216, c1o216, c1o216, c1o216, c1o216, c8o27 };
+const double WEIGTH[] = { c8o27,  
+                          c2o27,  c2o27,  c2o27,  c2o27,  c2o27,  c2o27,  
+                          c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,  c1o54,
+                          c1o216, c1o216, c1o216, c1o216, c1o216, c1o216, c1o216, c1o216 };
 
-const int INVDIR[] = { INV_E,   INV_W,   INV_N,   INV_S,   INV_T,   INV_B,   INV_NE,  INV_SW, INV_SE,
+const int INVDIR[] = { DIR_000, INV_E,   INV_W,   INV_N,   INV_S,   INV_T,   INV_B,   INV_NE,  INV_SW, INV_SE,
                        INV_NW,  INV_TE,  INV_BW,  INV_BE,  INV_TW,  INV_TN,  INV_BS,  INV_BN, INV_TS,
                        INV_TNE, INV_TNW, INV_TSE, INV_TSW, INV_BNE, INV_BNW, INV_BSE, INV_BSW };
 
+// index             0   1   2   3   4   5  6   7   8    9  10  11  12  13  14  15  16  17  18
+// direction:        E,  W,  N,  S,  T,  B, NE, SW, SE, NW, TE, BW, BE, TW, TN, BS, BN, TS, TNE TNW TSE TSW BNE BNW BSE
+// BSW
+//const int EX1[] = { 0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1 };
+//const int EX2[] = { 0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, -1 };
+//const int EX3[] = { 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1 };
+
+//////////////////////////////////////////////////////////////////////////
 
 
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index edfb903c1a5748b13d086bdd72c16b75368d3beb..5e80ee54f89366fc9d28b9294bf001a36f8252c9 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -64,6 +64,13 @@ extern const double WEIGTH[ENDDIR + 1];
 
 extern const double cNorm[3][ENDDIR];
 
+static const int MINLEVEL = 0;
+static const int MAXLEVEL = 25;
+
+extern const int EX1[ENDDIR + 1];
+extern const int EX2[ENDDIR + 1];
+extern const int EX3[ENDDIR + 1];
+
 //static const int E    = 0;
 //static const int W    = 1;
 //static const int N    = 2;
@@ -92,7 +99,35 @@ extern const double cNorm[3][ENDDIR];
 //static const int BSW  = 25;
 //static const int REST = 26;
 
-static constexpr int REST = 0;
+//static constexpr int REST = 0;
+//static constexpr int E = 1;
+//static constexpr int W = 2;
+//static constexpr int N = 3;
+//static constexpr int S = 4;
+//static constexpr int T = 5;
+//static constexpr int B = 6;
+//static constexpr int NE = 7;
+//static constexpr int SW = 8;
+//static constexpr int SE = 9;
+//static constexpr int NW = 10;
+//static constexpr int TE = 11;
+//static constexpr int BW = 12;
+//static constexpr int BE = 13;
+//static constexpr int TW = 14;
+//static constexpr int TN = 15;
+//static constexpr int BS = 16;
+//static constexpr int BN = 17;
+//static constexpr int TS = 18;
+//static constexpr int TNE = 19;
+//static constexpr int TNW = 20;
+//static constexpr int TSE = 21;
+//static constexpr int TSW = 22;
+//static constexpr int BNE = 23;
+//static constexpr int BNW = 24;
+//static constexpr int BSE = 25;
+//static constexpr int BSW = 26;
+
+static constexpr int DIR_000 = 0;
 static constexpr int E = 1;
 static constexpr int W = 2;
 static constexpr int N = 3;
@@ -176,6 +211,180 @@ static const int ET_BNW = 11;
 static const int ET_TSW = 12;
 static const int ET_BNE = 12;
 
+//////////////////////////////////////////////////////////////////////////
+inline std::string getDirectionString(int direction)
+{
+    switch (direction) {
+        case E:
+            return "E";
+        case W:
+            return "W";
+        case N:
+            return "N";
+        case S:
+            return "S";
+        case T:
+            return "T";
+        case B:
+            return "B";
+        case NE:
+            return "NE";
+        case NW:
+            return "NW";
+        case SE:
+            return "SE";
+        case SW:
+            return "SW";
+        case TE:
+            return "TE";
+        case TW:
+            return "TW";
+        case BE:
+            return "BE";
+        case BW:
+            return "BW";
+        case TN:
+            return "TN";
+        case TS:
+            return "TS";
+        case BN:
+            return "BN";
+        case BS:
+            return "BS";
+        case TNE:
+            return "TNE";
+        case TNW:
+            return "TNW";
+        case TSE:
+            return "TSE";
+        case TSW:
+            return "TSW";
+        case BNE:
+            return "BNE";
+        case BNW:
+            return "BNW";
+        case BSE:
+            return "BSE";
+        case BSW:
+            return "BSW";
+        default:
+            return "Cell3DSystem::getDrectionString(...) - unknown dir";
+    }
+}
+//////////////////////////////////////////////////////////////////////////
+static inline void setNeighborCoordinatesForDirection(int &x1, int &x2, int &x3, const int &direction)
+{
+    switch (direction) {
+        case D3Q27System::E:
+            x1++;
+            break;
+        case D3Q27System::N:
+            x2++;
+            break;
+        case D3Q27System::T:
+            x3++;
+            break;
+        case D3Q27System::W:
+            x1--;
+            break;
+        case D3Q27System::S:
+            x2--;
+            break;
+        case D3Q27System::B:
+            x3--;
+            break;
+        case D3Q27System::NE:
+            x1++;
+            x2++;
+            break;
+        case D3Q27System::NW:
+            x1--;
+            x2++;
+            break;
+        case D3Q27System::SW:
+            x1--;
+            x2--;
+            break;
+        case D3Q27System::SE:
+            x1++;
+            x2--;
+            break;
+        case D3Q27System::TE:
+            x1++;
+            x3++;
+            break;
+        case D3Q27System::BW:
+            x1--;
+            x3--;
+            break;
+        case D3Q27System::BE:
+            x1++;
+            x3--;
+            break;
+        case D3Q27System::TW:
+            x1--;
+            x3++;
+            break;
+        case D3Q27System::TN:
+            x2++;
+            x3++;
+            break;
+        case D3Q27System::BS:
+            x2--;
+            x3--;
+            break;
+        case D3Q27System::BN:
+            x2++;
+            x3--;
+            break;
+        case D3Q27System::TS:
+            x2--;
+            x3++;
+            break;
+        case D3Q27System::TNE:
+            x1++;
+            x2++;
+            x3++;
+            break;
+        case D3Q27System::TNW:
+            x1--;
+            x2++;
+            x3++;
+            break;
+        case D3Q27System::TSE:
+            x1++;
+            x2--;
+            x3++;
+            break;
+        case D3Q27System::TSW:
+            x1--;
+            x2--;
+            x3++;
+            break;
+        case D3Q27System::BNE:
+            x1++;
+            x2++;
+            x3--;
+            break;
+        case D3Q27System::BNW:
+            x1--;
+            x2++;
+            x3--;
+            break;
+        case D3Q27System::BSE:
+            x1++;
+            x2--;
+            x3--;
+            break;
+        case D3Q27System::BSW:
+            x1--;
+            x2--;
+            x3--;
+            break;
+        default:
+            throw UbException(UB_EXARGS, "no direction ...");
+    }
+}
 
 //////////////////////////////////////////////////////////////////////////
 // MACROSCOPIC VALUES
@@ -197,7 +406,7 @@ static void calcDensity(const LBMReal *const &f /*[27]*/, LBMReal &rho)
     rho = ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
           (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
            ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
-          ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+          ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000];
 }
 /*=====================================================================*/
 static void calcIncompVelocityX1(const LBMReal *const &f /*[27]*/, LBMReal &vx1)
@@ -289,7 +498,7 @@ static LBMReal getCompFeqForDirection(const int &direction, const LBMReal &drho,
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
     LBMReal rho   = drho + UbMath::one;
     switch (direction) {
-        case REST:
+        case DIR_000:
             return REAL_CAST(UbMath::c8o27 * (drho + rho * (-cu_sq)));
         case E:
             return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx1) + UbMath::c9o2 * (vx1) * (vx1)-cu_sq)));
@@ -382,7 +591,7 @@ static void calcCompFeq(LBMReal *const &feq /*[27]*/, const LBMReal &drho, const
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
     LBMReal rho   = drho + UbMath::one;
 
-    feq[REST] = UbMath::c8o27 * (drho + rho * (-cu_sq));
+    feq[DIR_000] = UbMath::c8o27 * (drho + rho * (-cu_sq));
     feq[E]    = UbMath::c2o27 * (drho + rho * (3.0 * (vx1) + UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
     feq[W]    = UbMath::c2o27 * (drho + rho * (3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq));
     feq[N]    = UbMath::c2o27 * (drho + rho * (3.0 * (vx2) + UbMath::c9o2 * (vx2) * (vx2)-cu_sq));
@@ -429,7 +638,7 @@ static LBMReal getIncompFeqForDirection(const int &direction, const LBMReal &drh
     LBMReal cu_sq = 1.5f * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
     switch (direction) {
-        case REST:
+        case DIR_000:
             return REAL_CAST(UbMath::c8o27 * (drho - cu_sq));
         case E:
             return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx1) + UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
@@ -513,7 +722,7 @@ static void calcIncompFeq(LBMReal *const &feq /*[27]*/, const LBMReal &drho, con
 {
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-    feq[REST] = UbMath::c8o27 * (drho - cu_sq);
+    feq[DIR_000] = UbMath::c8o27 * (drho - cu_sq);
     feq[E]    = UbMath::c2o27 * (drho + 3.0 * (vx1) + UbMath::c9o2 * (vx1) * (vx1)-cu_sq);
     feq[W]    = UbMath::c2o27 * (drho + 3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq);
     feq[N]    = UbMath::c2o27 * (drho + 3.0 * (vx2) + UbMath::c9o2 * (vx2) * (vx2)-cu_sq);
@@ -839,7 +1048,7 @@ static inline LBMReal getShearRate(const LBMReal *const f, LBMReal collFactorF)
     LBMReal mfaca = f[BNW];
     LBMReal mfcca = f[BNE];
 
-    LBMReal mfbbb = f[REST];
+    LBMReal mfbbb = f[DIR_000];
 
     LBMReal m0, m1, m2;
 
@@ -1143,7 +1352,7 @@ static void calcMultiphaseFeq(LBMReal *const &feq /*[27]*/, const LBMReal &rho,
     using namespace UbMath;
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-    feq[REST] = c8o27 * (p1 + rho * c1o3 * (-cu_sq));
+    feq[DIR_000] = c8o27 * (p1 + rho * c1o3 * (-cu_sq));
     feq[E]    = c2o27 * (p1 + rho * c1o3 * (3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq));
     feq[W]    = c2o27 * (p1 + rho * c1o3 * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
     feq[N]    = c2o27 * (p1 + rho * c1o3 * (3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq));
@@ -1186,7 +1395,7 @@ static void calcMultiphaseFeqVB(LBMReal *const &feq /*[27]*/, const LBMReal &p1,
     using namespace UbMath;
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-    feq[REST] = p1 + c8o27 * (-cu_sq);
+    feq[DIR_000] = p1 + c8o27 * (-cu_sq);
     feq[E]    = c2o27 * ((3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq));
     feq[W]    = c2o27 * ((3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
     feq[N]    = c2o27 * ((3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq));
@@ -1221,7 +1430,7 @@ static void calcMultiphaseHeq(LBMReal *const &heq /*[27]*/, const LBMReal &phi,
     using namespace UbMath;
     LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-    heq[REST] = c8o27 * phi * (1.0 - cu_sq);
+    heq[DIR_000] = c8o27 * phi * (1.0 - cu_sq);
     heq[E]    = c2o27 * phi * (1.0 + 3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq);
     heq[W]    = c2o27 * phi * (1.0 + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
     heq[N]    = c2o27 * phi * (1.0 + 3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq);
@@ -1249,7 +1458,6 @@ static void calcMultiphaseHeq(LBMReal *const &heq /*[27]*/, const LBMReal &phi,
     heq[BSE] = c1o216 * phi * (1.0 + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
     heq[TNW] = c1o216 * phi * (1.0 + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
 }
-//////////////////////////////////////////////////////////////////////////
 
 } // namespace D3Q27System
 #endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
index 15e6f1dddb88f31e7bf57e2d3235e04b48da1080..b4ccea0c996ba29d099ca831d6b2b2bad9ade52e 100644
--- a/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
@@ -569,7 +569,7 @@ void IncompressibleOffsetInterpolationProcessor::calcInterpolatedNode(LBMReal* f
    f[BSW]  = f_TNE  + xs*x_TNE  + ys*y_TNE  + zs*z_TNE  + xs*ys*xy_TNE  + xs*zs*xz_TNE  + ys*zs*yz_TNE  + feq[BSW];
    f[BSE]  = f_TNW  + xs*x_TNW  + ys*y_TNW  + zs*z_TNW  + xs*ys*xy_TNW  + xs*zs*xz_TNW  + ys*zs*yz_TNW  + feq[BSE];
    f[BNW]  = f_TSE  + xs*x_TSE  + ys*y_TSE  + zs*z_TSE  + xs*ys*xy_TSE  + xs*zs*xz_TSE  + ys*zs*yz_TSE  + feq[BNW];
-   f[REST] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[REST];
+   f[DIR_000] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 //Position SWB -0.25, -0.25, -0.25
@@ -763,7 +763,7 @@ void IncompressibleOffsetInterpolationProcessor::calcInterpolatedNodeFC(LBMReal*
    f[BNW]  = f_TSE  + feq[BNW];
    f[BSE]  = f_TNW  + feq[BSE];
    f[BSW]  = f_TNE  + feq[BSW];
-   f[REST] = f_ZERO + feq[REST];
+   f[DIR_000] = f_ZERO + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void IncompressibleOffsetInterpolationProcessor::calcInterpolatedVelocity(LBMReal x, LBMReal y, LBMReal z, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3)
diff --git a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
index 554fc6614c11e2e8b96f5829d81d1c27e9365870..56074586c1da94da2382b439b5ccf6c50aa88187 100644
--- a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
@@ -894,7 +894,7 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                //////////////////////////////////////////////////////////////////////////
                //read distribution
                ////////////////////////////////////////////////////////////////////////////
-               f[REST] = (*this->zeroDistributions)(x1, x2, x3);
+               f[DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
 
                f[E] = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
                f[N] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
@@ -928,7 +928,7 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                drho = ((f[TNE]+f[BSW])+(f[TSE]+f[BNW]))+((f[BSE]+f[TNW])+(f[TSW]+f[BNE]))
                   +(((f[NE]+f[SW])+(f[SE]+f[NW]))+((f[TE]+f[BW])+(f[BE]+f[TW]))
                      +((f[BN]+f[TS])+(f[TN]+f[BS])))+((f[E]+f[W])+(f[N]+f[S])
-                        +(f[T]+f[B]))+f[REST];
+                        +(f[T]+f[B]))+f[DIR_000];
 
                //vx1 = ((((f[TNE]-f[BSW])+(f[TSE]-f[BNW]))+((f[BSE]-f[TNW])+(f[BNE]-f[TSW])))+
                //   (((f[BE]-f[TW])+(f[TE]-f[BW]))+((f[SE]-f[NW])+(f[NE]-f[SW])))+
@@ -956,7 +956,7 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
 
                LBMReal cu_sq = 1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
 
-               feq[REST] = c8o27*(drho-cu_sq);
+               feq[DIR_000] = c8o27*(drho-cu_sq);
                feq[E] = c2o27*(drho+3.0*(vx1)+c9o2*(vx1)*(vx1)-cu_sq);
                feq[W] = c2o27*(drho+3.0*(-vx1)+c9o2*(-vx1)*(-vx1)-cu_sq);
                feq[N] = c2o27*(drho+3.0*(vx2)+c9o2*(vx2)*(vx2)-cu_sq);
@@ -985,7 +985,7 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                feq[TNW] = c1o216*(drho+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
 
                //Relaxation
-               f[REST] += (feq[REST]-f[REST])*collFactor;
+               f[DIR_000] += (feq[DIR_000]-f[DIR_000])*collFactor;
                f[E] += (feq[E]-f[E])*collFactor;
                f[W] += (feq[W]-f[W])*collFactor;
                f[N] += (feq[N]-f[N])*collFactor;
@@ -1061,7 +1061,7 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[D3Q27System::INV_BNW];
                (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = f[D3Q27System::INV_BNE];
 
-               (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+               (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::DIR_000];
                //////////////////////////////////////////////////////////////////////////
 
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp b/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
index 6076eb018097dc77afcf37af2d14206325be463c..e48e846a4b69ed6d074b0729696dc321531d7808 100644
--- a/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
@@ -88,7 +88,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                //////////////////////////////////////////////////////////////////////////
                //read distribution
                ////////////////////////////////////////////////////////////////////////////
-               f[REST] = (*this->zeroDistributions)(x1,x2,x3);
+               f[DIR_000] = (*this->zeroDistributions)(x1,x2,x3);
 
                f[E] = (*this->localDistributions)(D3Q27System::ET_E, x1,x2,x3);
                f[N] = (*this->localDistributions)(D3Q27System::ET_N,x1,x2,x3);  
@@ -119,7 +119,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                f[BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE,x1,x2,x3p);
                //////////////////////////////////////////////////////////////////////////
 
-               drho = f[REST] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] 
+               drho = f[DIR_000] + f[E] + f[W] + f[N] + f[S] + f[T] + f[B] 
                + f[NE] + f[SW] + f[SE] + f[NW] + f[TE] + f[BW] + f[BE]
                + f[TW] + f[TN] + f[BS] + f[BN] + f[TS] + f[TNE] + f[TSW]
                + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW];
@@ -138,7 +138,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
 
                LBMReal cu_sq=1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
 
-               feq[REST] =  c8o27*(drho-cu_sq);
+               feq[DIR_000] =  c8o27*(drho-cu_sq);
                feq[E] =  c2o27*(drho+3.0*( vx1   )+c9o2*( vx1   )*( vx1   )-cu_sq);
                feq[W] =  c2o27*(drho+3.0*(-vx1   )+c9o2*(-vx1   )*(-vx1   )-cu_sq);
                feq[N] =  c2o27*(drho+3.0*(    vx2)+c9o2*(    vx2)*(    vx2)-cu_sq);
@@ -167,7 +167,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                feq[TNW]= c1o216*(drho+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
 
                //Relaxation
-               f[REST] += (feq[REST]-f[REST])*collFactor;
+               f[DIR_000] += (feq[DIR_000]-f[DIR_000])*collFactor;
                f[E] += (feq[E]-f[E])*collFactor;
                f[W] += (feq[W]-f[W])*collFactor;
                f[N] += (feq[N]-f[N])*collFactor;
@@ -208,7 +208,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                   forcingX2 = muForcingX2.Eval();
                   forcingX3 = muForcingX3.Eval();
 
-                  f[REST] +=                   0.0                        ;
+                  f[DIR_000] +=                   0.0                        ;
                   f[E  ] +=  3.0*c2o27  *  (forcingX1)                    ;
                   f[W  ] +=  3.0*c2o27  *  (-forcingX1)                   ;
                   f[N  ] +=  3.0*c2o27  *             (forcingX2)         ;
@@ -283,7 +283,7 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                (*this->nonLocalDistributions)(D3Q27System::ET_BNW,x1p,x2,  x3p) = f[D3Q27System::INV_BNW];
                (*this->nonLocalDistributions)(D3Q27System::ET_BNE,x1,  x2,  x3p) = f[D3Q27System::INV_BNE];
 
-               (*this->zeroDistributions)(x1,x2,x3) = f[D3Q27System::REST];
+               (*this->zeroDistributions)(x1,x2,x3) = f[D3Q27System::DIR_000];
                //////////////////////////////////////////////////////////////////////////
 
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
index 5cfee04e4a26129d0c27bc0d98c70a5ae64c3825..9de46c4f97be387f4c17a109f6630875a7bd1255 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
@@ -235,13 +235,13 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         LBMReal dX3_phi = gradX3_phi();
 
                         LBMReal denom = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi) + 1e-9;
-                        collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+                        collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
-                        LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                        LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
                         if (withForcing) {
                             // muX1 = static_cast<double>(x1-1+ix1*maxX1);
@@ -308,9 +308,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                                                (DX3[dir]) * (fac1 * dX3_phi + gamma * (mu * dX3_phi + forcingX3));
                         }
 
-                        LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
-                        LBMReal fac1      = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
-                        forcingTerm[REST] = (-ux) * (fac1 * dX1_phi + gamma * (mu * dX1_phi + forcingX1)) +
+                        LBMReal gamma = WEIGTH[DIR_000] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
+                        LBMReal fac1      = (gamma - WEIGTH[DIR_000]) * c1o3 * rhoToPhi;
+                        forcingTerm[DIR_000] = (-ux) * (fac1 * dX1_phi + gamma * (mu * dX1_phi + forcingX1)) +
                                             (-uy) * (fac1 * dX2_phi + gamma * (mu * dX2_phi + forcingX2)) +
                                             (-uz) * (fac1 * dX3_phi + gamma * (mu * dX3_phi + forcingX3));
 
@@ -342,7 +342,7 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         mfcaa = 3.0 * (mfcaa + 0.5 * forcingTerm[BSE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BSE];
                         mfaca = 3.0 * (mfaca + 0.5 * forcingTerm[BNW]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNW];
                         mfcca = 3.0 * (mfcca + 0.5 * forcingTerm[BNE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNE];
-                        mfbbb = 3.0 * (mfbbb + 0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST];
+                        mfbbb = 3.0 * (mfbbb + 0.5 * forcingTerm[DIR_000]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST];
 
                         LBMReal rho1 = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca) +
                                        (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) +
@@ -1050,7 +1050,7 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         mfcaa = rho * c1o3 * (mfcaa) + 0.5 * forcingTerm[BSE];
                         mfaca = rho * c1o3 * (mfaca) + 0.5 * forcingTerm[BNW];
                         mfcca = rho * c1o3 * (mfcca) + 0.5 * forcingTerm[BNE];
-                        mfbbb = rho * c1o3 * (mfbbb) + 0.5 * forcingTerm[REST];
+                        mfbbb = rho * c1o3 * (mfbbb) + 0.5 * forcingTerm[DIR_000];
 
                         //////////////////////////////////////////////////////////////////////////
                         // write distribution for F
@@ -1119,24 +1119,24 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         h[BNW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p);
                         h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                        h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
+                        h[DIR_000] = (*this->zeroDistributionsH)(x1, x2, x3);
 
                         for (int dir = STARTF; dir < (ENDF + 1); dir++) {
                             LBMReal velProd = DX1[dir] * ux + DX2[dir] * uy + DX3[dir] * uz;
                             LBMReal velSq1  = velProd * velProd;
                             LBMReal hEq; //, gEq;
 
-                            if (dir != REST) {
+                            if (dir != DIR_000) {
                                 LBMReal dirGrad_phi = (phi[dir] - phi[INVDIR[dir]]) / 2.0;
-                                LBMReal hSource     = (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST]) * (dirGrad_phi) / denom; 
-                                hEq = phi[REST] * WEIGTH[dir] * (1.0 + 3.0 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)) +                                 hSource * WEIGTH[dir];
+                                LBMReal hSource     = (tauH - 0.5) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * (dirGrad_phi) / denom; 
+                                hEq = phi[DIR_000] * WEIGTH[dir] * (1.0 + 3.0 * velProd + 4.5 * velSq1 - 1.5 * (ux2 + uy2 + uz2)) +                                 hSource * WEIGTH[dir];
 
                                 // This corresponds with the collision factor of 1.0 which equals (tauH + 0.5).
                                 h[dir] = h[dir] - (h[dir] - hEq) / (tauH); 
 
                             } else {
-                                hEq = phi[REST] * WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
-                                h[REST] = h[REST] - (h[REST] - hEq) / (tauH); 
+                                hEq = phi[DIR_000] * WEIGTH[DIR_000] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
+                                h[DIR_000] = h[DIR_000] - (h[DIR_000] - hEq) / (tauH); 
                             }
                         }
 
@@ -1168,7 +1168,7 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p)  = h[D3Q27System::INV_BNW];
                         (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p)   = h[D3Q27System::INV_BNE];
 
-                        (*this->zeroDistributionsH)(x1, x2, x3) = h[D3Q27System::REST];
+                        (*this->zeroDistributionsH)(x1, x2, x3) = h[D3Q27System::DIR_000];
 
                         /////////////////////   END OF OLD BGK SOLVER ///////////////////////////////
                     }
@@ -1215,7 +1215,7 @@ LBMReal MultiphaseCumulantLBMKernel::nabla2_phi()
     using namespace D3Q27System;
     LBMReal sum = 0.0;
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
-        sum += WEIGTH[k] * (phi[k] - phi[REST]);
+        sum += WEIGTH[k] * (phi[k] - phi[DIR_000]);
     }
     return 6.0 * sum;
 }
@@ -1270,7 +1270,7 @@ void MultiphaseCumulantLBMKernel::computePhasefield()
                     h[BNW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p);
                     h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                    h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
+                    h[DIR_000] = (*this->zeroDistributionsH)(x1, x2, x3);
                 }
             }
         }
@@ -1284,7 +1284,7 @@ void MultiphaseCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal, IndexerX3X2X1
 
     SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-    phi[REST] = (*ph)(x1, x2, x3);
+    phi[DIR_000] = (*ph)(x1, x2, x3);
 
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
index 01fb23625afbf4700c337cf6ee7194a46f0bc895..50885d680735b065beeb87912bb4647de2a58b8d 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
@@ -383,13 +383,13 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::calculate(int step)
 
 
 
-					collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+					collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-					LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+					LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
 					//----------- Calculating Macroscopic Values -------------
-					LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH); //Incompressible
+					LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH); //Incompressible
 					//LBMReal rho = rhoL + (rhoH - rhoL) * phi[REST] + (one - phi[REST]) * (*pressure)(x1, x2, x3) * three; //compressible
 
 					LBMReal m0, m1, m2;
@@ -457,7 +457,7 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::calculate(int step)
 					}
 
 					//Viscosity increase by pressure gradient
-					LBMReal errPhi = (((1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale)- denom);
+					LBMReal errPhi = (((1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale)- denom);
 					//LBMReal limVis = 0.0000001*10;//0.01;
 					// collFactorM =collFactorM/(c1+limVis*(errPhi*errPhi)*collFactorM);
 					// collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
@@ -1333,11 +1333,11 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::calculate(int step)
 
 						// collision of 1st order moments
 						cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-							normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 						cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-							normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 						cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-							normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 						cx2 = cx * cx;
 						cy2 = cy * cy;
@@ -1524,17 +1524,17 @@ LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::nabla2_phi()
 {
 	using namespace D3Q27System;
 	LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-		(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+		(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-		((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+		((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
 
 	return 6.0 * sum;
@@ -1590,7 +1590,7 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::computePhasefield()
 					h[BNW] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
 					h[BNE] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-					h[REST] = (*this->zeroDistributionsH1)(x1, x2, x3);
+					h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
 				}
 			}
 		}
@@ -1604,7 +1604,7 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::findNeighbors(CbArray3D<L
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi[REST] = (*ph)(x1, x2, x3);
+	phi[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
@@ -1624,7 +1624,7 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::findNeighbors2(CbArray3D<
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi2[REST] = (*ph)(x1, x2, x3);
+	phi2[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
index 64b3cf1685561f3ef04ffc2662268eea83694fbf..2b68a072da69ffb97a9c597675961d676c0ff2ed 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
@@ -375,17 +375,17 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 					LBMReal normX2 = dX2_phi / denom;
 					LBMReal normX3 = dX3_phi / denom;
 
-					dX1_phi = normX1 * (1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale;
-                    dX2_phi = normX2 * (1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale;
-                    dX3_phi = normX3 * (1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale;
+					dX1_phi = normX1 * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale;
+                    dX2_phi = normX2 * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale;
+                    dX3_phi = normX3 * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale;
 
-					collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+					collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-					LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+					LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
 					//----------- Calculating Macroscopic Values -------------
-					LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+					LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
 					LBMReal m0, m1, m2;
 					LBMReal rhoRef=c1;
@@ -452,7 +452,7 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 					}
 
 					//Viscosity increase by pressure gradient
-					LBMReal errPhi = (((1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale)- denom);
+					LBMReal errPhi = (((1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale)- denom);
 					//LBMReal limVis = 0.0000001*10;//0.01;
 					// collFactorM =collFactorM/(c1+limVis*(errPhi*errPhi)*collFactorM);
 					// collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
@@ -1482,11 +1482,11 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
 
 						// collision of 1st order moments
 						cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-							normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 						cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-							normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 						cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-							normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+							normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 						cx2 = cx * cx;
 						cy2 = cy * cy;
@@ -1649,17 +1649,17 @@ LBMReal MultiphasePressureFilterLBMKernel::nabla2_phi()
 {
 	using namespace D3Q27System;
 	LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-		(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+		(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-		((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+		((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
 
 	return 6.0 * sum;
@@ -1715,7 +1715,7 @@ void MultiphasePressureFilterLBMKernel::computePhasefield()
 					h[BNW] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
 					h[BNE] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-					h[REST] = (*this->zeroDistributionsH1)(x1, x2, x3);
+					h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
 				}
 			}
 		}
@@ -1729,7 +1729,7 @@ void MultiphasePressureFilterLBMKernel::findNeighbors(CbArray3D<LBMReal, Indexer
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi[REST] = (*ph)(x1, x2, x3);
+	phi[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index 77695179c9d3750eed80f0bf56e57103a691df9c..f251dfafa1657ef3b3ab66c367bf678a6ce35b2c 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -431,13 +431,13 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 
 						///!test
 
-						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 						//collFactorM = phi[REST] - phiL < (phiH - phiL) * 0.05 ? collFactorG : collFactorL;
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
-                        LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                        LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
 						if (withForcing) {
 							// muX1 = static_cast<double>(x1-1+ix1*maxX1);
@@ -588,9 +588,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 				  // vvxF = vvxF;
 			   //}
 			   LBMReal weightGrad =  1.0-denom*denom/(denom*denom+0.0001*0.001);
-			   LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX1;
-			   LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX2;
-			   LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX3;
+			   LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX1;
+			   LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX2;
+			   LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX3;
 
 			   //dX1_phiF *= 1.2;
 			   //dX2_phiF *= 1.2;
@@ -646,9 +646,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 
 			   }
 
-			   LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
-			   LBMReal fac1 = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
-			   forcingTerm[REST] = (-vvxF) * (fac1 * dX1_phiF ) +
+			   LBMReal gamma = WEIGTH[DIR_000] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
+			   LBMReal fac1 = (gamma - WEIGTH[DIR_000]) * c1o3 * rhoToPhi;
+			   forcingTerm[DIR_000] = (-vvxF) * (fac1 * dX1_phiF ) +
 				   (-vvyF) * (fac1 * dX2_phiF ) +
 				   (-vvzF) * (fac1 * dX3_phiF );
 
@@ -774,7 +774,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   mfcaa +=3.0 * ( 0.5 * forcingTerm[BSE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BSE];
 			   mfaca +=3.0 * ( 0.5 * forcingTerm[BNW]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNW];
 			   mfcca +=3.0 * ( 0.5 * forcingTerm[BNE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNE];
-			   mfbbb +=3.0 * ( 0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+			   mfbbb +=3.0 * ( 0.5 * forcingTerm[DIR_000]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
 
 			   //--------------------------------------------------------
 
@@ -1162,7 +1162,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 
 			   ////relax unfiltered
 			   //! divergenceFilter 10.05.2021
-			   LBMReal divMag= (1.0 - phi[REST]) * (phi[REST])*10*5*sqrt(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
+			   LBMReal divMag= (1.0 - phi[DIR_000]) * (phi[DIR_000])*10*5*sqrt(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
 			  // LBMReal divMag = 500 *500* 50*(fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))))* (fabs((OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz))));
 			   //LBMReal divMag = (dX1_phi * dxux) > 0 ? (dX1_phi * dxux) : 0;
 			   //divMag += (dX2_phi * dyuy) > 0 ? (dX2_phi * dyuy) : 0;
@@ -1645,7 +1645,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   mfcaa += 3.0 * (0.5 * forcingTerm[BSE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BSE];
 			   mfaca += 3.0 * (0.5 * forcingTerm[BNW]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNW];
 			   mfcca += 3.0 * (0.5 * forcingTerm[BNE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNE];
-			   mfbbb += 3.0 * (0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+			   mfbbb += 3.0 * (0.5 * forcingTerm[DIR_000]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
 
 
 
@@ -2671,9 +2671,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   LBMReal Mccb = mfccb - mfaab * c1o9;
 
 			   // collision of 1st order moments
-			   cx = cx * (c1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
-			   cy = cy * (c1 - omegaD) + omegaD * vvy * concentration + normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
-			   cz = cz * (c1 - omegaD) + omegaD * vvz * concentration + normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+			   cx = cx * (c1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
+			   cy = cy * (c1 - omegaD) + omegaD * vvy * concentration + normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
+			   cz = cz * (c1 - omegaD) + omegaD * vvz * concentration + normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 			   //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
 			   //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -2961,17 +2961,17 @@ LBMReal MultiphaseScratchCumulantLBMKernel::nabla2_phi()
 {
     using namespace D3Q27System;
     LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-			(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+			(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-			((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3029,7 +3029,7 @@ void MultiphaseScratchCumulantLBMKernel::computePhasefield()
                     h[BNW] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNW, x1p, x2, x3p);
                     h[BNE] = (*this->nonLocalDistributionsH)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                    h[REST] = (*this->zeroDistributionsH)(x1, x2, x3);
+                    h[DIR_000] = (*this->zeroDistributionsH)(x1, x2, x3);
                 }
             }
         }
@@ -3043,7 +3043,7 @@ void MultiphaseScratchCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal, Indexe
 
     SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-    phi[REST] = (*ph)(x1, x2, x3);
+    phi[DIR_000] = (*ph)(x1, x2, x3);
 
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
index f226b8ecfa3c5654801561248f09e0251d2fa2ac..4212683ed33f444ac56a0d5e7049e7a19868e9b8 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
@@ -346,13 +346,13 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 
 
 
-						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
-                        LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                        LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
                             			   ////Incompressible Kernal
 
@@ -446,9 +446,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 				  // vvxF = vvxF;
 			   //}
 			   LBMReal weightGrad = 1.0;// -denom * denom / (denom * denom + 0.0001 * 0.001);
-			   LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX1;
-			   LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX2;
-			   LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[REST]) * (phi[REST]) * normX3;
+			   LBMReal dX1_phiF = dX1_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX1;
+			   LBMReal dX2_phiF = dX2_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX2;
+			   LBMReal dX3_phiF = dX3_phi * weightGrad + (1.0 - weightGrad) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * normX3;
 
 			   //dX1_phiF *= 1.2;
 			   //dX2_phiF *= 1.2;
@@ -498,9 +498,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 
 			   }
 
-			   LBMReal gamma = WEIGTH[REST] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
-			   LBMReal fac1 = (gamma - WEIGTH[REST]) * c1o3 * rhoToPhi;
-			   forcingTerm[REST] =	 (-vvxF) * (fac1 * (dX1_phiF * rhoH + dX2_phi2 * rhoL)) +
+			   LBMReal gamma = WEIGTH[DIR_000] * (1.0 - 1.5 * (ux2 + uy2 + uz2));
+			   LBMReal fac1 = (gamma - WEIGTH[DIR_000]) * c1o3 * rhoToPhi;
+			   forcingTerm[DIR_000] =	 (-vvxF) * (fac1 * (dX1_phiF * rhoH + dX2_phi2 * rhoL)) +
 				   (-vvyF) * (fac1 * (dX2_phiF * rhoH + dX2_phi2 * rhoL)) +
 				   (-vvzF) * (fac1 * (dX3_phiF * rhoH + dX3_phi2 * rhoL));
 
@@ -626,7 +626,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   mfcaa += 3.0 * (0.5 * forcingTerm[BSE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BSE];
 			   mfaca += 3.0 * (0.5 * forcingTerm[BNW]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNW];
 			   mfcca += 3.0 * (0.5 * forcingTerm[BNE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNE];
-			   mfbbb += 3.0 * (0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+			   mfbbb += 3.0 * (0.5 * forcingTerm[DIR_000]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
 
 			   //--------------------------------------------------------
 
@@ -1370,7 +1370,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   mfcaa += 3.0 * (0.5 * forcingTerm[BSE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BSE];
 			   mfaca += 3.0 * (0.5 * forcingTerm[BNW]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNW];
 			   mfcca += 3.0 * (0.5 * forcingTerm[BNE]) / rho;  //-(3.0*p1 - rho)*WEIGTH[BNE];
-			   mfbbb += 3.0 * (0.5 * forcingTerm[REST]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
+			   mfbbb += 3.0 * (0.5 * forcingTerm[DIR_000]) / rho; //- (3.0*p1 - rho)*WEIGTH[REST]
 
 
 
@@ -2422,11 +2422,11 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   
 
                cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 			   //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
 			   //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -2735,11 +2735,11 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 
    // collision of 1st order moments
    cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-	   normX1 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX1 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-	   normX2 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX2 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-	   normX3 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX3 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
    //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
    //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -3067,17 +3067,17 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::nabla2_phi()
 {
     using namespace D3Q27System;
     LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-			(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+			(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-			((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3135,7 +3135,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::computePhasefield()
                     h[BNW] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
                     h[BNE] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                    h[REST] = (*this->zeroDistributionsH1)(x1, x2, x3);
+                    h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
                 }
             }
         }
@@ -3149,7 +3149,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::findNeighbors(CbArray3D<LBMReal,
 
     SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-    phi[REST] = (*ph)(x1, x2, x3);
+    phi[DIR_000] = (*ph)(x1, x2, x3);
 
 
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
@@ -3169,7 +3169,7 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::findNeighbors2(CbArray3D<LBMReal
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi2[REST] = (*ph)(x1, x2, x3);
+	phi2[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
index 179fc6f4d398be89481336af0c753b699828d2ee..b2eddb8746303897b9d0c88087cb4c0cc5bc2b94 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
@@ -567,13 +567,13 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 
 
 
-						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
-                        LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                        LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
 						//! variable density -> TRANSFER!
 						//LBMReal rho = rhoH * ((*phaseField)(x1, x2, x3)) + rhoL * ((*phaseField2)(x1, x2, x3));
@@ -712,7 +712,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 			   }
 
 			   //Viscosity increase by pressure gradient
-			   LBMReal errPhi = (((1.0 - phi[REST]) * (phi[REST]) * oneOverInterfaceScale)- denom);
+			   LBMReal errPhi = (((1.0 - phi[DIR_000]) * (phi[DIR_000]) * oneOverInterfaceScale)- denom);
 			   //LBMReal limVis = 0.0000001*10;//0.01;
 			  // collFactorM =collFactorM/(c1+limVis*(errPhi*errPhi)*collFactorM);
 			  // collFactorM = (collFactorM < 1.8) ? 1.8 : collFactorM;
@@ -2780,11 +2780,11 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 			   
 
                cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 			   //cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
 				  // normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST])*(phi[REST]+phi2[REST]) * c1o3 * oneOverInterfaceScale;
@@ -3100,11 +3100,11 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::calculate(int step)
 
    // collision of 1st order moments
    cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-	   normX1 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX1 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-	   normX2 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX2 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-	   normX3 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX3 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
    //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
    //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -3432,17 +3432,17 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::nabla2_phi()
 {
     using namespace D3Q27System;
     LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-			(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+			(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-			((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3500,7 +3500,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::computePhasefield()
                     h[BNW] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
                     h[BNE] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                    h[REST] = (*this->zeroDistributionsH1)(x1, x2, x3);
+                    h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
                 }
             }
         }
@@ -3514,7 +3514,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::findNeighbors(CbArray3D<LB
 
     SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-    phi[REST] = (*ph)(x1, x2, x3);
+    phi[DIR_000] = (*ph)(x1, x2, x3);
 
 
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
@@ -3534,7 +3534,7 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::findNeighbors2(CbArray3D<L
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi2[REST] = (*ph)(x1, x2, x3);
+	phi2[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
index 294568b4f8b845e3c510e77e50eb780c1a739e4b..b10c528bd63e08e37a0399c3b7a919970781c1a3 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
@@ -295,7 +295,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::calculate(int step)
 							+ (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
 							+ (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
 
-						LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+						LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 						(*pressure)(x1, x2, x3) = (*pressure)(x1, x2, x3) + rho * c1o3 * drho;
 
 						////!!!!!! relplace by pointer swap!
@@ -545,13 +545,13 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::calculate(int step)
 
 
 
-						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
+						collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[DIR_000] - phiH) / (phiH - phiL);
 
 
-                        LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+                        LBMReal mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
                         //----------- Calculating Macroscopic Values -------------
-                        LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+                        LBMReal rho = rhoH + rhoToPhi * (phi[DIR_000] - phiH);
 
                             			   ////Incompressible Kernal
 
@@ -2719,11 +2719,11 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::calculate(int step)
 			   
 
                cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX1 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX2 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
                cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[REST]) * (phi[REST]) * c1o3 * oneOverInterfaceScale;
+                    normX3 * (c1 - 0.5 * omegaD) * (1.0 - phi[DIR_000]) * (phi[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
 			   //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
 			   //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -3032,11 +3032,11 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::calculate(int step)
 
    // collision of 1st order moments
    cx = cx * (c1 - omegaD) + omegaD * vvx * concentration +
-	   normX1 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX1 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cy = cy * (c1 - omegaD) + omegaD * vvy * concentration +
-	   normX2 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX2 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
    cz = cz * (c1 - omegaD) + omegaD * vvz * concentration +
-	   normX3 * (c1 - 0.5 * omegaD) * ( phi[REST]) * (phi2[REST]) * c1o3 * oneOverInterfaceScale;
+	   normX3 * (c1 - 0.5 * omegaD) * ( phi[DIR_000]) * (phi2[DIR_000]) * c1o3 * oneOverInterfaceScale;
 
    //mhx = (ux * phi[REST] + normX1 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhx;
    //mhy = (uy * phi[REST] + normX2 * (tauH - 0.5) * (1.0 - phi[REST]) * (phi[REST])) / tauH + (1.0 - 1.0 / tauH) * mhy;
@@ -3364,17 +3364,17 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::nabla2_phi()
 {
     using namespace D3Q27System;
     LBMReal sum = 0.0;
-	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[REST]) + (phi[BSW] - phi[REST])) + ((phi[TSW] - phi[REST]) + (phi[BNE] - phi[REST])))
-		+ (((phi[TNW] - phi[REST]) + (phi[BSE] - phi[REST])) + ((phi[TSE] - phi[REST]) + (phi[BNW] - phi[REST]))));
+	sum += WEIGTH[TNE] * ((((phi[TNE] - phi[DIR_000]) + (phi[BSW] - phi[DIR_000])) + ((phi[TSW] - phi[DIR_000]) + (phi[BNE] - phi[DIR_000])))
+		+ (((phi[TNW] - phi[DIR_000]) + (phi[BSE] - phi[DIR_000])) + ((phi[TSE] - phi[DIR_000]) + (phi[BNW] - phi[DIR_000]))));
 	sum += WEIGTH[TN] * (
-			(((phi[TN] - phi[REST]) + (phi[BS] - phi[REST])) + ((phi[TS] - phi[REST]) + (phi[BN] - phi[REST])))
-		+	(((phi[TE] - phi[REST]) + (phi[BW] - phi[REST])) + ((phi[TW] - phi[REST]) + (phi[BE] - phi[REST])))
-		+	(((phi[NE] - phi[REST]) + (phi[SW] - phi[REST])) + ((phi[NW] - phi[REST]) + (phi[SE] - phi[REST])))
+			(((phi[TN] - phi[DIR_000]) + (phi[BS] - phi[DIR_000])) + ((phi[TS] - phi[DIR_000]) + (phi[BN] - phi[DIR_000])))
+		+	(((phi[TE] - phi[DIR_000]) + (phi[BW] - phi[DIR_000])) + ((phi[TW] - phi[DIR_000]) + (phi[BE] - phi[DIR_000])))
+		+	(((phi[NE] - phi[DIR_000]) + (phi[SW] - phi[DIR_000])) + ((phi[NW] - phi[DIR_000]) + (phi[SE] - phi[DIR_000])))
 		);
 	sum += WEIGTH[T] * (
-			((phi[T] - phi[REST]) + (phi[B] - phi[REST]))
-		+	((phi[N] - phi[REST]) + (phi[S] - phi[REST]))
-		+	((phi[E] - phi[REST]) + (phi[W] - phi[REST]))
+			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
+		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
+		+	((phi[E] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3432,7 +3432,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::computePhasefield()
                     h[BNW] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
                     h[BNE] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
 
-                    h[REST] = (*this->zeroDistributionsH1)(x1, x2, x3);
+                    h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
                 }
             }
         }
@@ -3446,7 +3446,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::findNeighbors(CbArray3D<
 
     SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-    phi[REST] = (*ph)(x1, x2, x3);
+    phi[DIR_000] = (*ph)(x1, x2, x3);
 
 
     for (int k = FSTARTDIR; k <= FENDDIR; k++) {
@@ -3466,7 +3466,7 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::findNeighbors2(CbArray3D
 
 	SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
 
-	phi2[REST] = (*ph)(x1, x2, x3);
+	phi2[DIR_000] = (*ph)(x1, x2, x3);
 
 
 	for (int k = FSTARTDIR; k <= FENDDIR; k++) {
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
index 7ee35063f0cfeb9379313a38b9eeb6f0e6388d49..fee069a10092ba2ba3b7dd5e3ff666a00f9dfa2d 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
@@ -469,7 +469,7 @@ void RheologyInterpolationProcessor::calcInterpolatedNode(LBMReal* f, /*LBMReal
    f[BSW]  = f_TNE  + xs*x_TNE  + ys*y_TNE  + zs*z_TNE  + xs*ys*xy_TNE  + xs*zs*xz_TNE  + ys*zs*yz_TNE  + feq[BSW];
    f[BSE]  = f_TNW  + xs*x_TNW  + ys*y_TNW  + zs*z_TNW  + xs*ys*xy_TNW  + xs*zs*xz_TNW  + ys*zs*yz_TNW  + feq[BSE];
    f[BNW]  = f_TSE  + xs*x_TSE  + ys*y_TSE  + zs*z_TSE  + xs*ys*xy_TSE  + xs*zs*xz_TSE  + ys*zs*yz_TSE  + feq[BNW];
-   f[REST] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[REST];
+   f[DIR_000] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO                                                 + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 //Position SWB -0.25, -0.25, -0.25
@@ -658,7 +658,7 @@ void RheologyInterpolationProcessor::calcInterpolatedNodeFC(LBMReal* f, LBMReal
    f[BNW]  = f_TSE  + feq[BNW];
    f[BSE]  = f_TNW  + feq[BSE];
    f[BSW]  = f_TNE  + feq[BSW];
-   f[REST] = f_ZERO + feq[REST];
+   f[DIR_000] = f_ZERO + feq[DIR_000];
 }
 //////////////////////////////////////////////////////////////////////////
 void RheologyInterpolationProcessor::calcInterpolatedVelocity(LBMReal x, LBMReal y, LBMReal z, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3)
diff --git a/src/cpu/VirtualFluidsCore/Visitors/BoundaryConditionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/BoundaryConditionsBlockVisitor.cpp
index a6372fc31712899dab0b8edaf919a141663991ca..f5c87b9fc695d81ad492f89113f2d9e5c56fa9a7 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/BoundaryConditionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/BoundaryConditionsBlockVisitor.cpp
@@ -39,7 +39,7 @@
 #include "D3Q27EsoTwist3DSplittedVector.h"
 #include "DataSet3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "BCAdapter.h"
 #include "Block3D.h"
 #include "BCArray3D.h"
@@ -52,7 +52,9 @@
 #include "ThixotropyVelocityWithDensityBCAlgorithm.h"
 
 
-BoundaryConditionsBlockVisitor::BoundaryConditionsBlockVisitor() : Block3DVisitor(0, Grid3DSystem::MAXLEVEL) {}
+BoundaryConditionsBlockVisitor::BoundaryConditionsBlockVisitor() : Block3DVisitor(0, D3Q27System::MAXLEVEL)
+{
+}
 //////////////////////////////////////////////////////////////////////////
 BoundaryConditionsBlockVisitor::~BoundaryConditionsBlockVisitor() = default;
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/ChangeBoundaryDensityBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/ChangeBoundaryDensityBlockVisitor.cpp
index b222d2c67af03cae39a385c91f25b29962565c39..e26b59729594fc3175e523e25d23ce7adc56d74e 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/ChangeBoundaryDensityBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/ChangeBoundaryDensityBlockVisitor.cpp
@@ -4,11 +4,11 @@
 #include "Block3D.h"
 #include "BoundaryConditions.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 
 ChangeBoundaryDensityBlockVisitor::ChangeBoundaryDensityBlockVisitor(float oldBoundaryDensity, float newBoundaryDensity)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), oldBoundaryDensity(oldBoundaryDensity),
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), oldBoundaryDensity(oldBoundaryDensity),
       newBoundaryDensity(newBoundaryDensity)
 {
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/CheckRatioBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/CheckRatioBlockVisitor.cpp
index 1f8a9c30e1f7a98812a06716cb461ee0e7d41aba..d329763a43d6985b8930ec0e73b7a06b991801d0 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/CheckRatioBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/CheckRatioBlockVisitor.cpp
@@ -1,10 +1,10 @@
 #include "CheckRatioBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 
 CheckRatioBlockVisitor::CheckRatioBlockVisitor(int levelDepth /*shut be maxGridLevel*/, bool includeNotActiveBlocks)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), levelDepth(levelDepth), includeNotActiveBlocks(includeNotActiveBlocks),
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), levelDepth(levelDepth), includeNotActiveBlocks(includeNotActiveBlocks),
       state(true)
 {
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/GenBlocksGridVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/GenBlocksGridVisitor.cpp
index a8270d40d7b2e193024056551591f6b0f3464b5e..29ea3bfda98c2ce191d1f7c5bc20691049dc2a04 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/GenBlocksGridVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/GenBlocksGridVisitor.cpp
@@ -35,7 +35,7 @@
 #include "Block3D.h"
 #include "CoordinateTransformation3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 
 #include <geometry3d/GbObject3D.h>
 
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index aba67e6749a47e1c06bf28f01d284799eb39c328..fd0b68655a1134fe2dc1035279f9da143c088ae3 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -38,10 +38,10 @@
 #include "DataSet3D.h"
 #include "EsoTwist3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 
-InitDistributionsBlockVisitor::InitDistributionsBlockVisitor() : Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+InitDistributionsBlockVisitor::InitDistributionsBlockVisitor() : Block3DVisitor(0, D3Q27System::MAXLEVEL)
 {
     this->setVx1(0.0);
     this->setVx2(0.0);
@@ -273,7 +273,7 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
                f[BNW]  = f_TSE  + feq[BNW];
                f[BSE]  = f_TNW  + feq[BSE];
                f[BSW]  = f_TNE  + feq[BSW];
-               f[REST] = f_ZERO + feq[REST];
+               f[DIR_000] = f_ZERO + feq[DIR_000];
 
                //calcFeqsFct(f,rho,vx1,vx2,vx3);
                distributions->setDistribution(f, ix1, ix2, ix3);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
index b10151a9d2926546c2807db05912f79f6815bf86..c20cccdf64a80a388889d23463b62d3b25c03f77 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsFromFileBlockVisitor.cpp
@@ -5,14 +5,14 @@
 #include "DataSet3D.h"
 #include "EsoTwist3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "InitDensityLBMKernel.h"
 #include "LBMKernel.h"
 #include <basics/utilities/UbFileInputASCII.h>
 
 InitDistributionsFromFileBlockVisitor::InitDistributionsFromFileBlockVisitor(/*LBMReal nu, */ LBMReal rho,
                                                                              std::string filename)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), /*nu(nu),*/ rho(rho)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), /*nu(nu),*/ rho(rho)
 {
     UbFileInputASCII in(filename);
     if (!in) {
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
index 343353c7b7caf58269e1b25799dbbc730fc227b5..567ce2e7ff5b40f3c8042bd404394a3fbf9ffee4 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
@@ -7,7 +7,7 @@
 #include "D3Q27EsoTwist3DSplittedVector.h"
 #include "DataSet3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "InterpolationProcessor.h"
 #include "LBMKernel.h"
 #include <CbArray2D.h>
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
index f9be891771838940294a5e3b348a0b4d8c31f366..d35b9d18fe2539affc528e27a41abded0fddb55d 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
@@ -34,7 +34,7 @@
 #include "InitThixotropyBlockVisitor.h"
 #include "LBMKernel.h"
 #include "BCProcessor.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "DataSet3D.h"
 #include "EsoTwist3D.h"
 #include "Grid3D.h"
@@ -42,7 +42,7 @@
 #include "BCArray3D.h"
 
 InitThixotropyBlockVisitor::InitThixotropyBlockVisitor()
-   : Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+   : Block3DVisitor(0, D3Q27System::MAXLEVEL)
 {
    //this->setVx1(0.0);
    //this->setVx2(0.0);
@@ -56,7 +56,7 @@ InitThixotropyBlockVisitor::InitThixotropyBlockVisitor()
 }
 //////////////////////////////////////////////////////////////////////////
 //InitThixotropyBlockVisitor::InitThixotropyBlockVisitor(LBMReal lambda /*LBMReal nu, LBMReal D, LBMReal rho, LBMReal vx1, LBMReal vx2, LBMReal vx3, LBMReal c, LBMReal f1, LBMReal f2, LBMReal f3*/)
-//	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+//	: Block3DVisitor(0, D3Q27System::MAXLEVEL)
 //{
 //	//this->setVx1(vx1);
 //	//this->setVx2(vx2);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseBoundaryConditionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseBoundaryConditionsBlockVisitor.cpp
index 34dba741103d8160304507b5e68a311c5149ddbd..003d5d31204fafc82f78a0fddb04897c2c60e77f 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseBoundaryConditionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseBoundaryConditionsBlockVisitor.cpp
@@ -39,14 +39,14 @@
 #include "D3Q27EsoTwist3DSplittedVector.h"
 #include "DataSet3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "BCAdapter.h"
 #include "Block3D.h"
 #include "BCArray3D.h"
 #include "LBMKernel.h"
 
 MultiphaseBoundaryConditionsBlockVisitor::MultiphaseBoundaryConditionsBlockVisitor() :
-Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+Block3DVisitor(0, D3Q27System::MAXLEVEL)
 {
 
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
index 484972355fa4685769b092721fc566d92ecbe71e..a6ae293a60ead3d34ba83bced8dc01cb7f25a21f 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
@@ -38,11 +38,11 @@
 #include "DataSet3D.h"
 #include "EsoTwist3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 
 MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor() 
-	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+	: Block3DVisitor(0, D3Q27System::MAXLEVEL)
 {
 	this->setVx1(0.0);
 	this->setVx2(0.0);
@@ -51,7 +51,7 @@ MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor
 }
 //////////////////////////////////////////////////////////////////////////
 MultiphaseInitDistributionsBlockVisitor::MultiphaseInitDistributionsBlockVisitor( LBMReal densityRatio, LBMReal vx1, LBMReal vx2, LBMReal vx3, LBMReal rho)
-	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), densityRatio(densityRatio) 
+	: Block3DVisitor(0, D3Q27System::MAXLEVEL), densityRatio(densityRatio) 
 {
 	this->setVx1(vx1);
 	this->setVx2(vx2);
@@ -253,7 +253,7 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 					f[BNW]  =  geq[BNW]  ;
 					f[BSE]  =  geq[BSE]  ;
 					f[BSW]  =  geq[BSW]  ;
-					f[REST] =  geq[REST] ;
+					f[DIR_000] =  geq[DIR_000] ;
 
 					distributionsF->setDistribution(f, ix1, ix2, ix3);
 					distributionsF->setDistributionInv(f, ix1, ix2, ix3);
@@ -284,7 +284,7 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 					f[BNW]  =  phi * feq[BNW]  / rho;
 					f[BSE]  =  phi * feq[BSE]  / rho;
 					f[BSW]  =  phi * feq[BSW]  / rho;
-					f[REST] =  phi * feq[REST] / rho;
+					f[DIR_000] =  phi * feq[DIR_000] / rho;
 
 					distributionsH->setDistribution(f, ix1, ix2, ix3);
 					distributionsH->setDistributionInv(f, ix1, ix2, ix3);
@@ -317,7 +317,7 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 						f[BNW]  = (1.-phi) * feq[BNW] / rho;
 						f[BSE]  = (1.-phi) * feq[BSE] / rho;
 						f[BSW]  = (1.-phi) * feq[BSW] / rho;
-						f[REST] = (1.-phi) * feq[REST] / rho;
+						f[DIR_000] = (1.-phi) * feq[DIR_000] / rho;
 
                         distributionsH2->setDistribution(f, ix1, ix2, ix3);
                         distributionsH2->setDistributionInv(f, ix1, ix2, ix3);                    
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
index afe288e8c70a7ab3b0ceadae721b68db6a0102f4..4990690e2d7d464cfbdc69f2966655568021e7d0 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
@@ -1,17 +1,17 @@
 #include "MultiphaseSetKernelBlockVisitor.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMSystem.h"
 #include "Block3D.h"
 #include "Grid3D.h"
 
 //SetKernelBlockVisitor::SetKernelBlockVisitor(LBMKernel3DPtr kernel, LBMReal nue) : 
-//                        Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nue(nue)
+//                        Block3DVisitor(0, D3Q27System::MAXLEVEL), kernel(kernel), nue(nue)
 //{
 //
 //}
 //////////////////////////////////////////////////////////////////////////
 //SetKernelBlockVisitor::SetKernelBlockVisitor( LBMKernel3DPtr kernel, LBMReal nue, double availMem, double needMem ) : 
-//                                              Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nue(nue)
+//                                              Block3DVisitor(0, D3Q27System::MAXLEVEL), kernel(kernel), nue(nue)
 //{
 //   if (needMem > availMem)
 //   {
@@ -20,7 +20,7 @@
 //}
 //////////////////////////////////////////////////////////////////////////
 MultiphaseSetKernelBlockVisitor::MultiphaseSetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nuL, LBMReal nuG, double availMem, double needMem, MultiphaseSetKernelBlockVisitor::Action action /*= SetKernelBlockVisitor::New*/) :
-	Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nuL(nuL), nuG(nuG), action(action), dataSetFlag(true)
+	Block3DVisitor(0, D3Q27System::MAXLEVEL), kernel(kernel), nuL(nuL), nuG(nuG), action(action), dataSetFlag(true)
 {
 	if (needMem > availMem)
 	{
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
index 8a201fa3e7236572ca94c83c72ee90149d676390..e10a9774d85bd1e6d32c2d7b320d3a4b591da94d 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
@@ -38,11 +38,11 @@
 #include "DataSet3D.h"
 #include "EsoTwist3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 
 MultiphaseVelocityFormInitDistributionsBlockVisitor::MultiphaseVelocityFormInitDistributionsBlockVisitor() 
-	: Block3DVisitor(0, Grid3DSystem::MAXLEVEL)
+	: Block3DVisitor(0, D3Q27System::MAXLEVEL)
 {
 	this->setVx1(0.0);
 	this->setVx2(0.0);
@@ -245,7 +245,7 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 					f[BNW]  =  geq[BNW]  ;
 					f[BSE]  =  geq[BSE]  ;
 					f[BSW]  =  geq[BSW]  ;
-					f[REST] =  geq[REST] ;
+					f[DIR_000] =  geq[DIR_000] ;
 
 					distributionsF->setDistribution(f, ix1, ix2, ix3);
 					distributionsF->setDistributionInv(f, ix1, ix2, ix3);
@@ -276,7 +276,7 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 					f[BNW]  =  phi * feq[BNW]  ;// / rho;
 					f[BSE]  =  phi * feq[BSE]  ;// / rho;
 					f[BSW]  =  phi * feq[BSW]  ;// / rho;
-					f[REST] =  phi * feq[REST] ;// / rho;
+					f[DIR_000] =  phi * feq[DIR_000] ;// / rho;
 
 					distributionsH->setDistribution(f, ix1, ix2, ix3);
 					distributionsH->setDistributionInv(f, ix1, ix2, ix3);
@@ -309,7 +309,7 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 						f[BNW]  = (1.-phi) * feq[BNW] ;// / rho;
 						f[BSE]  = (1.-phi) * feq[BSE] ;// / rho;
 						f[BSW]  = (1.-phi) * feq[BSW] ;// / rho;
-						f[REST] = (1.-phi) * feq[REST];//  / rho;
+						f[DIR_000] = (1.-phi) * feq[DIR_000];//  / rho;
 
                         distributionsH2->setDistribution(f, ix1, ix2, ix3);
                         distributionsH2->setDistributionInv(f, ix1, ix2, ix3);                    
diff --git a/src/cpu/VirtualFluidsCore/Visitors/OverlapBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/OverlapBlockVisitor.cpp
index 338c6de90499db39836910c47a7b60f00ee7c675..54bbeda59663bf16173abcdb302f123a97b776ac 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/OverlapBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/OverlapBlockVisitor.cpp
@@ -1,10 +1,10 @@
 #include "OverlapBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 
 OverlapBlockVisitor::OverlapBlockVisitor(int levelDepth /*shut be maxGridLevel*/, bool includeNotActiveBlocks)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), levelDepth(levelDepth), includeNotActiveBlocks(includeNotActiveBlocks)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), levelDepth(levelDepth), includeNotActiveBlocks(includeNotActiveBlocks)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/RatioBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/RatioBlockVisitor.cpp
index bf25b8876f540c6ff5f678b32ad67ddd736c145d..137c737e6dd85853da328b0d0cc9ccf1bca1f878 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/RatioBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/RatioBlockVisitor.cpp
@@ -1,10 +1,10 @@
 #include "RatioBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 
 RatioBlockVisitor::RatioBlockVisitor(int levelDepth, bool includeNotActiveBlocks)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), maxLevelRatio(1), expandBlocks(true), levelDepth(levelDepth),
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), maxLevelRatio(1), expandBlocks(true), levelDepth(levelDepth),
       includeNotActiveBlocks(includeNotActiveBlocks)
 {
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/RatioSmoothBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/RatioSmoothBlockVisitor.cpp
index 3bcdd29299a8183a2cda60253101df89b681a51a..b48aee7cd8dcd3d3cc33c949e8a71281214c8f19 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/RatioSmoothBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/RatioSmoothBlockVisitor.cpp
@@ -1,10 +1,10 @@
 #include "RatioSmoothBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 
 RatioSmoothBlockVisitor::RatioSmoothBlockVisitor(int levelDepth, bool includeNotActiveBlocks)
-    : Block3DVisitor(Grid3DSystem::MAXLEVEL, 0), maxLevelRatio(1), expandBlocks(true), levelDepth(levelDepth),
+    : Block3DVisitor(D3Q27System::MAXLEVEL, 0), maxLevelRatio(1), expandBlocks(true), levelDepth(levelDepth),
       includeNotActiveBlocks(includeNotActiveBlocks)
 {
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/RenumberBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/RenumberBlockVisitor.cpp
index 538fd95118b0c07069ba854c2fbb8264713907d3..b2a4d9337ad49edc4e508555929a6ca64e6f9243 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/RenumberBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/RenumberBlockVisitor.cpp
@@ -1,12 +1,12 @@
 #include "RenumberBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMSystem.h"
 
 int RenumberBlockVisitor::counter = 0;
 
-RenumberBlockVisitor::RenumberBlockVisitor() : Block3DVisitor(0, Grid3DSystem::MAXLEVEL) {}
+RenumberBlockVisitor::RenumberBlockVisitor() : Block3DVisitor(0, D3Q27System::MAXLEVEL) {}
 //////////////////////////////////////////////////////////////////////////
 void RenumberBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
 {
diff --git a/src/cpu/VirtualFluidsCore/Visitors/RenumberGridVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/RenumberGridVisitor.cpp
index fc9c5c203c5d631ae7e125f75d72d70e8502890d..ed9a3ee59c87ab755416eecd5468a4cc763837e4 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/RenumberGridVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/RenumberGridVisitor.cpp
@@ -1,7 +1,7 @@
 #include "RenumberGridVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 //#include <mpi.h>
 
 RenumberGridVisitor::RenumberGridVisitor(std::shared_ptr<vf::mpi::Communicator> com) : comm(com) {}
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp
index bbccaa785583fdc810865337af46fca8a9872a65..de3924453dadb396a5eef4a8fb23e85f850d760f 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetBcBlocksBlockVisitor.cpp
@@ -35,11 +35,11 @@
 
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "Interactor3D.h"
 
 SetBcBlocksBlockVisitor::SetBcBlocksBlockVisitor(SPtr<Interactor3D> interactor)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(interactor)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), interactor(interactor)
 {
 }
 
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
index 91f18e7292e8ef28cac159fb45aae485529e3f79..c47c275cf50d9a7789f7c6fd302f1c2b2177f5f6 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
@@ -38,7 +38,7 @@
 
 #include "Block3DVisitor.h"
 #include "D3Q27System.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "Grid3D.h"
 #include "CreateTransmittersHelper.h"
 #include <mpi/Communicator.h>
@@ -75,7 +75,7 @@ protected:
 
 template <class T1, class T2>
 SetConnectorsBlockVisitor<T1, T2>::SetConnectorsBlockVisitor(std::shared_ptr<vf::mpi::Communicator> comm)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), comm(comm)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), comm(comm)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetForcingBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetForcingBlockVisitor.cpp
index 679b63de44dd451f030aaf866bc259579efab8ef..abf828a06e0ec83b492ff9107be4a9a3c4445674 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetForcingBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetForcingBlockVisitor.cpp
@@ -1,18 +1,18 @@
 #include "SetForcingBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMSystem.h"
 
 SetForcingBlockVisitor::SetForcingBlockVisitor(LBMReal forcingX1, LBMReal forcingX2, LBMReal forcingX3)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), forcingX1(forcingX1), forcingX2(forcingX2), forcingX3(forcingX3)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), forcingX1(forcingX1), forcingX2(forcingX2), forcingX3(forcingX3)
 {
     ftype = 0;
 }
 //////////////////////////////////////////////////////////////////////////
 SetForcingBlockVisitor::SetForcingBlockVisitor(const mu::Parser &muForcingX1, const mu::Parser &muForcingX2,
                                                const mu::Parser &muForcingX3)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), muForcingX1(muForcingX1), muForcingX2(muForcingX2),
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), muForcingX1(muForcingX1), muForcingX2(muForcingX2),
       muForcingX3(muForcingX3)
 
 {
@@ -21,7 +21,7 @@ SetForcingBlockVisitor::SetForcingBlockVisitor(const mu::Parser &muForcingX1, co
 //////////////////////////////////////////////////////////////////////////
 SetForcingBlockVisitor::SetForcingBlockVisitor(const std::string &sForcingX1, const std::string &sForcingX2,
                                                const std::string &sForcingX3)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), sForcingX1(sForcingX1), sForcingX2(sForcingX2), sForcingX3(sForcingX3)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), sForcingX1(sForcingX1), sForcingX2(sForcingX2), sForcingX3(sForcingX3)
 
 {
     ftype = 2;
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationConnectorsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationConnectorsBlockVisitor.cpp
index 6a55ee5af55df96b4c1335976728ca7e08ee8ece..4badec1f4d1a46afce4bfdb0b0172af80b606961 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationConnectorsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationConnectorsBlockVisitor.cpp
@@ -36,14 +36,14 @@
 #include "FineToCoarseVectorConnector.h"
 #include "TwoDistributionsFullDirectConnector.h"
 #include "TwoDistributionsFullVectorConnector.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include <basics/transmitter/TbTransmitterLocal.h>
 
 #include <mpi/Communicator.h>
 #include "InterpolationProcessor.h"
 
 SetInterpolationConnectorsBlockVisitor::SetInterpolationConnectorsBlockVisitor(std::shared_ptr<vf::mpi::Communicator> comm, LBMReal nue, SPtr<InterpolationProcessor> iProcessor) :
-Block3DVisitor(0, Grid3DSystem::MAXLEVEL), 
+Block3DVisitor(0, D3Q27System::MAXLEVEL), 
 	comm(comm),
 	nue(nue),
 	iProcessor(iProcessor)
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
index bb1ae79620179f65abd51672bd9958f471c398c7..c4b542f1c5277397427bb9f3d0520cb2ba131d3f 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
@@ -1,11 +1,11 @@
 #include "SetInterpolationDirsBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include <D3Q27System.h>
 
 SetInterpolationDirsBlockVisitor::SetInterpolationDirsBlockVisitor(std::vector<int> &dirs)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), dirs(dirs)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), dirs(dirs)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
index 5c813d28951b24269735b0dbf1f84fdfb360cf31..8cf73ca6f80de2fe6521c65cdb4a46634167ad87 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
@@ -38,7 +38,7 @@
 #include "Block3D.h"
 #include "DataSet3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMKernel.h"
 #include "LBMSystem.h"
 #include <utility>
@@ -46,7 +46,7 @@
 //////////////////////////////////////////////////////////////////////////
 SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, double availMem, double needMem,
                                              SetKernelBlockVisitor::Action action)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true)
 {
     if (needMem > availMem) {
         throw UbException(UB_EXARGS, "SetKernelBlockVisitor: Not enough memory!!!");
@@ -55,7 +55,7 @@ SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue
 
 SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int &numberOfProcesses,
                                              SetKernelBlockVisitor::Action action)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true),
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true),
       numberOfProcesses(numberOfProcesses)
 {
 }
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp
index 3354755f22f18df700523d795c8fced0d0f19628..e78300c5af6590da72bec3b28516143e4ec18a1d 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetSolidBlocksBlockVisitor.cpp
@@ -37,11 +37,11 @@
 
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "Interactor3D.h"
 
 SetSolidBlocksBlockVisitor::SetSolidBlocksBlockVisitor(SPtr<Interactor3D> interactor)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), interactor(std::move(interactor))
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), interactor(std::move(interactor))
 {
 }
 
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetSpongeLayerBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetSpongeLayerBlockVisitor.cpp
index 040f54b0fd645a206952ce2c992d4abeea3ace85..789b3e3dc932462155a897be1c2c998400adb7f3 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetSpongeLayerBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetSpongeLayerBlockVisitor.cpp
@@ -1,5 +1,5 @@
 #include "SetSpongeLayerBlockVisitor.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMSystem.h"
 
 #include "Block3D.h"
@@ -7,7 +7,7 @@
 #include "LBMKernel.h"
 
 SetSpongeLayerBlockVisitor::SetSpongeLayerBlockVisitor(const mu::Parser &spongeLayer)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), spongeLayer(spongeLayer)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), spongeLayer(spongeLayer)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetUndefinedNodesBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetUndefinedNodesBlockVisitor.cpp
index 15c8b82bae5c95f6783f73c7e7a70004f3c76574..6e3e6823de9bd73362c273cae32c09cd61a5220c 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetUndefinedNodesBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetUndefinedNodesBlockVisitor.cpp
@@ -5,11 +5,11 @@
 #include "BoundaryConditions.h"
 #include "D3Q27System.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "ILBMKernel.h"
 
 SetUndefinedNodesBlockVisitor::SetUndefinedNodesBlockVisitor(bool twoTypeOfConnectorsCheck)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), twoTypeOfConnectorsCheck(twoTypeOfConnectorsCheck)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), twoTypeOfConnectorsCheck(twoTypeOfConnectorsCheck)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
index c9bbc78cdd575f115d86d3e470cefe50dd636ba6..d763fe680b3ddbec7a20e1fa248b5035e1f94a4b 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
@@ -1,5 +1,5 @@
 #include "SpongeLayerBlockVisitor.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "LBMSystem.h"
 
 #include "BCArray3D.h"
@@ -15,7 +15,7 @@ using namespace std;
 
 SpongeLayerBlockVisitor::SpongeLayerBlockVisitor(SPtr<GbCuboid3D> boundingBox, SPtr<LBMKernel> kernel, double nue,
                                                  int dir)
-    : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), boundingBox(boundingBox), kernel(kernel), nue(nue), dir(dir)
+    : Block3DVisitor(0, D3Q27System::MAXLEVEL), boundingBox(boundingBox), kernel(kernel), nue(nue), dir(dir)
 {
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Visitors/ViscosityBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/ViscosityBlockVisitor.cpp
index 67d185d6ac401909d85b99b74ef1ede2d0054a6a..311a8bf19786198e85b00eb500f6e7c90d2d5106 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/ViscosityBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/ViscosityBlockVisitor.cpp
@@ -1,11 +1,11 @@
 #include "ViscosityBlockVisitor.h"
 #include "Block3D.h"
 #include "Grid3D.h"
-#include "Grid3DSystem.h"
+#include "D3Q27System.h"
 #include "ILBMKernel.h"
 #include "LBMSystem.h"
 
-ViscosityBlockVisitor::ViscosityBlockVisitor(LBMReal nu) : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), nu(nu) {}
+ViscosityBlockVisitor::ViscosityBlockVisitor(LBMReal nu) : Block3DVisitor(0, D3Q27System::MAXLEVEL), nu(nu) {}
 //////////////////////////////////////////////////////////////////////////
 void ViscosityBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
 {