From 7031961c5d8a6bf7b53b36a03df8d37198df5e2c Mon Sep 17 00:00:00 2001
From: alena <akaranchuk@list.ru>
Date: Thu, 18 Aug 2022 16:33:24 +0200
Subject: [PATCH] Changing N, S, E, W, T, B,... constans to DIR_00P,
 DIR_00M,... values.

---
 ...tiphaseNonReflectingOutflowBCAlgorithm.cpp |  30 ++---
 .../NonReflectingOutflowBCAlgorithm.cpp       |  70 ++++++------
 ...xotropyNonReflectingOutflowBCAlgorithm.cpp | 104 +++++++++---------
 .../CoProcessors/ShearStressCoProcessor.cpp   |  14 +--
 .../TurbulenceIntensityCoProcessor.cpp        |   6 +-
 .../WriteMultiphaseQuantitiesCoProcessor.cpp  |  60 +++++-----
 ...mpressibleOffsetInterpolationProcessor.cpp |  30 ++---
 .../LBM/InitDensityLBMKernel.cpp              |  42 +++----
 .../LBM/LBMKernelETD3Q27BGK.cpp               |  60 +++++-----
 .../LBM/MultiphaseCumulantLBMKernel.cpp       |  48 ++++----
 ...PressureFilterCompressibleAirLBMKernel.cpp |  46 ++++----
 .../LBM/MultiphasePressureFilterLBMKernel.cpp |  34 +++---
 .../MultiphaseScratchCumulantLBMKernel.cpp    |  58 +++++-----
 ...tiphaseTwoPhaseFieldsCumulantLBMKernel.cpp |  70 ++++++------
 ...eTwoPhaseFieldsPressureFilterLBMKernel.cpp |  46 ++++----
 ...woPhaseFieldsVelocityCumulantLBMKernel.cpp |  46 ++++----
 .../LBM/RheologyInterpolationProcessor.cpp    |  30 ++---
 .../InitDistributionsBlockVisitor.cpp         |  12 +-
 .../Visitors/MetisPartitioningGridVisitor.cpp |   4 +-
 ...ultiphaseInitDistributionsBlockVisitor.cpp |  36 +++---
 ...ocityFormInitDistributionsBlockVisitor.cpp |  36 +++---
 .../SetInterpolationDirsBlockVisitor.cpp      |  42 +++----
 22 files changed, 462 insertions(+), 462 deletions(-)

diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNonReflectingOutflowBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNonReflectingOutflowBCAlgorithm.cpp
index bc9a73dee..fc6019244 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNonReflectingOutflowBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNonReflectingOutflowBCAlgorithm.cpp
@@ -82,7 +82,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
    else if (bcPtr->hasDensityBoundaryFlag(DIR_0P0)) { nx2 += 1; direction = DIR_0P0; }
    else if (bcPtr->hasDensityBoundaryFlag(DIR_0M0)) { nx2 -= 1; direction = DIR_0M0; }
    else if (bcPtr->hasDensityBoundaryFlag(DIR_00P)) { nx3 += 1; direction = DIR_00P; }
-   else if (bcPtr->hasDensityBoundaryFlag(B)) { nx3 -= 1; direction = B; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_00M)) { nx3 -= 1; direction = DIR_00M; }
    else UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on density boundary..."));
 
    distributions->getDistribution(f, x1, x2, x3);
@@ -182,8 +182,8 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  distributionsH->setDistributionInvForDirection(h[DIR_MMM], x1+DX1[DIR_PPP], x2+DX2[DIR_PPP], x3+DX3[DIR_PPP], DIR_PPP);
 
       break;
-   case N:
-      f[N]   = ftemp[N]   * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[N]   ;
+   case DIR_0P0:
+      f[DIR_0P0]   = ftemp[DIR_0P0]   * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_0P0]   ;
       f[DIR_PP0]  = ftemp[DIR_PP0]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_PP0]  ;
       f[DIR_MP0]  = ftemp[DIR_MP0]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_MP0]  ;
       f[DIR_0PP]  = ftemp[DIR_0PP]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_0PP]  ;
@@ -193,7 +193,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_PPM] ;
       f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*f[DIR_MPM] ;
 
-      distributions->setDistributionInvForDirection(f[N],   x1+DX1[DIR_0M0],   x2+DX2[DIR_0M0],   x3+DX3[DIR_0M0],     DIR_0M0);
+      distributions->setDistributionInvForDirection(f[DIR_0P0],   x1+DX1[DIR_0M0],   x2+DX2[DIR_0M0],   x3+DX3[DIR_0M0],     DIR_0M0);
       distributions->setDistributionInvForDirection(f[DIR_PP0],  x1+DX1[DIR_MM0],  x2+DX2[DIR_MM0],  x3+DX3[DIR_MM0],   DIR_MM0);
       distributions->setDistributionInvForDirection(f[DIR_MP0],  x1+DX1[DIR_PM0],  x2+DX2[DIR_PM0],  x3+DX3[DIR_PM0],   DIR_PM0);
       distributions->setDistributionInvForDirection(f[DIR_0PP],  x1+DX1[DIR_0MM],  x2+DX2[DIR_0MM],  x3+DX3[DIR_0MM],   DIR_0MM);
@@ -203,7 +203,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PPM], x1+DX1[DIR_MMP], x2+DX2[DIR_MMP], x3+DX3[DIR_MMP], DIR_MMP);
       distributions->setDistributionInvForDirection(f[DIR_MPM], x1+DX1[DIR_PMP], x2+DX2[DIR_PMP], x3+DX3[DIR_PMP], DIR_PMP);
 
-	  h[N]   = htemp[N]   * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[N]   ;
+	  h[DIR_0P0]   = htemp[DIR_0P0]   * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_0P0]   ;
 	  h[DIR_PP0]  = htemp[DIR_PP0]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_PP0]  ;
 	  h[DIR_MP0]  = htemp[DIR_MP0]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_MP0]  ;
 	  h[DIR_0PP]  = htemp[DIR_0PP]  * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_0PP]  ;
@@ -213,7 +213,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  h[DIR_PPM] = htemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_PPM] ;
 	  h[DIR_MPM] = htemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2)*h[DIR_MPM] ;
 
-	  distributionsH->setDistributionInvForDirection(h[N],   x1+DX1[DIR_0M0],   x2+DX2[DIR_0M0],   x3+DX3[DIR_0M0],     DIR_0M0);
+	  distributionsH->setDistributionInvForDirection(h[DIR_0P0],   x1+DX1[DIR_0M0],   x2+DX2[DIR_0M0],   x3+DX3[DIR_0M0],     DIR_0M0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_PP0],  x1+DX1[DIR_MM0],  x2+DX2[DIR_MM0],  x3+DX3[DIR_MM0],   DIR_MM0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_MP0],  x1+DX1[DIR_PM0],  x2+DX2[DIR_PM0],  x3+DX3[DIR_PM0],   DIR_PM0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_0PP],  x1+DX1[DIR_0MM],  x2+DX2[DIR_0MM],  x3+DX3[DIR_0MM],   DIR_0MM);
@@ -235,7 +235,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2)*f[DIR_PMM] ;
       f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2)*f[DIR_MMM] ;
 
-      distributions->setDistributionInvForDirection(f[DIR_0M0],   x1+DX1[N],   x2+DX2[N],   x3+DX3[N],     N);
+      distributions->setDistributionInvForDirection(f[DIR_0M0],   x1+DX1[DIR_0P0],   x2+DX2[DIR_0P0],   x3+DX3[DIR_0P0], DIR_0P0);
       distributions->setDistributionInvForDirection(f[DIR_PM0],  x1+DX1[DIR_MP0],  x2+DX2[DIR_MP0],  x3+DX3[DIR_MP0],   DIR_MP0);
       distributions->setDistributionInvForDirection(f[DIR_MM0],  x1+DX1[DIR_PP0],  x2+DX2[DIR_PP0],  x3+DX3[DIR_PP0],   DIR_PP0);
       distributions->setDistributionInvForDirection(f[DIR_0MP],  x1+DX1[DIR_0PM],  x2+DX2[DIR_0PM],  x3+DX3[DIR_0PM],   DIR_0PM);
@@ -255,7 +255,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  h[DIR_PMM] = htemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2)*h[DIR_PMM] ;
 	  h[DIR_MMM] = htemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2)*h[DIR_MMM] ;
 
-	  distributionsH->setDistributionInvForDirection(h[DIR_0M0],   x1+DX1[N],   x2+DX2[N],   x3+DX3[N],     N);
+	  distributionsH->setDistributionInvForDirection(h[DIR_0M0],   x1+DX1[DIR_0P0],   x2+DX2[DIR_0P0],   x3+DX3[DIR_0P0], DIR_0P0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_PM0],  x1+DX1[DIR_MP0],  x2+DX2[DIR_MP0],  x3+DX3[DIR_MP0],   DIR_MP0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_MM0],  x1+DX1[DIR_PP0],  x2+DX2[DIR_PP0],  x3+DX3[DIR_PP0],   DIR_PP0);
 	  distributionsH->setDistributionInvForDirection(h[DIR_0MP],  x1+DX1[DIR_0PM],  x2+DX2[DIR_0PM],  x3+DX3[DIR_0PM],   DIR_0PM);
@@ -277,7 +277,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3)*f[DIR_PMP] ;
       f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3)*f[DIR_MMP] ;
 
-      distributions->setDistributionInvForDirection(f[DIR_00P],   x1+DX1[B],   x2+DX2[B],   x3+DX3[B],     B);
+      distributions->setDistributionInvForDirection(f[DIR_00P],   x1+DX1[DIR_00M],   x2+DX2[DIR_00M],   x3+DX3[DIR_00M], DIR_00M);
       distributions->setDistributionInvForDirection(f[DIR_P0P],  x1+DX1[DIR_M0M],  x2+DX2[DIR_M0M],  x3+DX3[DIR_M0M],   DIR_M0M);
       distributions->setDistributionInvForDirection(f[DIR_M0P],  x1+DX1[DIR_P0M],  x2+DX2[DIR_P0M],  x3+DX3[DIR_P0M],   DIR_P0M);
       distributions->setDistributionInvForDirection(f[DIR_0PP],  x1+DX1[DIR_0MM],  x2+DX2[DIR_0MM],  x3+DX3[DIR_0MM],   DIR_0MM);
@@ -297,7 +297,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  h[DIR_PMP] = htemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3)*h[DIR_PMP] ;
 	  h[DIR_MMP] = htemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3)*h[DIR_MMP] ;
 
-	  distributionsH->setDistributionInvForDirection(h[DIR_00P],   x1+DX1[B],   x2+DX2[B],   x3+DX3[B],     B);
+	  distributionsH->setDistributionInvForDirection(h[DIR_00P],   x1+DX1[DIR_00M],   x2+DX2[DIR_00M],   x3+DX3[DIR_00M], DIR_00M);
 	  distributionsH->setDistributionInvForDirection(h[DIR_P0P],  x1+DX1[DIR_M0M],  x2+DX2[DIR_M0M],  x3+DX3[DIR_M0M],   DIR_M0M);
 	  distributionsH->setDistributionInvForDirection(h[DIR_M0P],  x1+DX1[DIR_P0M],  x2+DX2[DIR_P0M],  x3+DX3[DIR_P0M],   DIR_P0M);
 	  distributionsH->setDistributionInvForDirection(h[DIR_0PP],  x1+DX1[DIR_0MM],  x2+DX2[DIR_0MM],  x3+DX3[DIR_0MM],   DIR_0MM);
@@ -308,8 +308,8 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  distributionsH->setDistributionInvForDirection(h[DIR_MMP], x1+DX1[DIR_PPM], x2+DX2[DIR_PPM], x3+DX3[DIR_PPM], DIR_PPM);
 
       break;
-   case B:
-      f[B]   = ftemp[B]   * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[B]   ;
+   case DIR_00M:
+      f[DIR_00M]   = ftemp[DIR_00M]   * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_00M]   ;
       f[DIR_P0M]  = ftemp[DIR_P0M]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_P0M]  ;
       f[DIR_M0M]  = ftemp[DIR_M0M]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_M0M]  ;
       f[DIR_0PM]  = ftemp[DIR_0PM]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_0PM]  ;
@@ -319,7 +319,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_PMM] ;
       f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*f[DIR_MMM] ;
 
-      distributions->setDistributionInvForDirection(f[B],   x1+DX1[DIR_00P],   x2+DX2[DIR_00P],   x3+DX3[DIR_00P],     DIR_00P);
+      distributions->setDistributionInvForDirection(f[DIR_00M],   x1+DX1[DIR_00P],   x2+DX2[DIR_00P],   x3+DX3[DIR_00P],     DIR_00P);
       distributions->setDistributionInvForDirection(f[DIR_P0M],  x1+DX1[DIR_M0P],  x2+DX2[DIR_M0P],  x3+DX3[DIR_M0P],   DIR_M0P);
       distributions->setDistributionInvForDirection(f[DIR_M0M],  x1+DX1[DIR_P0P],  x2+DX2[DIR_P0P],  x3+DX3[DIR_P0P],   DIR_P0P);
       distributions->setDistributionInvForDirection(f[DIR_0PM],  x1+DX1[DIR_0MP],  x2+DX2[DIR_0MP],  x3+DX3[DIR_0MP],   DIR_0MP);
@@ -329,7 +329,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PMM], x1+DX1[DIR_MPP], x2+DX2[DIR_MPP], x3+DX3[DIR_MPP], DIR_MPP);
       distributions->setDistributionInvForDirection(f[DIR_MMM], x1+DX1[DIR_PPP], x2+DX2[DIR_PPP], x3+DX3[DIR_PPP], DIR_PPP);
 
-	  h[B]   = htemp[B]   * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[B]   ;
+	  h[DIR_00M]   = htemp[DIR_00M]   * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_00M]   ;
 	  h[DIR_P0M]  = htemp[DIR_P0M]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_P0M]  ;
 	  h[DIR_M0M]  = htemp[DIR_M0M]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_M0M]  ;
 	  h[DIR_0PM]  = htemp[DIR_0PM]  * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_0PM]  ;
@@ -339,7 +339,7 @@ void MultiphaseNonReflectingOutflowBCAlgorithm::applyBC()
 	  h[DIR_PMM] = htemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_PMM] ;
 	  h[DIR_MMM] = htemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3)*h[DIR_MMM] ;
 
-	  distributionsH->setDistributionInvForDirection(h[B],   x1+DX1[DIR_00P],   x2+DX2[DIR_00P],   x3+DX3[DIR_00P],     DIR_00P);
+	  distributionsH->setDistributionInvForDirection(h[DIR_00M],   x1+DX1[DIR_00P],   x2+DX2[DIR_00P],   x3+DX3[DIR_00P],     DIR_00P);
 	  distributionsH->setDistributionInvForDirection(h[DIR_P0M],  x1+DX1[DIR_M0P],  x2+DX2[DIR_M0P],  x3+DX3[DIR_M0P],   DIR_M0P);
 	  distributionsH->setDistributionInvForDirection(h[DIR_M0M],  x1+DX1[DIR_P0P],  x2+DX2[DIR_P0P],  x3+DX3[DIR_P0P],   DIR_P0P);
 	  distributionsH->setDistributionInvForDirection(h[DIR_0PM],  x1+DX1[DIR_0MP],  x2+DX2[DIR_0MP],  x3+DX3[DIR_0MP],   DIR_0MP);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
index f92357cdb..6fa4c7b5d 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
@@ -69,24 +69,24 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
     int direction = -1;
 
     // flag points in direction of fluid
-    if (bcPtr->hasDensityBoundaryFlag(E)) {
+    if (bcPtr->hasDensityBoundaryFlag(DIR_P00)) {
         nx1 += 1;
-        direction = E;
-    } else if (bcPtr->hasDensityBoundaryFlag(W)) {
+        direction = DIR_P00;
+    } else if (bcPtr->hasDensityBoundaryFlag(DIR_M00)) {
         nx1 -= 1;
-        direction = W;
-    } else if (bcPtr->hasDensityBoundaryFlag(N)) {
+        direction = DIR_M00;
+    } else if (bcPtr->hasDensityBoundaryFlag(DIR_0P0)) {
         nx2 += 1;
-        direction = N;
-    } else if (bcPtr->hasDensityBoundaryFlag(S)) {
+        direction = DIR_0P0;
+    } else if (bcPtr->hasDensityBoundaryFlag(DIR_0M0)) {
         nx2 -= 1;
-        direction = S;
-    } else if (bcPtr->hasDensityBoundaryFlag(T)) {
+        direction = DIR_0M0;
+    } else if (bcPtr->hasDensityBoundaryFlag(DIR_00P)) {
         nx3 += 1;
-        direction = T;
-    } else if (bcPtr->hasDensityBoundaryFlag(B)) {
+        direction = DIR_00P;
+    } else if (bcPtr->hasDensityBoundaryFlag(DIR_00M)) {
         nx3 -= 1;
-        direction = B;
+        direction = DIR_00M;
     } else
         UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on density boundary..."));
 
@@ -97,9 +97,9 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
     calcMacrosFct(f, rho, vx1, vx2, vx3);
 
     switch (direction) {
-        case E:
+        case DIR_P00:
             f[DIR_P00]   = ftemp[DIR_P00] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P00];
-            f[NE]  = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+            f[DIR_PP0]  = ftemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PP0];
             f[DIR_PM0]  = ftemp[DIR_PM0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PM0];
             f[DIR_P0P]  = ftemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P0P];
             f[DIR_P0M]  = ftemp[DIR_P0M] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P0M];
@@ -108,8 +108,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_PPM] = ftemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PPM];
             f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PMM];
 
-            distributions->setDistributionInvForDirection(f[DIR_P00], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
-            distributions->setDistributionInvForDirection(f[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+            distributions->setDistributionInvForDirection(f[DIR_P00], x1 + DX1[DIR_M00], x2 + DX2[DIR_M00], x3 + DX3[DIR_M00], DIR_M00);
+            distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
             distributions->setDistributionInvForDirection(f[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
             distributions->setDistributionInvForDirection(f[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
             distributions->setDistributionInvForDirection(f[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
@@ -118,8 +118,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
             distributions->setDistributionInvForDirection(f[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
             break;
-        case W:
-            f[W]   = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+        case DIR_M00:
+            f[DIR_M00]   = ftemp[DIR_M00] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_M00];
             f[DIR_MP0]  = ftemp[DIR_MP0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MP0];
             f[DIR_MM0]  = ftemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MM0];
             f[DIR_M0P]  = ftemp[DIR_M0P] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_M0P];
@@ -129,9 +129,9 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_MPM] = ftemp[DIR_MPM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MPM];
             f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MMM];
 
-            distributions->setDistributionInvForDirection(f[W], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], E);
+            distributions->setDistributionInvForDirection(f[DIR_M00], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], DIR_P00);
             distributions->setDistributionInvForDirection(f[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
-            distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+            distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
             distributions->setDistributionInvForDirection(f[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
             distributions->setDistributionInvForDirection(f[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
             distributions->setDistributionInvForDirection(f[DIR_MPP], x1 + DX1[DIR_PMM], x2 + DX2[DIR_PMM], x3 + DX3[DIR_PMM], DIR_PMM);
@@ -139,9 +139,9 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
             distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
             break;
-        case N:
-            f[N]   = ftemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[N];
-            f[NE]  = ftemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NE];
+        case DIR_0P0:
+            f[DIR_0P0]   = ftemp[DIR_0P0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0P0];
+            f[DIR_PP0]  = ftemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_PP0];
             f[DIR_MP0]  = ftemp[DIR_MP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_MP0];
             f[DIR_0PP]  = ftemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0PP];
             f[DIR_0PM]  = ftemp[DIR_0PM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0PM];
@@ -150,8 +150,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_PPM] = ftemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_PPM];
             f[DIR_MPM] = ftemp[DIR_MPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_MPM];
 
-            distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
-            distributions->setDistributionInvForDirection(f[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+            distributions->setDistributionInvForDirection(f[DIR_0P0], x1 + DX1[DIR_0M0], x2 + DX2[DIR_0M0], x3 + DX3[DIR_0M0], DIR_0M0);
+            distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
             distributions->setDistributionInvForDirection(f[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
             distributions->setDistributionInvForDirection(f[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
             distributions->setDistributionInvForDirection(f[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
@@ -160,8 +160,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
             distributions->setDistributionInvForDirection(f[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
             break;
-        case S:
-            f[S]   = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+        case DIR_0M0:
+            f[DIR_0M0]   = ftemp[DIR_0M0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_0M0];
             f[DIR_PM0]  = ftemp[DIR_PM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_PM0];
             f[DIR_MM0]  = ftemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_MM0];
             f[DIR_0MP]  = ftemp[DIR_0MP] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_0MP];
@@ -171,9 +171,9 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_PMM];
             f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_MMM];
 
-            distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
+            distributions->setDistributionInvForDirection(f[DIR_0M0], x1 + DX1[DIR_0P0], x2 + DX2[DIR_0P0], x3 + DX3[DIR_0P0], DIR_0P0);
             distributions->setDistributionInvForDirection(f[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
-            distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+            distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
             distributions->setDistributionInvForDirection(f[DIR_0MP], x1 + DX1[DIR_0PM], x2 + DX2[DIR_0PM], x3 + DX3[DIR_0PM], DIR_0PM);
             distributions->setDistributionInvForDirection(f[DIR_0MM], x1 + DX1[DIR_0PP], x2 + DX2[DIR_0PP], x3 + DX3[DIR_0PP], DIR_0PP);
             distributions->setDistributionInvForDirection(f[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
@@ -181,8 +181,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
             distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
             break;
-        case T:
-            f[T]   = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+        case DIR_00P:
+            f[DIR_00P]   = ftemp[DIR_00P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_00P];
             f[DIR_P0P]  = ftemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_P0P];
             f[DIR_M0P]  = ftemp[DIR_M0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_M0P];
             f[DIR_0PP]  = ftemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_0PP];
@@ -192,7 +192,7 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_PMP] = ftemp[DIR_PMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_PMP];
             f[DIR_MMP] = ftemp[DIR_MMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_MMP];
 
-            distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
+            distributions->setDistributionInvForDirection(f[DIR_00P], x1 + DX1[DIR_00M], x2 + DX2[DIR_00M], x3 + DX3[DIR_00M], DIR_00M);
             distributions->setDistributionInvForDirection(f[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
             distributions->setDistributionInvForDirection(f[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
             distributions->setDistributionInvForDirection(f[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
@@ -202,8 +202,8 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
             distributions->setDistributionInvForDirection(f[DIR_MMP], x1 + DX1[DIR_PPM], x2 + DX2[DIR_PPM], x3 + DX3[DIR_PPM], DIR_PPM);
             break;
-        case B:
-            f[B]   = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+        case DIR_00M:
+            f[DIR_00M]   = ftemp[DIR_00M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_00M];
             f[DIR_P0M]  = ftemp[DIR_P0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_P0M];
             f[DIR_M0M]  = ftemp[DIR_M0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_M0M];
             f[DIR_0PM]  = ftemp[DIR_0PM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_0PM];
@@ -213,7 +213,7 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_PMM];
             f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_MMM];
 
-            distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
+            distributions->setDistributionInvForDirection(f[DIR_00M], x1 + DX1[DIR_00P], x2 + DX2[DIR_00P], x3 + DX3[DIR_00P], DIR_00P);
             distributions->setDistributionInvForDirection(f[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
             distributions->setDistributionInvForDirection(f[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
             distributions->setDistributionInvForDirection(f[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNonReflectingOutflowBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNonReflectingOutflowBCAlgorithm.cpp
index 27d9bd363..ed90cc759 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNonReflectingOutflowBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNonReflectingOutflowBCAlgorithm.cpp
@@ -79,12 +79,12 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
    int direction = -1;
 
    //flag points in direction of fluid
-   if (bcPtr->hasDensityBoundaryFlag(E)) { nx1 += 1; direction = E; }
-   else if (bcPtr->hasDensityBoundaryFlag(W)) { nx1 -= 1; direction = W; }
-   else if (bcPtr->hasDensityBoundaryFlag(N)) { nx2 += 1; direction = N; }
-   else if (bcPtr->hasDensityBoundaryFlag(S)) { nx2 -= 1; direction = S; }
-   else if (bcPtr->hasDensityBoundaryFlag(T)) { nx3 += 1; direction = T; }
-   else if (bcPtr->hasDensityBoundaryFlag(B)) { nx3 -= 1; direction = B; }
+   if (bcPtr->hasDensityBoundaryFlag(DIR_P00)) { nx1 += 1; direction = DIR_P00; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_M00)) { nx1 -= 1; direction = DIR_M00; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_0P0)) { nx2 += 1; direction = DIR_0P0; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_0M0)) { nx2 -= 1; direction = DIR_0M0; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_00P)) { nx3 += 1; direction = DIR_00P; }
+   else if (bcPtr->hasDensityBoundaryFlag(DIR_00M)) { nx3 -= 1; direction = DIR_00M; }
    else UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on density boundary..."));
 
    distributions->getDistribution(f, x1, x2, x3);
@@ -95,9 +95,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
 
    switch (direction)
    {
-   case E:
+   case DIR_P00:
       f[DIR_P00] = ftemp[DIR_P00] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P00];
-      f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+      f[DIR_PP0] = ftemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PP0];
       f[DIR_PM0] = ftemp[DIR_PM0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PM0];
       f[DIR_P0P] = ftemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P0P];
       f[DIR_P0M] = ftemp[DIR_P0M] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_P0M];
@@ -106,8 +106,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PPM] = ftemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PPM];
       f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[DIR_PMM];
 
-      distributions->setDistributionInvForDirection(f[DIR_P00], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
-      distributions->setDistributionInvForDirection(f[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+      distributions->setDistributionInvForDirection(f[DIR_P00], x1 + DX1[DIR_M00], x2 + DX2[DIR_M00], x3 + DX3[DIR_M00], DIR_M00);
+      distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
       distributions->setDistributionInvForDirection(f[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
       distributions->setDistributionInvForDirection(f[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
       distributions->setDistributionInvForDirection(f[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
@@ -116,8 +116,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
       distributions->setDistributionInvForDirection(f[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
       break;
-   case W:
-      f[W] = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+   case DIR_M00:
+      f[DIR_M00] = ftemp[DIR_M00] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_M00];
       f[DIR_MP0] = ftemp[DIR_MP0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MP0];
       f[DIR_MM0] = ftemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MM0];
       f[DIR_M0P] = ftemp[DIR_M0P] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_M0P];
@@ -127,9 +127,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_MPM] = ftemp[DIR_MPM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MPM];
       f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[DIR_MMM];
 
-      distributions->setDistributionInvForDirection(f[W], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], E);
+      distributions->setDistributionInvForDirection(f[DIR_M00], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], DIR_P00);
       distributions->setDistributionInvForDirection(f[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
-      distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+      distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
       distributions->setDistributionInvForDirection(f[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
       distributions->setDistributionInvForDirection(f[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
       distributions->setDistributionInvForDirection(f[DIR_MPP], x1 + DX1[DIR_PMM], x2 + DX2[DIR_PMM], x3 + DX3[DIR_PMM], DIR_PMM);
@@ -137,9 +137,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
       distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
       break;
-   case N:
-      f[N] = ftemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[N];
-      f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NE];
+   case DIR_0P0:
+      f[DIR_0P0] = ftemp[DIR_0P0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0P0];
+      f[DIR_PP0] = ftemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_PP0];
       f[DIR_MP0] = ftemp[DIR_MP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_MP0];
       f[DIR_0PP] = ftemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0PP];
       f[DIR_0PM] = ftemp[DIR_0PM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_0PM];
@@ -148,8 +148,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PPM] = ftemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_PPM];
       f[DIR_MPM] = ftemp[DIR_MPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[DIR_MPM];
 
-      distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
-      distributions->setDistributionInvForDirection(f[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+      distributions->setDistributionInvForDirection(f[DIR_0P0], x1 + DX1[DIR_0M0], x2 + DX2[DIR_0M0], x3 + DX3[DIR_0M0], DIR_0M0);
+      distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
       distributions->setDistributionInvForDirection(f[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
       distributions->setDistributionInvForDirection(f[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
       distributions->setDistributionInvForDirection(f[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
@@ -158,8 +158,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
       distributions->setDistributionInvForDirection(f[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
       break;
-   case S:
-      f[S] = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+   case DIR_0M0:
+      f[DIR_0M0] = ftemp[DIR_0M0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_0M0];
       f[DIR_PM0] = ftemp[DIR_PM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_PM0];
       f[DIR_MM0] = ftemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_MM0];
       f[DIR_0MP] = ftemp[DIR_0MP] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_0MP];
@@ -169,9 +169,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_PMM];
       f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[DIR_MMM];
 
-      distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
+      distributions->setDistributionInvForDirection(f[DIR_0M0], x1 + DX1[DIR_0P0], x2 + DX2[DIR_0P0], x3 + DX3[DIR_0P0], DIR_0P0);
       distributions->setDistributionInvForDirection(f[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
-      distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+      distributions->setDistributionInvForDirection(f[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
       distributions->setDistributionInvForDirection(f[DIR_0MP], x1 + DX1[DIR_0PM], x2 + DX2[DIR_0PM], x3 + DX3[DIR_0PM], DIR_0PM);
       distributions->setDistributionInvForDirection(f[DIR_0MM], x1 + DX1[DIR_0PP], x2 + DX2[DIR_0PP], x3 + DX3[DIR_0PP], DIR_0PP);
       distributions->setDistributionInvForDirection(f[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
@@ -179,8 +179,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
       distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
       break;
-   case T:
-      f[T] = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+   case DIR_00P:
+      f[DIR_00P] = ftemp[DIR_00P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_00P];
       f[DIR_P0P] = ftemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_P0P];
       f[DIR_M0P] = ftemp[DIR_M0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_M0P];
       f[DIR_0PP] = ftemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_0PP];
@@ -190,7 +190,7 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMP] = ftemp[DIR_PMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_PMP];
       f[DIR_MMP] = ftemp[DIR_MMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[DIR_MMP];
 
-      distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
+      distributions->setDistributionInvForDirection(f[DIR_00P], x1 + DX1[DIR_00M], x2 + DX2[DIR_00M], x3 + DX3[DIR_00M], DIR_00M);
       distributions->setDistributionInvForDirection(f[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
       distributions->setDistributionInvForDirection(f[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
       distributions->setDistributionInvForDirection(f[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
@@ -200,8 +200,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
       distributions->setDistributionInvForDirection(f[DIR_MMP], x1 + DX1[DIR_PPM], x2 + DX2[DIR_PPM], x3 + DX3[DIR_PPM], DIR_PPM);
       break;
-   case B:
-      f[B] = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+   case DIR_00M:
+      f[DIR_00M] = ftemp[DIR_00M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_00M];
       f[DIR_P0M] = ftemp[DIR_P0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_P0M];
       f[DIR_M0M] = ftemp[DIR_M0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_M0M];
       f[DIR_0PM] = ftemp[DIR_0PM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_0PM];
@@ -211,7 +211,7 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       f[DIR_PMM] = ftemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_PMM];
       f[DIR_MMM] = ftemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[DIR_MMM];
 
-      distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
+      distributions->setDistributionInvForDirection(f[DIR_00M], x1 + DX1[DIR_00P], x2 + DX2[DIR_00P], x3 + DX3[DIR_00P], DIR_00P);
       distributions->setDistributionInvForDirection(f[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
       distributions->setDistributionInvForDirection(f[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
       distributions->setDistributionInvForDirection(f[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
@@ -239,9 +239,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
 
    switch (direction)
    {
-   case E:
+   case DIR_P00:
       h[DIR_P00]  = htemp[DIR_P00] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_P00];
-      h[NE] = htemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[NE];
+      h[DIR_PP0] = htemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_PP0];
       h[DIR_PM0] = htemp[DIR_PM0] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_PM0];
       h[DIR_P0P] = htemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_P0P];
       h[DIR_P0M] = htemp[DIR_P0M] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_P0M];
@@ -250,8 +250,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_PPM] = htemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_PPM];
       h[DIR_PMM] = htemp[DIR_PMM] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[DIR_PMM];
 
-      distributionsH->setDistributionInvForDirection(h[DIR_P00], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
-      distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+      distributionsH->setDistributionInvForDirection(h[DIR_P00], x1 + DX1[DIR_M00], x2 + DX2[DIR_M00], x3 + DX3[DIR_M00], DIR_M00);
+      distributionsH->setDistributionInvForDirection(h[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
       distributionsH->setDistributionInvForDirection(h[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
       distributionsH->setDistributionInvForDirection(h[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
       distributionsH->setDistributionInvForDirection(h[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
@@ -260,8 +260,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
       distributionsH->setDistributionInvForDirection(h[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
       break;
-   case W:
-      h[W] = htemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[W];
+   case DIR_M00:
+      h[DIR_M00] = htemp[DIR_M00] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_M00];
       h[DIR_MP0] = htemp[DIR_MP0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_MP0];
       h[DIR_MM0] = htemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_MM0];
       h[DIR_M0P] = htemp[DIR_M0P] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_M0P];
@@ -271,9 +271,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_MPM] = htemp[DIR_MPM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_MPM];
       h[DIR_MMM] = htemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[DIR_MMM];
 
-      distributionsH->setDistributionInvForDirection(h[W], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], E);
+      distributionsH->setDistributionInvForDirection(h[DIR_M00], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], DIR_P00);
       distributionsH->setDistributionInvForDirection(h[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
-      distributionsH->setDistributionInvForDirection(h[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+      distributionsH->setDistributionInvForDirection(h[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
       distributionsH->setDistributionInvForDirection(h[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
       distributionsH->setDistributionInvForDirection(h[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
       distributionsH->setDistributionInvForDirection(h[DIR_MPP], x1 + DX1[DIR_PMM], x2 + DX2[DIR_PMM], x3 + DX3[DIR_PMM], DIR_PMM);
@@ -281,9 +281,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
       distributionsH->setDistributionInvForDirection(h[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
       break;
-   case N:
-      h[N] = htemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[N];
-      h[NE] = htemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[NE];
+   case DIR_0P0:
+      h[DIR_0P0] = htemp[DIR_0P0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_0P0];
+      h[DIR_PP0] = htemp[DIR_PP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_PP0];
       h[DIR_MP0] = htemp[DIR_MP0] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_MP0];
       h[DIR_0PP] = htemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_0PP];
       h[DIR_0PM] = htemp[DIR_0PM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_0PM];
@@ -292,8 +292,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_PPM] = htemp[DIR_PPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_PPM];
       h[DIR_MPM] = htemp[DIR_MPM] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[DIR_MPM];
 
-      distributionsH->setDistributionInvForDirection(h[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
-      distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
+      distributionsH->setDistributionInvForDirection(h[DIR_0P0], x1 + DX1[DIR_0M0], x2 + DX2[DIR_0M0], x3 + DX3[DIR_0M0], DIR_0M0);
+      distributionsH->setDistributionInvForDirection(h[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
       distributionsH->setDistributionInvForDirection(h[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
       distributionsH->setDistributionInvForDirection(h[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
       distributionsH->setDistributionInvForDirection(h[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
@@ -302,8 +302,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[DIR_PPM], x1 + DX1[DIR_MMP], x2 + DX2[DIR_MMP], x3 + DX3[DIR_MMP], DIR_MMP);
       distributionsH->setDistributionInvForDirection(h[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
       break;
-   case S:
-      h[S] = htemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[S];
+   case DIR_0M0:
+      h[DIR_0M0] = htemp[DIR_0M0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_0M0];
       h[DIR_PM0] = htemp[DIR_PM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_PM0];
       h[DIR_MM0] = htemp[DIR_MM0] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_MM0];
       h[DIR_0MP] = htemp[DIR_0MP] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_0MP];
@@ -313,9 +313,9 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_PMM] = htemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_PMM];
       h[DIR_MMM] = htemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[DIR_MMM];
 
-      distributionsH->setDistributionInvForDirection(h[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
+      distributionsH->setDistributionInvForDirection(h[DIR_0M0], x1 + DX1[DIR_0P0], x2 + DX2[DIR_0P0], x3 + DX3[DIR_0P0], DIR_0P0);
       distributionsH->setDistributionInvForDirection(h[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
-      distributionsH->setDistributionInvForDirection(h[DIR_MM0], x1 + DX1[NE], x2 + DX2[NE], x3 + DX3[NE], NE);
+      distributionsH->setDistributionInvForDirection(h[DIR_MM0], x1 + DX1[DIR_PP0], x2 + DX2[DIR_PP0], x3 + DX3[DIR_PP0], DIR_PP0);
       distributionsH->setDistributionInvForDirection(h[DIR_0MP], x1 + DX1[DIR_0PM], x2 + DX2[DIR_0PM], x3 + DX3[DIR_0PM], DIR_0PM);
       distributionsH->setDistributionInvForDirection(h[DIR_0MM], x1 + DX1[DIR_0PP], x2 + DX2[DIR_0PP], x3 + DX3[DIR_0PP], DIR_0PP);
       distributionsH->setDistributionInvForDirection(h[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
@@ -323,8 +323,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
       distributionsH->setDistributionInvForDirection(h[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
       break;
-   case T:
-      h[T] = htemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[T];
+   case DIR_00P:
+      h[DIR_00P] = htemp[DIR_00P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_00P];
       h[DIR_P0P] = htemp[DIR_P0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_P0P];
       h[DIR_M0P] = htemp[DIR_M0P] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_M0P];
       h[DIR_0PP] = htemp[DIR_0PP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_0PP];
@@ -334,7 +334,7 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_PMP] = htemp[DIR_PMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_PMP];
       h[DIR_MMP] = htemp[DIR_MMP] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[DIR_MMP];
 
-      distributionsH->setDistributionInvForDirection(h[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
+      distributionsH->setDistributionInvForDirection(h[DIR_00P], x1 + DX1[DIR_00M], x2 + DX2[DIR_00M], x3 + DX3[DIR_00M], DIR_00M);
       distributionsH->setDistributionInvForDirection(h[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
       distributionsH->setDistributionInvForDirection(h[DIR_M0P], x1 + DX1[DIR_P0M], x2 + DX2[DIR_P0M], x3 + DX3[DIR_P0M], DIR_P0M);
       distributionsH->setDistributionInvForDirection(h[DIR_0PP], x1 + DX1[DIR_0MM], x2 + DX2[DIR_0MM], x3 + DX3[DIR_0MM], DIR_0MM);
@@ -344,8 +344,8 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[DIR_PMP], x1 + DX1[DIR_MPM], x2 + DX2[DIR_MPM], x3 + DX3[DIR_MPM], DIR_MPM);
       distributionsH->setDistributionInvForDirection(h[DIR_MMP], x1 + DX1[DIR_PPM], x2 + DX2[DIR_PPM], x3 + DX3[DIR_PPM], DIR_PPM);
       break;
-   case B:
-      h[B] = htemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[B];
+   case DIR_00M:
+      h[DIR_00M] = htemp[DIR_00M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_00M];
       h[DIR_P0M] = htemp[DIR_P0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_P0M];
       h[DIR_M0M] = htemp[DIR_M0M] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_M0M];
       h[DIR_0PM] = htemp[DIR_0PM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_0PM];
@@ -355,7 +355,7 @@ void ThixotropyNonReflectingOutflowBCAlgorithm::applyBC()
       h[DIR_PMM] = htemp[DIR_PMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_PMM];
       h[DIR_MMM] = htemp[DIR_MMM] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[DIR_MMM];
 
-      distributionsH->setDistributionInvForDirection(h[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
+      distributionsH->setDistributionInvForDirection(h[DIR_00M], x1 + DX1[DIR_00P], x2 + DX2[DIR_00P], x3 + DX3[DIR_00P], DIR_00P);
       distributionsH->setDistributionInvForDirection(h[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
       distributionsH->setDistributionInvForDirection(h[DIR_M0M], x1 + DX1[DIR_P0P], x2 + DX2[DIR_P0P], x3 + DX3[DIR_P0P], DIR_P0P);
       distributionsH->setDistributionInvForDirection(h[DIR_0PM], x1 + DX1[DIR_0MP], x2 + DX2[DIR_0MP], x3 + DX3[DIR_0MP], DIR_0MP);
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/ShearStressCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/ShearStressCoProcessor.cpp
index 9d865e22b..64ecc177f 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/ShearStressCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/ShearStressCoProcessor.cpp
@@ -173,17 +173,17 @@ void ShearStressCoProcessor::calculateShearStress(double timeStep)
                     // compute velocity
                     //////////////////////////////////////////////////////////////////////////
                     vx = ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_PMM] - f[DIR_MPP]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                          (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[NE] - f[DIR_MM0]))) + (f[DIR_P00] - f[W]));
+                          (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_P00] - f[DIR_M00]));
 
                     vy = ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_MPM] - f[DIR_PMP])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                          (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[NE] - f[DIR_MM0]))) + (f[N] - f[S]));
+                          (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_0P0] - f[DIR_0M0]));
 
                     vz = ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_MMP] - f[DIR_PPM]))) +
-                          (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[T] - f[B]));
+                          (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[DIR_00P] - f[DIR_00M]));
 
                     sxy = 3.0 * collFactor / (collFactor - 1.0) *
                           (((f[DIR_PPP] + f[DIR_MMM]) - (f[DIR_PMP] + f[DIR_MPM])) + (-(f[DIR_PMM] + f[DIR_MPP]) + (f[DIR_MMP] + f[DIR_PPM])) +
-                           (((f[NE] + f[DIR_MM0]) - (f[DIR_PM0] + f[DIR_MP0]))) - vx * vy);
+                           (((f[DIR_PP0] + f[DIR_MM0]) - (f[DIR_PM0] + f[DIR_MP0]))) - vx * vy);
 
                     sxz = 3.0 * collFactor / (collFactor - 1.0) *
                           (((f[DIR_PPP] + f[DIR_MMM]) + (f[DIR_PMP] + f[DIR_MPM])) + (-(f[DIR_PMM] + f[DIR_MPP]) - (f[DIR_MMP] + f[DIR_PPM])) +
@@ -195,11 +195,11 @@ void ShearStressCoProcessor::calculateShearStress(double timeStep)
 
                     LBMReal dxxMyy = 3.0 / 2.0 * collFactor / (collFactor - 1.0) *
                                      (((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) - ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM])) +
-                                      ((f[DIR_P00] + f[W]) - (f[N] + f[S])) - vx * vx + vy * vy);
+                                      ((f[DIR_P00] + f[DIR_M00]) - (f[DIR_0P0] + f[DIR_0M0])) - vx * vx + vy * vy);
 
                     LBMReal dxxMzz = 3.0 / 2.0 * collFactor / (collFactor - 1.0) *
-                                     ((((f[NE] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) - ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
-                                      ((f[DIR_P00] + f[W]) - (f[T] + f[B])) - vx * vx + vz * vz);
+                                     ((((f[DIR_PP0] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) - ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
+                                      ((f[DIR_P00] + f[DIR_M00]) - (f[DIR_00P] + f[DIR_00M])) - vx * vx + vz * vz);
 
                     // LBMReal dyyMzz =3.0/2.0 *collFactor/(collFactor-1.0)*((((f[NE] + f[SW]) + (f[SE] +
                     // f[NW]))-((f[TE] + f[BW])+(f[BE]+ f[TW])))
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/TurbulenceIntensityCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/TurbulenceIntensityCoProcessor.cpp
index 5de46c899..6a06a20d4 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/TurbulenceIntensityCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/TurbulenceIntensityCoProcessor.cpp
@@ -215,13 +215,13 @@ void TurbulenceIntensityCoProcessor::calculateAverageValues(double timeStep)
                                 //////////////////////////////////////////////////////////////////////////
                                 // compute velocity
                                 //////////////////////////////////////////////////////////////////////////
-                                vx = f[DIR_P00] - f[W] + f[NE] - f[DIR_MM0] + f[DIR_PM0] - f[DIR_MP0] + f[DIR_P0P] - f[DIR_M0M] + f[DIR_P0M] - f[DIR_M0P] +
+                                vx = f[DIR_P00] - f[DIR_M00] + f[DIR_PP0] - f[DIR_MM0] + f[DIR_PM0] - f[DIR_MP0] + f[DIR_P0P] - f[DIR_M0M] + f[DIR_P0M] - f[DIR_M0P] +
                                      f[DIR_PPP] - f[DIR_MMP] + f[DIR_PMP] - f[DIR_MPP] + f[DIR_PPM] - f[DIR_MMM] + f[DIR_PMM] - f[DIR_MPM];
 
-                                vy = f[N] - f[S] + f[NE] - f[DIR_MM0] - f[DIR_PM0] + f[DIR_MP0] + f[DIR_0PP] - f[DIR_0MM] + f[DIR_0PM] - f[DIR_0MP] +
+                                vy = f[DIR_0P0] - f[DIR_0M0] + f[DIR_PP0] - f[DIR_MM0] - f[DIR_PM0] + f[DIR_MP0] + f[DIR_0PP] - f[DIR_0MM] + f[DIR_0PM] - f[DIR_0MP] +
                                      f[DIR_PPP] - f[DIR_MMP] - f[DIR_PMP] + f[DIR_MPP] + f[DIR_PPM] - f[DIR_MMM] - f[DIR_PMM] + f[DIR_MPM];
 
-                                vz = f[T] - f[B] + f[DIR_P0P] - f[DIR_M0M] - f[DIR_P0M] + f[DIR_M0P] + f[DIR_0PP] - f[DIR_0MM] - f[DIR_0PM] + f[DIR_0MP] +
+                                vz = f[DIR_00P] - f[DIR_00M] + f[DIR_P0P] - f[DIR_M0M] - f[DIR_P0M] + f[DIR_M0P] + f[DIR_0PP] - f[DIR_0MM] - f[DIR_0PM] + f[DIR_0MP] +
                                      f[DIR_PPP] + f[DIR_MMP] + f[DIR_PMP] + f[DIR_MPP] - f[DIR_PPM] - f[DIR_MMM] - f[DIR_PMM] - f[DIR_MPM];
                                 //////////////////////////////////////////////////////////////////////////
                                 // compute average values
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
index 73bc94c73..0298c1dbe 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -226,16 +226,16 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     distributionsH->getDistribution(f, ix1, ix2, ix3);
                     (*phaseField)(ix1, ix2, ix3) =
                         ((f[DIR_PPP] + f[DIR_MMM]) + (f[DIR_PMP] + f[DIR_MPM])) + ((f[DIR_PMM] + f[DIR_MPP]) + (f[DIR_MMP] + f[DIR_PPM])) +
-                        (((f[NE] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
+                        (((f[DIR_PP0] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
                         ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
-                            ((f[DIR_P00] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000];
+                            ((f[DIR_P00] + f[DIR_M00]) + (f[DIR_0P0] + f[DIR_0M0]) + (f[DIR_00P] + f[DIR_00M])) + f[DIR_000];
                     if (distributionsH2) {
                     distributionsH2->getDistribution(f, ix1, ix2, ix3);
                     (*phaseField2)(ix1, ix2, ix3) =
                         ((f[DIR_PPP] + f[DIR_MMM]) + (f[DIR_PMP] + f[DIR_MPM])) + ((f[DIR_PMM] + f[DIR_MPP]) + (f[DIR_MMP] + f[DIR_PPM])) +
-                        (((f[NE] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
+                        (((f[DIR_PP0] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
                         ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
-                            ((f[DIR_P00] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000];
+                            ((f[DIR_P00] + f[DIR_M00]) + (f[DIR_0P0] + f[DIR_0M0]) + (f[DIR_00P] + f[DIR_00M])) + f[DIR_000];
                 }
                     else { (*phaseField2)(ix1, ix2, ix3) = 999.0; }
                     
@@ -285,24 +285,24 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                         // vx2 = 0.0;
                         // vx3 = 0.0;
                     } else {
-                        phi[DIR_P00]   = (*phaseField)(ix1 + DX1[DIR_P00], ix2 + DX2[DIR_P00], ix3 + DX3[DIR_P00]);
-                        phi[N]   = (*phaseField)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
-                        phi[T]   = (*phaseField)(ix1 + DX1[T], ix2 + DX2[T], ix3 + DX3[T]);
-                        phi[W]   = (*phaseField)(ix1 + DX1[W], ix2 + DX2[W], ix3 + DX3[W]);
-                        phi[S]   = (*phaseField)(ix1 + DX1[S], ix2 + DX2[S], ix3 + DX3[S]);
-                        phi[B]   = (*phaseField)(ix1 + DX1[B], ix2 + DX2[B], ix3 + DX3[B]);
-                        phi[NE]  = (*phaseField)(ix1 + DX1[NE], ix2 + DX2[NE], ix3 + DX3[NE]);
-                        phi[DIR_MP0]  = (*phaseField)(ix1 + DX1[DIR_MP0], ix2 + DX2[DIR_MP0], ix3 + DX3[DIR_MP0]);
-                        phi[DIR_P0P]  = (*phaseField)(ix1 + DX1[DIR_P0P], ix2 + DX2[DIR_P0P], ix3 + DX3[DIR_P0P]);
-                        phi[DIR_M0P]  = (*phaseField)(ix1 + DX1[DIR_M0P], ix2 + DX2[DIR_M0P], ix3 + DX3[DIR_M0P]);
-                        phi[DIR_0PP]  = (*phaseField)(ix1 + DX1[DIR_0PP], ix2 + DX2[DIR_0PP], ix3 + DX3[DIR_0PP]);
-                        phi[DIR_0MP]  = (*phaseField)(ix1 + DX1[DIR_0MP], ix2 + DX2[DIR_0MP], ix3 + DX3[DIR_0MP]);
-                        phi[DIR_MM0]  = (*phaseField)(ix1 + DX1[DIR_MM0], ix2 + DX2[DIR_MM0], ix3 + DX3[DIR_MM0]);
-                        phi[DIR_PM0]  = (*phaseField)(ix1 + DX1[DIR_PM0], ix2 + DX2[DIR_PM0], ix3 + DX3[DIR_PM0]);
-                        phi[DIR_M0M]  = (*phaseField)(ix1 + DX1[DIR_M0M], ix2 + DX2[DIR_M0M], ix3 + DX3[DIR_M0M]);
-                        phi[DIR_P0M]  = (*phaseField)(ix1 + DX1[DIR_P0M], ix2 + DX2[DIR_P0M], ix3 + DX3[DIR_P0M]);
-                        phi[DIR_0MM]  = (*phaseField)(ix1 + DX1[DIR_0MM], ix2 + DX2[DIR_0MM], ix3 + DX3[DIR_0MM]);
-                        phi[DIR_0PM]  = (*phaseField)(ix1 + DX1[DIR_0PM], ix2 + DX2[DIR_0PM], ix3 + DX3[DIR_0PM]);
+                        phi[DIR_P00] = (*phaseField)(ix1 + DX1[DIR_P00], ix2 + DX2[DIR_P00], ix3 + DX3[DIR_P00]);
+                        phi[DIR_0P0] = (*phaseField)(ix1 + DX1[DIR_0P0], ix2 + DX2[DIR_0P0], ix3 + DX3[DIR_0P0]);
+                        phi[DIR_00P] = (*phaseField)(ix1 + DX1[DIR_00P], ix2 + DX2[DIR_00P], ix3 + DX3[DIR_00P]);
+                        phi[DIR_M00] = (*phaseField)(ix1 + DX1[DIR_M00], ix2 + DX2[DIR_M00], ix3 + DX3[DIR_M00]);
+                        phi[DIR_0M0] = (*phaseField)(ix1 + DX1[DIR_0M0], ix2 + DX2[DIR_0M0], ix3 + DX3[DIR_0M0]);
+                        phi[DIR_00M] = (*phaseField)(ix1 + DX1[DIR_00M], ix2 + DX2[DIR_00M], ix3 + DX3[DIR_00M]);
+                        phi[DIR_PP0] = (*phaseField)(ix1 + DX1[DIR_PP0], ix2 + DX2[DIR_PP0], ix3 + DX3[DIR_PP0]);
+                        phi[DIR_MP0] = (*phaseField)(ix1 + DX1[DIR_MP0], ix2 + DX2[DIR_MP0], ix3 + DX3[DIR_MP0]);
+                        phi[DIR_P0P] = (*phaseField)(ix1 + DX1[DIR_P0P], ix2 + DX2[DIR_P0P], ix3 + DX3[DIR_P0P]);
+                        phi[DIR_M0P] = (*phaseField)(ix1 + DX1[DIR_M0P], ix2 + DX2[DIR_M0P], ix3 + DX3[DIR_M0P]);
+                        phi[DIR_0PP] = (*phaseField)(ix1 + DX1[DIR_0PP], ix2 + DX2[DIR_0PP], ix3 + DX3[DIR_0PP]);
+                        phi[DIR_0MP] = (*phaseField)(ix1 + DX1[DIR_0MP], ix2 + DX2[DIR_0MP], ix3 + DX3[DIR_0MP]);
+                        phi[DIR_MM0] = (*phaseField)(ix1 + DX1[DIR_MM0], ix2 + DX2[DIR_MM0], ix3 + DX3[DIR_MM0]);
+                        phi[DIR_PM0] = (*phaseField)(ix1 + DX1[DIR_PM0], ix2 + DX2[DIR_PM0], ix3 + DX3[DIR_PM0]);
+                        phi[DIR_M0M] = (*phaseField)(ix1 + DX1[DIR_M0M], ix2 + DX2[DIR_M0M], ix3 + DX3[DIR_M0M]);
+                        phi[DIR_P0M] = (*phaseField)(ix1 + DX1[DIR_P0M], ix2 + DX2[DIR_P0M], ix3 + DX3[DIR_P0M]);
+                        phi[DIR_0MM] = (*phaseField)(ix1 + DX1[DIR_0MM], ix2 + DX2[DIR_0MM], ix3 + DX3[DIR_0MM]);
+                        phi[DIR_0PM] = (*phaseField)(ix1 + DX1[DIR_0PM], ix2 + DX2[DIR_0PM], ix3 + DX3[DIR_0PM]);
                         phi[DIR_MMM] = (*phaseField)(ix1 + DX1[DIR_MMM], ix2 + DX2[DIR_MMM], ix3 + DX3[DIR_MMM]);
                         phi[DIR_PMM] = (*phaseField)(ix1 + DX1[DIR_PMM], ix2 + DX2[DIR_PMM], ix3 + DX3[DIR_PMM]);
                         phi[DIR_MPM] = (*phaseField)(ix1 + DX1[DIR_MPM], ix2 + DX2[DIR_MPM], ix3 + DX3[DIR_MPM]);
@@ -363,42 +363,42 @@ void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     if (pressure) {
                         vx1 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_PMM] - f[DIR_MPP]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                            (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[NE] - f[DIR_MM0]))) + (f[DIR_P00] - f[W])) ;
+                            (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_P00] - f[DIR_M00])) ;
 
                         vx2 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_MPM] - f[DIR_PMP])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                            (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[NE] - f[DIR_MM0]))) + (f[N] - f[S])) ;
+                            (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_0P0] - f[DIR_0M0])) ;
 
                         vx3 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_MMP] - f[DIR_PPM]))) +
-                            (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[T] - f[B]));
+                            (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[DIR_00P] - f[DIR_00M]));
 
                     }
                     else {
                         vx1 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_PMM] - f[DIR_MPP]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                            (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[NE] - f[DIR_MM0]))) + (f[DIR_P00] - f[W])) /
+                            (((f[DIR_P0M] - f[DIR_M0P]) + (f[DIR_P0P] - f[DIR_M0M])) + ((f[DIR_PM0] - f[DIR_MP0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_P00] - f[DIR_M00])) /
                                 (rho * c1o3) +
                             mu * dX1_phi / (2 * rho);
 
                         vx2 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_MPM] - f[DIR_PMP])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_PPM] - f[DIR_MMP]))) +
-                            (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[NE] - f[DIR_MM0]))) + (f[N] - f[S])) /
+                            (((f[DIR_0PM] - f[DIR_0MP]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_MP0] - f[DIR_PM0]) + (f[DIR_PP0] - f[DIR_MM0]))) + (f[DIR_0P0] - f[DIR_0M0])) /
                                 (rho * c1o3) +
                             mu * dX2_phi / (2 * rho);
 
                         vx3 =
                             ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_MMP] - f[DIR_PPM]))) +
-                            (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[T] - f[B])) /
+                            (((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[DIR_00P] - f[DIR_00M])) /
                                 (rho * c1o3) +
                             mu * dX3_phi / (2 * rho);
 
                     }
 
                     p1 = (((f[DIR_PPP] + f[DIR_MMM]) + (f[DIR_PMP] + f[DIR_MPM])) + ((f[DIR_PMM] + f[DIR_MPP]) + (f[DIR_MMP] + f[DIR_PPM])) +
-                          (((f[NE] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
+                          (((f[DIR_PP0] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) + ((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) +
                            ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
-                          ((f[DIR_P00] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[DIR_000]) +
+                          ((f[DIR_P00] + f[DIR_M00]) + (f[DIR_0P0] + f[DIR_0M0]) + (f[DIR_00P] + f[DIR_00M])) + f[DIR_000]) +
                          (vx1 * rhoToPhi * dX1_phi * c1o3 + vx2 * rhoToPhi * dX2_phi * c1o3 +
                           vx3 * rhoToPhi * dX3_phi * c1o3) /
                              2.0;
diff --git a/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
index 2d8f71843..39b83f72a 100644
--- a/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/IncompressibleOffsetInterpolationProcessor.cpp
@@ -82,11 +82,11 @@ void IncompressibleOffsetInterpolationProcessor::calcMoments(const LBMReal* cons
    //press = D3Q27System::calcPress(f,rho,vx1,vx2,vx3);
    press = rho; //interpolate rho!
 
-   kxy   = -3.*omega*((((f[DIR_MMP]+f[DIR_PPM])-(f[DIR_MPP]+f[DIR_PMM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_MPM]+f[DIR_PMP])))+((f[DIR_MM0]+f[NE])-(f[DIR_MP0]+f[DIR_PM0]))-(vx1*vx2));// might not be optimal MG 25.2.13
+   kxy   = -3.*omega*((((f[DIR_MMP]+f[DIR_PPM])-(f[DIR_MPP]+f[DIR_PMM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_MPM]+f[DIR_PMP])))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_MP0]+f[DIR_PM0]))-(vx1*vx2));// might not be optimal MG 25.2.13
    kyz   = -3.*omega*((((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMP]+f[DIR_MPM]))+((f[DIR_PMM]+f[DIR_MPP])-(f[DIR_MMP]+f[DIR_PPM])))+((f[DIR_0MM]+f[DIR_0PP])-(f[DIR_0MP]+f[DIR_0PM]))-(vx2*vx3));
    kxz   = -3.*omega*((((f[DIR_MPM]+f[DIR_PMP])-(f[DIR_MMP]+f[DIR_PPM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMM]+f[DIR_MPP])))+((f[DIR_M0M]+f[DIR_P0P])-(f[DIR_M0P]+f[DIR_P0M]))-(vx1*vx3));
-   kxxMyy = -3./2.*omega*((((f[D3Q27System::DIR_M0M]+f[DIR_P0P])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_M0P]+f[DIR_P0M])-(f[DIR_0MP]+f[DIR_0PM])))+((f[W]+f[DIR_P00])-(f[S]+f[N]))-(vx1*vx1-vx2*vx2));
-   kxxMzz = -3./2.*omega*((((f[DIR_MP0]+f[DIR_PM0])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_MM0]+f[NE])-(f[DIR_0MP]+f[DIR_0PM])))+((f[W]+f[DIR_P00])-(f[B]+f[T]))-(vx1*vx1-vx3*vx3));
+   kxxMyy = -3./2.*omega*((((f[D3Q27System::DIR_M0M]+f[DIR_P0P])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_M0P]+f[DIR_P0M])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_0M0]+f[DIR_0P0]))-(vx1*vx1-vx2*vx2));
+   kxxMzz = -3./2.*omega*((((f[DIR_MP0]+f[DIR_PM0])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_00M]+f[DIR_00P]))-(vx1*vx1-vx3*vx3));
    //kxxMzz = -3./2.*omega*(((((f[NW]+f[SE])-(f[BS]+f[TN]))+((f[SW]+f[NE])-(f[17]+f[BN])))+((f[W]+f[DIR_P00])-(f[B]+f[T])))-(vx1*vx1-vx3*vx3));
 
    //UBLOG(logINFO, "t1 = "<<(((f[NW]+f[SE])-(f[BS]+f[TN]))+((f[SW]+f[NE])-(f[17]+f[BN])))+((f[W]+f[DIR_P00])-(f[B]+f[T])));
@@ -544,12 +544,12 @@ void IncompressibleOffsetInterpolationProcessor::calcInterpolatedNode(LBMReal* f
    D3Q27System::calcIncompFeq(feq,rho,vx1,vx2,vx3);
 
    f[DIR_P00]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[DIR_P00];
-   f[W]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[W];
-   f[N]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[N];
-   f[S]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[S];
-   f[T]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[T];
-   f[B]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[B];
-   f[NE]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[NE];
+   f[DIR_M00]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[DIR_M00];
+   f[DIR_0P0]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[DIR_0P0];
+   f[DIR_0M0]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[DIR_0M0];
+   f[DIR_00P]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[DIR_00P];
+   f[DIR_00M]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[DIR_00M];
+   f[DIR_PP0]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[DIR_PP0];
    f[DIR_MM0]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[DIR_MM0];
    f[DIR_PM0]   = f_SE   + xs*x_SE   + ys*y_SE   + zs*z_SE   + xs*ys*xy_SE   + xs*zs*xz_SE   + ys*zs*yz_SE   + feq[DIR_PM0];
    f[DIR_MP0]   = f_SE   + xs*x_SE   + ys*y_SE   + zs*z_SE   + xs*ys*xy_SE   + xs*zs*xz_SE   + ys*zs*yz_SE   + feq[DIR_MP0];
@@ -738,12 +738,12 @@ void IncompressibleOffsetInterpolationProcessor::calcInterpolatedNodeFC(LBMReal*
    f_TNW = eps_new*((ay + az + bx - bz + cx - cy+kxyAverage+kxzAverage-kyzAverage)/(72.*o));
 
    f[DIR_P00]    = f_E    + feq[DIR_P00];
-   f[W]    = f_E    + feq[W];
-   f[N]    = f_N    + feq[N];
-   f[S]    = f_N    + feq[S];
-   f[T]    = f_T    + feq[T];
-   f[B]    = f_T    + feq[B];
-   f[NE]   = f_NE   + feq[NE];
+   f[DIR_M00]    = f_E    + feq[DIR_M00];
+   f[DIR_0P0]    = f_N    + feq[DIR_0P0];
+   f[DIR_0M0]    = f_N    + feq[DIR_0M0];
+   f[DIR_00P]    = f_T    + feq[DIR_00P];
+   f[DIR_00M]    = f_T    + feq[DIR_00M];
+   f[DIR_PP0]   = f_NE   + feq[DIR_PP0];
    f[DIR_MM0]   = f_NE   + feq[DIR_MM0];
    f[DIR_PM0]   = f_SE   + feq[DIR_PM0];
    f[DIR_MP0]   = f_SE   + feq[DIR_MP0];
diff --git a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
index 9a5396636..c37571337 100644
--- a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
@@ -897,9 +897,9 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                f[DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
 
                f[DIR_P00] = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
-               f[N] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
-               f[T] = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
-               f[NE] = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
+               f[DIR_0P0] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
+               f[DIR_00P] = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
+               f[DIR_PP0] = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
                f[DIR_MP0] = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
                f[DIR_P0P] = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
                f[DIR_M0P] = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -910,9 +910,9 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                f[DIR_PMP] = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
                f[DIR_MMP] = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-               f[W] = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
-               f[S] = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
-               f[B] = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
+               f[DIR_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
+               f[DIR_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
+               f[DIR_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
                f[DIR_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
                f[DIR_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
                f[DIR_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
@@ -926,9 +926,9 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                //////////////////////////////////////////////////////////////////////////
 
                drho = ((f[DIR_PPP]+f[DIR_MMM])+(f[DIR_PMP]+f[DIR_MPM]))+((f[DIR_PMM]+f[DIR_MPP])+(f[DIR_MMP]+f[DIR_PPM]))
-                  +(((f[NE]+f[DIR_MM0])+(f[DIR_PM0]+f[DIR_MP0]))+((f[DIR_P0P]+f[DIR_M0M])+(f[DIR_P0M]+f[DIR_M0P]))
-                     +((f[DIR_0PM]+f[DIR_0MP])+(f[DIR_0PP]+f[DIR_0MM])))+((f[DIR_P00]+f[W])+(f[N]+f[S])
-                        +(f[T]+f[B]))+f[DIR_000];
+                  +(((f[DIR_PP0]+f[DIR_MM0])+(f[DIR_PM0]+f[DIR_MP0]))+((f[DIR_P0P]+f[DIR_M0M])+(f[DIR_P0M]+f[DIR_M0P]))
+                     +((f[DIR_0PM]+f[DIR_0MP])+(f[DIR_0PP]+f[DIR_0MM])))+((f[DIR_P00]+f[DIR_M00])+(f[DIR_0P0]+f[DIR_0M0])
+                        +(f[DIR_00P]+f[DIR_00M]))+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])))+
@@ -958,12 +958,12 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
 
                feq[DIR_000] = c8o27*(drho-cu_sq);
                feq[DIR_P00] = 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);
-               feq[S] = c2o27*(drho+3.0*(-vx2)+c9o2*(-vx2)*(-vx2)-cu_sq);
-               feq[T] = c2o27*(drho+3.0*(vx3)+c9o2*(vx3)*(vx3)-cu_sq);
-               feq[B] = c2o27*(drho+3.0*(-vx3)+c9o2*(-vx3)*(-vx3)-cu_sq);
-               feq[NE] = c1o54*(drho+3.0*(vx1+vx2)+c9o2*(vx1+vx2)*(vx1+vx2)-cu_sq);
+               feq[DIR_M00] = c2o27*(drho+3.0*(-vx1)+c9o2*(-vx1)*(-vx1)-cu_sq);
+               feq[DIR_0P0] = c2o27*(drho+3.0*(vx2)+c9o2*(vx2)*(vx2)-cu_sq);
+               feq[DIR_0M0] = c2o27*(drho+3.0*(-vx2)+c9o2*(-vx2)*(-vx2)-cu_sq);
+               feq[DIR_00P] = c2o27*(drho+3.0*(vx3)+c9o2*(vx3)*(vx3)-cu_sq);
+               feq[DIR_00M] = c2o27*(drho+3.0*(-vx3)+c9o2*(-vx3)*(-vx3)-cu_sq);
+               feq[DIR_PP0] = c1o54*(drho+3.0*(vx1+vx2)+c9o2*(vx1+vx2)*(vx1+vx2)-cu_sq);
                feq[DIR_MM0] = c1o54*(drho+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
                feq[DIR_PM0] = c1o54*(drho+3.0*(vx1-vx2)+c9o2*(vx1-vx2)*(vx1-vx2)-cu_sq);
                feq[DIR_MP0] = c1o54*(drho+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
@@ -987,12 +987,12 @@ void InitDensityLBMKernel::calculate(int  /*step*/)
                //Relaxation
                f[DIR_000] += (feq[DIR_000]-f[DIR_000])*collFactor;
                f[DIR_P00] += (feq[DIR_P00]-f[DIR_P00])*collFactor;
-               f[W] += (feq[W]-f[W])*collFactor;
-               f[N] += (feq[N]-f[N])*collFactor;
-               f[S] += (feq[S]-f[S])*collFactor;
-               f[T] += (feq[T]-f[T])*collFactor;
-               f[B] += (feq[B]-f[B])*collFactor;
-               f[NE] += (feq[NE]-f[NE])*collFactor;
+               f[DIR_M00] += (feq[DIR_M00]-f[DIR_M00])*collFactor;
+               f[DIR_0P0] += (feq[DIR_0P0]-f[DIR_0P0])*collFactor;
+               f[DIR_0M0] += (feq[DIR_0M0]-f[DIR_0M0])*collFactor;
+               f[DIR_00P] += (feq[DIR_00P]-f[DIR_00P])*collFactor;
+               f[DIR_00M] += (feq[DIR_00M]-f[DIR_00M])*collFactor;
+               f[DIR_PP0] += (feq[DIR_PP0]-f[DIR_PP0])*collFactor;
                f[DIR_MM0] += (feq[DIR_MM0]-f[DIR_MM0])*collFactor;
                f[DIR_PM0] += (feq[DIR_PM0]-f[DIR_PM0])*collFactor;
                f[DIR_MP0] += (feq[DIR_MP0]-f[DIR_MP0])*collFactor;
diff --git a/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp b/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
index 5a4034ec0..1fcdf118f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/LBMKernelETD3Q27BGK.cpp
@@ -91,9 +91,9 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                f[DIR_000] = (*this->zeroDistributions)(x1,x2,x3);
 
                f[DIR_P00] = (*this->localDistributions)(D3Q27System::ET_E, x1,x2,x3);
-               f[N] = (*this->localDistributions)(D3Q27System::ET_N,x1,x2,x3);  
-               f[T] = (*this->localDistributions)(D3Q27System::ET_T,x1,x2,x3);
-               f[NE] = (*this->localDistributions)(D3Q27System::ET_NE,x1,x2,x3);
+               f[DIR_0P0] = (*this->localDistributions)(D3Q27System::ET_N,x1,x2,x3);
+               f[DIR_00P] = (*this->localDistributions)(D3Q27System::ET_T,x1,x2,x3);
+               f[DIR_PP0] = (*this->localDistributions)(D3Q27System::ET_NE,x1,x2,x3);
                f[DIR_MP0] = (*this->localDistributions)(D3Q27System::ET_NW,x1p,x2,x3);
                f[DIR_P0P] = (*this->localDistributions)(D3Q27System::ET_TE,x1,x2,x3);
                f[DIR_M0P] = (*this->localDistributions)(D3Q27System::ET_TW, x1p,x2,x3);
@@ -104,9 +104,9 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                f[DIR_PMP] = (*this->localDistributions)(D3Q27System::ET_TSE,x1,x2p,x3);
                f[DIR_MMP] = (*this->localDistributions)(D3Q27System::ET_TSW,x1p,x2p,x3);
 
-               f[W ] = (*this->nonLocalDistributions)(D3Q27System::ET_W,x1p,x2,x3  );
-               f[S ] = (*this->nonLocalDistributions)(D3Q27System::ET_S,x1,x2p,x3  );
-               f[B ] = (*this->nonLocalDistributions)(D3Q27System::ET_B,x1,x2,x3p  );
+               f[DIR_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_W,x1p,x2,x3  );
+               f[DIR_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_S,x1,x2p,x3  );
+               f[DIR_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_B,x1,x2,x3p  );
                f[DIR_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SW,x1p,x2p,x3 );
                f[DIR_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SE,x1,x2p,x3 );
                f[DIR_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BW,x1p,x2,x3p );
@@ -119,20 +119,20 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                f[DIR_PPM] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE,x1,x2,x3p);
                //////////////////////////////////////////////////////////////////////////
 
-               drho = f[DIR_000] + f[DIR_P00] + f[W] + f[N] + f[S] + f[T] + f[B] 
-               + f[NE] + f[DIR_MM0] + f[DIR_PM0] + f[DIR_MP0] + f[DIR_P0P] + f[DIR_M0M] + f[DIR_P0M]
+               drho = f[DIR_000] + f[DIR_P00] + f[DIR_M00] + f[DIR_0P0] + f[DIR_0M0] + f[DIR_00P] + f[DIR_00M]
+               + f[DIR_PP0] + f[DIR_MM0] + f[DIR_PM0] + f[DIR_MP0] + f[DIR_P0P] + f[DIR_M0M] + f[DIR_P0M]
                + f[DIR_M0P] + f[DIR_0PP] + f[DIR_0MM] + f[DIR_0PM] + f[DIR_0MP] + f[DIR_PPP] + f[DIR_MMP]
                + f[DIR_PMP] + f[DIR_MPP] + f[DIR_PPM] + f[DIR_MMM] + f[DIR_PMM] + f[DIR_MPM];
 
-               vx1 = f[DIR_P00] - f[W] + f[NE] - f[DIR_MM0] + f[DIR_PM0] - f[DIR_MP0] + f[DIR_P0P] - f[DIR_M0M]
+               vx1 = f[DIR_P00] - f[DIR_M00] + f[DIR_PP0] - f[DIR_MM0] + f[DIR_PM0] - f[DIR_MP0] + f[DIR_P0P] - f[DIR_M0M]
                + f[DIR_P0M] - f[DIR_M0P] + f[DIR_PPP] - f[DIR_MMP] + f[DIR_PMP] - f[DIR_MPP] + f[DIR_PPM] - f[DIR_MMM]
                + f[DIR_PMM] - f[DIR_MPM]; 
 
-               vx2 = f[N] - f[S] + f[NE] - f[DIR_MM0] - f[DIR_PM0] + f[DIR_MP0] + f[DIR_0PP] - f[DIR_0MM] + f[DIR_0PM]
+               vx2 = f[DIR_0P0] - f[DIR_0M0] + f[DIR_PP0] - f[DIR_MM0] - f[DIR_PM0] + f[DIR_MP0] + f[DIR_0PP] - f[DIR_0MM] + f[DIR_0PM]
                - f[DIR_0MP] + f[DIR_PPP] - f[DIR_MMP] - f[DIR_PMP] + f[DIR_MPP] + f[DIR_PPM] - f[DIR_MMM] - f[DIR_PMM] 
                + f[DIR_MPM]; 
 
-               vx3 = f[T] - f[B] + f[DIR_P0P] - f[DIR_M0M] - f[DIR_P0M] + f[DIR_M0P] + f[DIR_0PP] - f[DIR_0MM] - f[DIR_0PM] 
+               vx3 = f[DIR_00P] - f[DIR_00M] + f[DIR_P0P] - f[DIR_M0M] - f[DIR_P0M] + f[DIR_M0P] + f[DIR_0PP] - f[DIR_0MM] - f[DIR_0PM]
                + f[DIR_0MP] + f[DIR_PPP] + f[DIR_MMP] + f[DIR_PMP] + f[DIR_MPP] - f[DIR_PPM] - f[DIR_MMM] - f[DIR_PMM] 
                - f[DIR_MPM];
 
@@ -140,12 +140,12 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
 
                feq[DIR_000] =  c8o27*(drho-cu_sq);
                feq[DIR_P00] =  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);
-               feq[S] =  c2o27*(drho+3.0*(   -vx2)+c9o2*(   -vx2)*(   -vx2)-cu_sq);
-               feq[T] =  c2o27*(drho+3.0*( vx3   )+c9o2*(    vx3)*(    vx3)-cu_sq);
-               feq[B] =  c2o27*(drho+3.0*(   -vx3)+c9o2*(   -vx3)*(   -vx3)-cu_sq);
-               feq[NE] = c1o54*(drho+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
+               feq[DIR_M00] =  c2o27*(drho+3.0*(-vx1   )+c9o2*(-vx1   )*(-vx1   )-cu_sq);
+               feq[DIR_0P0] =  c2o27*(drho+3.0*(    vx2)+c9o2*(    vx2)*(    vx2)-cu_sq);
+               feq[DIR_0M0] =  c2o27*(drho+3.0*(   -vx2)+c9o2*(   -vx2)*(   -vx2)-cu_sq);
+               feq[DIR_00P] =  c2o27*(drho+3.0*( vx3   )+c9o2*(    vx3)*(    vx3)-cu_sq);
+               feq[DIR_00M] =  c2o27*(drho+3.0*(   -vx3)+c9o2*(   -vx3)*(   -vx3)-cu_sq);
+               feq[DIR_PP0] = c1o54*(drho+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
                feq[DIR_MM0] = c1o54*(drho+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
                feq[DIR_PM0] = c1o54*(drho+3.0*( vx1-vx2)+c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq);
                feq[DIR_MP0] = c1o54*(drho+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
@@ -169,12 +169,12 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                //Relaxation
                f[DIR_000] += (feq[DIR_000]-f[DIR_000])*collFactor;
                f[DIR_P00] += (feq[DIR_P00]-f[DIR_P00])*collFactor;
-               f[W] += (feq[W]-f[W])*collFactor;
-               f[N] += (feq[N]-f[N])*collFactor;
-               f[S] += (feq[S]-f[S])*collFactor;
-               f[T] += (feq[T]-f[T])*collFactor;
-               f[B] += (feq[B]-f[B])*collFactor;
-               f[NE] += (feq[NE]-f[NE])*collFactor;
+               f[DIR_M00] += (feq[DIR_M00]-f[DIR_M00])*collFactor;
+               f[DIR_0P0] += (feq[DIR_0P0]-f[DIR_0P0])*collFactor;
+               f[DIR_0M0] += (feq[DIR_0M0]-f[DIR_0M0])*collFactor;
+               f[DIR_00P] += (feq[DIR_00P]-f[DIR_00P])*collFactor;
+               f[DIR_00M] += (feq[DIR_00M]-f[DIR_00M])*collFactor;
+               f[DIR_PP0] += (feq[DIR_PP0]-f[DIR_PP0])*collFactor;
                f[DIR_MM0] += (feq[DIR_MM0]-f[DIR_MM0])*collFactor;
                f[DIR_PM0] += (feq[DIR_PM0]-f[DIR_PM0])*collFactor;
                f[DIR_MP0] += (feq[DIR_MP0]-f[DIR_MP0])*collFactor;
@@ -209,13 +209,13 @@ void LBMKernelETD3Q27BGK::calculate(int  /*step*/)
                   forcingX3 = muForcingX3.Eval();
 
                   f[DIR_000] +=                   0.0                        ;
-                  f[E  ] +=  3.0*c2o27  *  (forcingX1)                    ;
-                  f[W  ] +=  3.0*c2o27  *  (-forcingX1)                   ;
-                  f[N  ] +=  3.0*c2o27  *             (forcingX2)         ;
-                  f[S  ] +=  3.0*c2o27  *             (-forcingX2)        ;
-                  f[T  ] +=  3.0*c2o27  *                     (forcingX3) ;
-                  f[B  ] +=  3.0*c2o27  *                     (-forcingX3);
-                  f[NE ] +=  3.0*c1o54 * ( forcingX1+forcingX2          ) ;
+                  f[DIR_P00] +=  3.0*c2o27  *  (forcingX1)                    ;
+                  f[DIR_M00] +=  3.0*c2o27  *  (-forcingX1)                   ;
+                  f[DIR_0P0] +=  3.0*c2o27  *             (forcingX2)         ;
+                  f[DIR_0M0] +=  3.0*c2o27  *             (-forcingX2)        ;
+                  f[DIR_00P] +=  3.0*c2o27  *                     (forcingX3) ;
+                  f[DIR_00M] +=  3.0*c2o27  *                     (-forcingX3);
+                  f[DIR_PP0] +=  3.0*c1o54 * ( forcingX1+forcingX2          ) ;
                   f[DIR_MM0 ] +=  3.0*c1o54 * (-forcingX1-forcingX2          ) ;
                   f[DIR_PM0 ] +=  3.0*c1o54 * ( forcingX1-forcingX2          ) ;
                   f[DIR_MP0 ] +=  3.0*c1o54 * (-forcingX1+forcingX2          ) ;
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
index b1270904a..10f1e395d 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseCumulantLBMKernel.cpp
@@ -317,9 +317,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         //--------------------------------------------------------
 
                         mfcbb = 3.0 * (mfcbb + 0.5 * forcingTerm[DIR_P00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[E  ];
-                        mfbcb = 3.0 * (mfbcb + 0.5 * forcingTerm[N]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
-                        mfbbc = 3.0 * (mfbbc + 0.5 * forcingTerm[T]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
-                        mfccb = 3.0 * (mfccb + 0.5 * forcingTerm[NE]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
+                        mfbcb = 3.0 * (mfbcb + 0.5 * forcingTerm[DIR_0P0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
+                        mfbbc = 3.0 * (mfbbc + 0.5 * forcingTerm[DIR_00P]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
+                        mfccb = 3.0 * (mfccb + 0.5 * forcingTerm[DIR_PP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
                         mfacb = 3.0 * (mfacb + 0.5 * forcingTerm[DIR_MP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NW ];
                         mfcbc = 3.0 * (mfcbc + 0.5 * forcingTerm[DIR_P0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TE ];
                         mfabc = 3.0 * (mfabc + 0.5 * forcingTerm[DIR_M0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TW ];
@@ -329,9 +329,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         mfacc = 3.0 * (mfacc + 0.5 * forcingTerm[DIR_MPP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TNW];
                         mfcac = 3.0 * (mfcac + 0.5 * forcingTerm[DIR_PMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSE];
                         mfaac = 3.0 * (mfaac + 0.5 * forcingTerm[DIR_MMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSW];
-                        mfabb = 3.0 * (mfabb + 0.5 * forcingTerm[W]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
-                        mfbab = 3.0 * (mfbab + 0.5 * forcingTerm[S]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
-                        mfbba = 3.0 * (mfbba + 0.5 * forcingTerm[B]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
+                        mfabb = 3.0 * (mfabb + 0.5 * forcingTerm[DIR_M00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
+                        mfbab = 3.0 * (mfbab + 0.5 * forcingTerm[DIR_0M0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
+                        mfbba = 3.0 * (mfbba + 0.5 * forcingTerm[DIR_00M]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
                         mfaab = 3.0 * (mfaab + 0.5 * forcingTerm[DIR_MM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SW ];
                         mfcab = 3.0 * (mfcab + 0.5 * forcingTerm[DIR_PM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SE ];
                         mfaba = 3.0 * (mfaba + 0.5 * forcingTerm[DIR_M0M]) / rho;   //-(3.0*p1 - rho)*WEIGTH[BW ];
@@ -1025,9 +1025,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
 #endif
 
                         mfcbb = rho * c1o3 * (mfcbb) + 0.5 * forcingTerm[DIR_P00];
-                        mfbcb = rho * c1o3 * (mfbcb) + 0.5 * forcingTerm[N];
-                        mfbbc = rho * c1o3 * (mfbbc) + 0.5 * forcingTerm[T];
-                        mfccb = rho * c1o3 * (mfccb) + 0.5 * forcingTerm[NE];
+                        mfbcb = rho * c1o3 * (mfbcb) + 0.5 * forcingTerm[DIR_0P0];
+                        mfbbc = rho * c1o3 * (mfbbc) + 0.5 * forcingTerm[DIR_00P];
+                        mfccb = rho * c1o3 * (mfccb) + 0.5 * forcingTerm[DIR_PP0];
                         mfacb = rho * c1o3 * (mfacb) + 0.5 * forcingTerm[DIR_MP0];
                         mfcbc = rho * c1o3 * (mfcbc) + 0.5 * forcingTerm[DIR_P0P];
                         mfabc = rho * c1o3 * (mfabc) + 0.5 * forcingTerm[DIR_M0P];
@@ -1037,9 +1037,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         mfacc = rho * c1o3 * (mfacc) + 0.5 * forcingTerm[DIR_MPP];
                         mfcac = rho * c1o3 * (mfcac) + 0.5 * forcingTerm[DIR_PMP];
                         mfaac = rho * c1o3 * (mfaac) + 0.5 * forcingTerm[DIR_MMP];
-                        mfabb = rho * c1o3 * (mfabb) + 0.5 * forcingTerm[W];
-                        mfbab = rho * c1o3 * (mfbab) + 0.5 * forcingTerm[S];
-                        mfbba = rho * c1o3 * (mfbba) + 0.5 * forcingTerm[B];
+                        mfabb = rho * c1o3 * (mfabb) + 0.5 * forcingTerm[DIR_M00];
+                        mfbab = rho * c1o3 * (mfbab) + 0.5 * forcingTerm[DIR_0M0];
+                        mfbba = rho * c1o3 * (mfbba) + 0.5 * forcingTerm[DIR_00M];
                         mfaab = rho * c1o3 * (mfaab) + 0.5 * forcingTerm[DIR_MM0];
                         mfcab = rho * c1o3 * (mfcab) + 0.5 * forcingTerm[DIR_PM0];
                         mfaba = rho * c1o3 * (mfaba) + 0.5 * forcingTerm[DIR_M0M];
@@ -1092,9 +1092,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         /////////////////////   PHASE-FIELD BGK SOLVER ///////////////////////////////
 
                         h[DIR_P00]   = (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3);
-                        h[N]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
-                        h[T]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
-                        h[NE]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
+                        h[DIR_0P0]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
+                        h[DIR_00P]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
+                        h[DIR_PP0]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
                         h[DIR_MP0]  = (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3);
                         h[DIR_P0P]  = (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3);
                         h[DIR_M0P]  = (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -1105,9 +1105,9 @@ void MultiphaseCumulantLBMKernel::calculate(int step)
                         h[DIR_PMP] = (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3);
                         h[DIR_MMP] = (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                        h[W]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
-                        h[S]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
-                        h[B]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
+                        h[DIR_M00]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
+                        h[DIR_0M0]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
+                        h[DIR_00M]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
                         h[DIR_MM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3);
                         h[DIR_PM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3);
                         h[DIR_M0M]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p);
@@ -1243,9 +1243,9 @@ void MultiphaseCumulantLBMKernel::computePhasefield()
                     int x3p = x3 + 1;
 
                     h[DIR_P00]   = (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3);
-                    h[N]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
-                    h[T]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
-                    h[NE]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
+                    h[DIR_0P0]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
+                    h[DIR_00P]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
+                    h[DIR_PP0]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
                     h[DIR_MP0]  = (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3);
                     h[DIR_P0P]  = (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3);
                     h[DIR_M0P]  = (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -1256,9 +1256,9 @@ void MultiphaseCumulantLBMKernel::computePhasefield()
                     h[DIR_PMP] = (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3);
                     h[DIR_MMP] = (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                    h[W]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
-                    h[S]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
-                    h[B]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
+                    h[DIR_M00]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
+                    h[DIR_0M0]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
+                    h[DIR_00M]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
                     h[DIR_MM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3);
                     h[DIR_PM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3);
                     h[DIR_M0M]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
index 04cfbdf2d..bd4df8aea 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterCompressibleAirLBMKernel.cpp
@@ -1476,48 +1476,48 @@ LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX1_phi()
 {
 	using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX2_phi()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX3_phi()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX1_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) + (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) + (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[NE] - phi2[DIR_MM0])))) +
-		+WEIGTH[N] * (phi2[DIR_P00] - phi2[W]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[DIR_PP0] - phi2[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_P00] - phi2[DIR_M00]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX2_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PPM] - phi2[DIR_MMP]) - (phi2[DIR_PMP] - phi2[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[NE] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
-		+WEIGTH[N] * (phi2[N] - phi2[S]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[DIR_PP0] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_0P0] - phi2[DIR_0M0]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::gradX3_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) - (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
-		+WEIGTH[N] * (phi2[T] - phi2[B]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_00P] - phi2[DIR_00M]));
 }
 
 LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::nabla2_phi()
@@ -1529,12 +1529,12 @@ LBMReal MultiphasePressureFilterCompressibleAirLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 		(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-		((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+		((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
 
 	return 6.0 * sum;
@@ -1563,9 +1563,9 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::computePhasefield()
 					int x3p = x3 + 1;
 
 					h[DIR_P00]   = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
-					h[N]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
-					h[T]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
-					h[NE]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+					h[DIR_0P0]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+					h[DIR_00P]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+					h[DIR_PP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
 					h[DIR_MP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
 					h[DIR_P0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
 					h[DIR_M0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -1576,9 +1576,9 @@ void MultiphasePressureFilterCompressibleAirLBMKernel::computePhasefield()
 					h[DIR_PMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
 					h[DIR_MMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-					h[W]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
-					h[S]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
-					h[B]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+					h[DIR_M00]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+					h[DIR_0M0]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+					h[DIR_00M]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
 					h[DIR_MM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
 					h[DIR_PM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
 					h[DIR_M0M]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
index e5fcc71d8..918a162e3 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
@@ -1625,24 +1625,24 @@ LBMReal MultiphasePressureFilterLBMKernel::gradX1_phi()
 {
 	using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
 }
 
 LBMReal MultiphasePressureFilterLBMKernel::gradX2_phi()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
 }
 
 LBMReal MultiphasePressureFilterLBMKernel::gradX3_phi()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
 }
 
 LBMReal MultiphasePressureFilterLBMKernel::nabla2_phi()
@@ -1654,12 +1654,12 @@ LBMReal MultiphasePressureFilterLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 		(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-		((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+		((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
 
 	return 6.0 * sum;
@@ -1688,9 +1688,9 @@ void MultiphasePressureFilterLBMKernel::computePhasefield()
 					int x3p = x3 + 1;
 
 					h[DIR_P00]   = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
-					h[N]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
-					h[T]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
-					h[NE]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+					h[DIR_0P0]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+					h[DIR_00P]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+					h[DIR_PP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
 					h[DIR_MP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
 					h[DIR_P0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
 					h[DIR_M0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -1701,9 +1701,9 @@ void MultiphasePressureFilterLBMKernel::computePhasefield()
 					h[DIR_PMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
 					h[DIR_MMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-					h[W]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
-					h[S]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
-					h[B]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+					h[DIR_M00]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+					h[DIR_0M0]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+					h[DIR_00M]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
 					h[DIR_MM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
 					h[DIR_PM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
 					h[DIR_M0M]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index caa9272f5..c1c2010e1 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -749,9 +749,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 
 
 			   mfcbb +=3.0 * ( 0.5 * forcingTerm[DIR_P00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[E  ];
-			   mfbcb +=3.0 * ( 0.5 * forcingTerm[N]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
-			   mfbbc +=3.0 * ( 0.5 * forcingTerm[T]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
-			   mfccb +=3.0 * ( 0.5 * forcingTerm[NE]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
+			   mfbcb +=3.0 * ( 0.5 * forcingTerm[DIR_0P0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
+			   mfbbc +=3.0 * ( 0.5 * forcingTerm[DIR_00P]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
+			   mfccb +=3.0 * ( 0.5 * forcingTerm[DIR_PP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
 			   mfacb +=3.0 * ( 0.5 * forcingTerm[DIR_MP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NW ];
 			   mfcbc +=3.0 * ( 0.5 * forcingTerm[DIR_P0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TE ];
 			   mfabc +=3.0 * ( 0.5 * forcingTerm[DIR_M0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TW ];
@@ -761,9 +761,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   mfacc +=3.0 * ( 0.5 * forcingTerm[DIR_MPP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TNW];
 			   mfcac +=3.0 * ( 0.5 * forcingTerm[DIR_PMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSE];
 			   mfaac +=3.0 * ( 0.5 * forcingTerm[DIR_MMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSW];
-			   mfabb +=3.0 * ( 0.5 * forcingTerm[W]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
-			   mfbab +=3.0 * ( 0.5 * forcingTerm[S]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
-			   mfbba +=3.0 * ( 0.5 * forcingTerm[B]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
+			   mfabb +=3.0 * ( 0.5 * forcingTerm[DIR_M00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
+			   mfbab +=3.0 * ( 0.5 * forcingTerm[DIR_0M0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
+			   mfbba +=3.0 * ( 0.5 * forcingTerm[DIR_00M]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
 			   mfaab +=3.0 * ( 0.5 * forcingTerm[DIR_MM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SW ];
 			   mfcab +=3.0 * ( 0.5 * forcingTerm[DIR_PM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SE ];
 			   mfaba +=3.0 * ( 0.5 * forcingTerm[DIR_M0M]) / rho;   //-(3.0*p1 - rho)*WEIGTH[BW ];
@@ -1620,9 +1620,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   /////classical source term 8.4.2021
 
 			   mfcbb += 3.0 * (0.5 * forcingTerm[DIR_P00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[E  ];
-			   mfbcb += 3.0 * (0.5 * forcingTerm[N]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
-			   mfbbc += 3.0 * (0.5 * forcingTerm[T]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
-			   mfccb += 3.0 * (0.5 * forcingTerm[NE]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
+			   mfbcb += 3.0 * (0.5 * forcingTerm[DIR_0P0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
+			   mfbbc += 3.0 * (0.5 * forcingTerm[DIR_00P]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
+			   mfccb += 3.0 * (0.5 * forcingTerm[DIR_PP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
 			   mfacb += 3.0 * (0.5 * forcingTerm[DIR_MP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NW ];
 			   mfcbc += 3.0 * (0.5 * forcingTerm[DIR_P0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TE ];
 			   mfabc += 3.0 * (0.5 * forcingTerm[DIR_M0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TW ];
@@ -1632,9 +1632,9 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
 			   mfacc += 3.0 * (0.5 * forcingTerm[DIR_MPP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TNW];
 			   mfcac += 3.0 * (0.5 * forcingTerm[DIR_PMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSE];
 			   mfaac += 3.0 * (0.5 * forcingTerm[DIR_MMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSW];
-			   mfabb += 3.0 * (0.5 * forcingTerm[W]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
-			   mfbab += 3.0 * (0.5 * forcingTerm[S]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
-			   mfbba += 3.0 * (0.5 * forcingTerm[B]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
+			   mfabb += 3.0 * (0.5 * forcingTerm[DIR_M00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
+			   mfbab += 3.0 * (0.5 * forcingTerm[DIR_0M0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
+			   mfbba += 3.0 * (0.5 * forcingTerm[DIR_00M]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
 			   mfaab += 3.0 * (0.5 * forcingTerm[DIR_MM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SW ];
 			   mfcab += 3.0 * (0.5 * forcingTerm[DIR_PM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SE ];
 			   mfaba += 3.0 * (0.5 * forcingTerm[DIR_M0M]) / rho;   //-(3.0*p1 - rho)*WEIGTH[BW ];
@@ -2922,8 +2922,8 @@ LBMReal MultiphaseScratchCumulantLBMKernel::gradX1_phi()
 {
     using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX1[k] * phi[k];
@@ -2935,8 +2935,8 @@ LBMReal MultiphaseScratchCumulantLBMKernel::gradX2_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX2[k] * phi[k];
@@ -2948,8 +2948,8 @@ LBMReal MultiphaseScratchCumulantLBMKernel::gradX3_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX3[k] * phi[k];
@@ -2966,12 +2966,12 @@ LBMReal MultiphaseScratchCumulantLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 			(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+			((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3002,9 +3002,9 @@ void MultiphaseScratchCumulantLBMKernel::computePhasefield()
                     int x3p = x3 + 1;
 
                     h[DIR_P00]   = (*this->localDistributionsH)(D3Q27System::ET_E, x1, x2, x3);
-                    h[N]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
-                    h[T]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
-                    h[NE]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
+                    h[DIR_0P0]   = (*this->localDistributionsH)(D3Q27System::ET_N, x1, x2, x3);
+                    h[DIR_00P]   = (*this->localDistributionsH)(D3Q27System::ET_T, x1, x2, x3);
+                    h[DIR_PP0]  = (*this->localDistributionsH)(D3Q27System::ET_NE, x1, x2, x3);
                     h[DIR_MP0]  = (*this->localDistributionsH)(D3Q27System::ET_NW, x1p, x2, x3);
                     h[DIR_P0P]  = (*this->localDistributionsH)(D3Q27System::ET_TE, x1, x2, x3);
                     h[DIR_M0P]  = (*this->localDistributionsH)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -3015,9 +3015,9 @@ void MultiphaseScratchCumulantLBMKernel::computePhasefield()
                     h[DIR_PMP] = (*this->localDistributionsH)(D3Q27System::ET_TSE, x1, x2p, x3);
                     h[DIR_MMP] = (*this->localDistributionsH)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                    h[W]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
-                    h[S]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
-                    h[B]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
+                    h[DIR_M00]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_W, x1p, x2, x3);
+                    h[DIR_0M0]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_S, x1, x2p, x3);
+                    h[DIR_00M]   = (*this->nonLocalDistributionsH)(D3Q27System::ET_B, x1, x2, x3p);
                     h[DIR_MM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SW, x1p, x2p, x3);
                     h[DIR_PM0]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_SE, x1, x2p, x3);
                     h[DIR_M0M]  = (*this->nonLocalDistributionsH)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
index 08c196105..229ff921a 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsCumulantLBMKernel.cpp
@@ -601,9 +601,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 
 
 			   mfcbb += 3.0 * (0.5 * forcingTerm[DIR_P00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[E  ];
-			   mfbcb += 3.0 * (0.5 * forcingTerm[N]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
-			   mfbbc += 3.0 * (0.5 * forcingTerm[T]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
-			   mfccb += 3.0 * (0.5 * forcingTerm[NE]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
+			   mfbcb += 3.0 * (0.5 * forcingTerm[DIR_0P0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
+			   mfbbc += 3.0 * (0.5 * forcingTerm[DIR_00P]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
+			   mfccb += 3.0 * (0.5 * forcingTerm[DIR_PP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
 			   mfacb += 3.0 * (0.5 * forcingTerm[DIR_MP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NW ];
 			   mfcbc += 3.0 * (0.5 * forcingTerm[DIR_P0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TE ];
 			   mfabc += 3.0 * (0.5 * forcingTerm[DIR_M0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TW ];
@@ -613,9 +613,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   mfacc += 3.0 * (0.5 * forcingTerm[DIR_MPP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TNW];
 			   mfcac += 3.0 * (0.5 * forcingTerm[DIR_PMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSE];
 			   mfaac += 3.0 * (0.5 * forcingTerm[DIR_MMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSW];
-			   mfabb += 3.0 * (0.5 * forcingTerm[W]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
-			   mfbab += 3.0 * (0.5 * forcingTerm[S]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
-			   mfbba += 3.0 * (0.5 * forcingTerm[B]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
+			   mfabb += 3.0 * (0.5 * forcingTerm[DIR_M00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
+			   mfbab += 3.0 * (0.5 * forcingTerm[DIR_0M0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
+			   mfbba += 3.0 * (0.5 * forcingTerm[DIR_00M]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
 			   mfaab += 3.0 * (0.5 * forcingTerm[DIR_MM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SW ];
 			   mfcab += 3.0 * (0.5 * forcingTerm[DIR_PM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SE ];
 			   mfaba += 3.0 * (0.5 * forcingTerm[DIR_M0M]) / rho;   //-(3.0*p1 - rho)*WEIGTH[BW ];
@@ -1345,9 +1345,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   /////classical source term 8.4.2021
 
 			   mfcbb += 3.0 * (0.5 * forcingTerm[DIR_P00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[E  ];
-			   mfbcb += 3.0 * (0.5 * forcingTerm[N]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
-			   mfbbc += 3.0 * (0.5 * forcingTerm[T]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
-			   mfccb += 3.0 * (0.5 * forcingTerm[NE]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
+			   mfbcb += 3.0 * (0.5 * forcingTerm[DIR_0P0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[N  ];
+			   mfbbc += 3.0 * (0.5 * forcingTerm[DIR_00P]) / rho;    //-(3.0*p1 - rho)*WEIGTH[T  ];
+			   mfccb += 3.0 * (0.5 * forcingTerm[DIR_PP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NE ];
 			   mfacb += 3.0 * (0.5 * forcingTerm[DIR_MP0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[NW ];
 			   mfcbc += 3.0 * (0.5 * forcingTerm[DIR_P0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TE ];
 			   mfabc += 3.0 * (0.5 * forcingTerm[DIR_M0P]) / rho;   //-(3.0*p1 - rho)*WEIGTH[TW ];
@@ -1357,9 +1357,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::calculate(int step)
 			   mfacc += 3.0 * (0.5 * forcingTerm[DIR_MPP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TNW];
 			   mfcac += 3.0 * (0.5 * forcingTerm[DIR_PMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSE];
 			   mfaac += 3.0 * (0.5 * forcingTerm[DIR_MMP]) / rho;  //-(3.0*p1 - rho)*WEIGTH[TSW];
-			   mfabb += 3.0 * (0.5 * forcingTerm[W]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
-			   mfbab += 3.0 * (0.5 * forcingTerm[S]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
-			   mfbba += 3.0 * (0.5 * forcingTerm[B]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
+			   mfabb += 3.0 * (0.5 * forcingTerm[DIR_M00]) / rho;    //-(3.0*p1 - rho)*WEIGTH[W  ];
+			   mfbab += 3.0 * (0.5 * forcingTerm[DIR_0M0]) / rho;    //-(3.0*p1 - rho)*WEIGTH[S  ];
+			   mfbba += 3.0 * (0.5 * forcingTerm[DIR_00M]) / rho;    //-(3.0*p1 - rho)*WEIGTH[B  ];
 			   mfaab += 3.0 * (0.5 * forcingTerm[DIR_MM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SW ];
 			   mfcab += 3.0 * (0.5 * forcingTerm[DIR_PM0]) / rho;   //-(3.0*p1 - rho)*WEIGTH[SE ];
 			   mfaba += 3.0 * (0.5 * forcingTerm[DIR_M0M]) / rho;   //-(3.0*p1 - rho)*WEIGTH[BW ];
@@ -2985,8 +2985,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX1_phi()
 {
     using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX1[k] * phi[k];
@@ -2998,8 +2998,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX2_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX2[k] * phi[k];
@@ -3011,8 +3011,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX3_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX3[k] * phi[k];
@@ -3024,8 +3024,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX1_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) + (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) + (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[NE] - phi2[DIR_MM0])))) +
-		+WEIGTH[N] * (phi2[DIR_P00] - phi2[W]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[DIR_PP0] - phi2[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_P00] - phi2[DIR_M00]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX1[k] * phi2[k];
@@ -3037,8 +3037,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX2_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PPM] - phi2[DIR_MMP]) - (phi2[DIR_PMP] - phi2[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[NE] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
-		+WEIGTH[N] * (phi2[N] - phi2[S]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[DIR_PP0] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_0P0] - phi2[DIR_0M0]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX2[k] * phi2[k];
@@ -3050,8 +3050,8 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::gradX3_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) - (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
-		+WEIGTH[N] * (phi2[T] - phi2[B]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_00P] - phi2[DIR_00M]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX3[k] * phi2[k];
@@ -3072,12 +3072,12 @@ LBMReal MultiphaseTwoPhaseFieldsCumulantLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 			(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+			((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3108,9 +3108,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::computePhasefield()
                     int x3p = x3 + 1;
 
                     h[DIR_P00]   = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
-                    h[N]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
-                    h[T]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
-                    h[NE]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+                    h[DIR_0P0]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+                    h[DIR_00P]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+                    h[DIR_PP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
                     h[DIR_MP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
                     h[DIR_P0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
                     h[DIR_M0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -3121,9 +3121,9 @@ void MultiphaseTwoPhaseFieldsCumulantLBMKernel::computePhasefield()
                     h[DIR_PMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
                     h[DIR_MMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                    h[W]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
-                    h[S]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
-                    h[B]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+                    h[DIR_M00]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+                    h[DIR_0M0]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+                    h[DIR_00M]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
                     h[DIR_MM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
                     h[DIR_PM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
                     h[DIR_M0M]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
index 2f9d1d394..3baddc4fe 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsPressureFilterLBMKernel.cpp
@@ -3350,8 +3350,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX1_phi()
 {
     using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX1[k] * phi[k];
@@ -3363,8 +3363,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX2_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX2[k] * phi[k];
@@ -3376,8 +3376,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX3_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX3[k] * phi[k];
@@ -3389,8 +3389,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX1_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) + (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) + (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[NE] - phi2[DIR_MM0])))) +
-		+WEIGTH[N] * (phi2[DIR_P00] - phi2[W]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[DIR_PP0] - phi2[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_P00] - phi2[DIR_M00]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX1[k] * phi2[k];
@@ -3402,8 +3402,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX2_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PPM] - phi2[DIR_MMP]) - (phi2[DIR_PMP] - phi2[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[NE] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
-		+WEIGTH[N] * (phi2[N] - phi2[S]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[DIR_PP0] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_0P0] - phi2[DIR_0M0]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX2[k] * phi2[k];
@@ -3415,8 +3415,8 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::gradX3_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) - (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
-		+WEIGTH[N] * (phi2[T] - phi2[B]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_00P] - phi2[DIR_00M]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX3[k] * phi2[k];
@@ -3437,12 +3437,12 @@ LBMReal MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 			(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+			((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3473,9 +3473,9 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::computePhasefield()
                     int x3p = x3 + 1;
 
                     h[DIR_P00]   = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
-                    h[N]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
-                    h[T]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
-                    h[NE]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+                    h[DIR_0P0]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+                    h[DIR_00P]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+                    h[DIR_PP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
                     h[DIR_MP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
                     h[DIR_P0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
                     h[DIR_M0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -3486,9 +3486,9 @@ void MultiphaseTwoPhaseFieldsPressureFilterLBMKernel::computePhasefield()
                     h[DIR_PMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
                     h[DIR_MMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                    h[W]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
-                    h[S]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
-                    h[B]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+                    h[DIR_M00]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+                    h[DIR_0M0]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+                    h[DIR_00M]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
                     h[DIR_MM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
                     h[DIR_PM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
                     h[DIR_M0M]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
index 7fdd5f2ee..ffed1483c 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel.cpp
@@ -3282,8 +3282,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX1_phi()
 {
     using namespace D3Q27System;
 	return 3.0* ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) + (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) + (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[NE] - phi[DIR_MM0])))) +
-		+WEIGTH[N] * (phi[DIR_P00] - phi[W]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) + (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_PM0] - phi[DIR_MP0]) + (phi[DIR_PP0] - phi[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_P00] - phi[DIR_M00]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX1[k] * phi[k];
@@ -3295,8 +3295,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX2_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PPM] - phi[DIR_MMP])- (phi[DIR_PMP] - phi[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[NE] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
-		+WEIGTH[N] * (phi[N] - phi[S]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_0PP] - phi[DIR_0MM]) + (phi[DIR_0PM] - phi[DIR_0MP])) + ((phi[DIR_PP0] - phi[DIR_MM0])- (phi[DIR_PM0] - phi[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_0P0] - phi[DIR_0M0]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX2[k] * phi[k];
@@ -3308,8 +3308,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX3_phi()
 {
     using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi[DIR_PPP] - phi[DIR_MMM]) - (phi[DIR_PMM] - phi[DIR_MPP])) + ((phi[DIR_PMP] - phi[DIR_MPM]) - (phi[DIR_PPM] - phi[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
-		+WEIGTH[N] * (phi[T] - phi[B]));
+		+ WEIGTH[DIR_PP0] * (((phi[DIR_P0P] - phi[DIR_M0M]) - (phi[DIR_P0M] - phi[DIR_M0P])) + ((phi[DIR_0MP] - phi[DIR_0PM]) + (phi[DIR_0PP] - phi[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi[DIR_00P] - phi[DIR_00M]));
     //LBMReal sum = 0.0;
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * DX3[k] * phi[k];
@@ -3321,8 +3321,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX1_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) + (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) + (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[NE] - phi2[DIR_MM0])))) +
-		+WEIGTH[N] * (phi2[DIR_P00] - phi2[W]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) + (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_PM0] - phi2[DIR_MP0]) + (phi2[DIR_PP0] - phi2[DIR_MM0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_P00] - phi2[DIR_M00]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX1[k] * phi2[k];
@@ -3334,8 +3334,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX2_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PPM] - phi2[DIR_MMP]) - (phi2[DIR_PMP] - phi2[DIR_MPM])))
-		+ WEIGTH[NE] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[NE] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
-		+WEIGTH[N] * (phi2[N] - phi2[S]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_0PP] - phi2[DIR_0MM]) + (phi2[DIR_0PM] - phi2[DIR_0MP])) + ((phi2[DIR_PP0] - phi2[DIR_MM0]) - (phi2[DIR_PM0] - phi2[DIR_MP0])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_0P0] - phi2[DIR_0M0]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX2[k] * phi2[k];
@@ -3347,8 +3347,8 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::gradX3_phi2()
 {
 	using namespace D3Q27System;
 	return 3.0 * ((WEIGTH[DIR_PPP] * (((phi2[DIR_PPP] - phi2[DIR_MMM]) - (phi2[DIR_PMM] - phi2[DIR_MPP])) + ((phi2[DIR_PMP] - phi2[DIR_MPM]) - (phi2[DIR_PPM] - phi2[DIR_MMP])))
-		+ WEIGTH[NE] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
-		+WEIGTH[N] * (phi2[T] - phi2[B]));
+		+ WEIGTH[DIR_PP0] * (((phi2[DIR_P0P] - phi2[DIR_M0M]) - (phi2[DIR_P0M] - phi2[DIR_M0P])) + ((phi2[DIR_0MP] - phi2[DIR_0PM]) + (phi2[DIR_0PP] - phi2[DIR_0MM])))) +
+		+WEIGTH[DIR_0P0] * (phi2[DIR_00P] - phi2[DIR_00M]));
 	//LBMReal sum = 0.0;
 	//for (int k = FSTARTDIR; k <= FENDDIR; k++) {
 	//    sum += WEIGTH[k] * DX3[k] * phi2[k];
@@ -3369,12 +3369,12 @@ LBMReal MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::nabla2_phi()
 	sum += WEIGTH[DIR_0PP] * (
 			(((phi[DIR_0PP] - phi[DIR_000]) + (phi[DIR_0MM] - phi[DIR_000])) + ((phi[DIR_0MP] - phi[DIR_000]) + (phi[DIR_0PM] - phi[DIR_000])))
 		+	(((phi[DIR_P0P] - phi[DIR_000]) + (phi[DIR_M0M] - phi[DIR_000])) + ((phi[DIR_M0P] - phi[DIR_000]) + (phi[DIR_P0M] - phi[DIR_000])))
-		+	(((phi[NE] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
+		+	(((phi[DIR_PP0] - phi[DIR_000]) + (phi[DIR_MM0] - phi[DIR_000])) + ((phi[DIR_MP0] - phi[DIR_000]) + (phi[DIR_PM0] - phi[DIR_000])))
 		);
-	sum += WEIGTH[T] * (
-			((phi[T] - phi[DIR_000]) + (phi[B] - phi[DIR_000]))
-		+	((phi[N] - phi[DIR_000]) + (phi[S] - phi[DIR_000]))
-		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[W] - phi[DIR_000]))
+	sum += WEIGTH[DIR_00P] * (
+			((phi[DIR_00P] - phi[DIR_000]) + (phi[DIR_00M] - phi[DIR_000]))
+		+	((phi[DIR_0P0] - phi[DIR_000]) + (phi[DIR_0M0] - phi[DIR_000]))
+		+	((phi[DIR_P00] - phi[DIR_000]) + (phi[DIR_M00] - phi[DIR_000]))
 		);
     //for (int k = FSTARTDIR; k <= FENDDIR; k++) {
     //    sum += WEIGTH[k] * (phi[k] - phi[REST]);
@@ -3405,9 +3405,9 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::computePhasefield()
                     int x3p = x3 + 1;
 
                     h[DIR_P00]   = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
-                    h[N]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
-                    h[T]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
-                    h[NE]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+                    h[DIR_0P0]   = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+                    h[DIR_00P]   = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+                    h[DIR_PP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
                     h[DIR_MP0]  = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
                     h[DIR_P0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
                     h[DIR_M0P]  = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
@@ -3418,9 +3418,9 @@ void MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel::computePhasefield()
                     h[DIR_PMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
                     h[DIR_MMP] = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
 
-                    h[W]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
-                    h[S]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
-                    h[B]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+                    h[DIR_M00]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+                    h[DIR_0M0]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+                    h[DIR_00M]   = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
                     h[DIR_MM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
                     h[DIR_PM0]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
                     h[DIR_M0M]  = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
index 30b9c875a..09cd40c8e 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyInterpolationProcessor.cpp
@@ -121,11 +121,11 @@ void RheologyInterpolationProcessor::calcMoments(const LBMReal* const f, LBMReal
 
    press = rho; //interpolate rho!
 
-   kxy   = -3.*omega*((((f[DIR_MMP]+f[DIR_PPM])-(f[DIR_MPP]+f[DIR_PMM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_MPM]+f[DIR_PMP])))+((f[DIR_MM0]+f[NE])-(f[DIR_MP0]+f[DIR_PM0]))-(vx1*vx2));// might not be optimal MG 25.2.13
+   kxy   = -3.*omega*((((f[DIR_MMP]+f[DIR_PPM])-(f[DIR_MPP]+f[DIR_PMM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_MPM]+f[DIR_PMP])))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_MP0]+f[DIR_PM0]))-(vx1*vx2));// might not be optimal MG 25.2.13
    kyz   = -3.*omega*((((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMP]+f[DIR_MPM]))+((f[DIR_PMM]+f[DIR_MPP])-(f[DIR_MMP]+f[DIR_PPM])))+((f[DIR_0MM]+f[DIR_0PP])-(f[DIR_0MP]+f[DIR_0PM]))-(vx2*vx3));
    kxz   = -3.*omega*((((f[DIR_MPM]+f[DIR_PMP])-(f[DIR_MMP]+f[DIR_PPM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMM]+f[DIR_MPP])))+((f[DIR_M0M]+f[DIR_P0P])-(f[DIR_M0P]+f[DIR_P0M]))-(vx1*vx3));
-   kxxMyy = -3./2.*omega*((((f[D3Q27System::DIR_M0M]+f[DIR_P0P])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_M0P]+f[DIR_P0M])-(f[DIR_0MP]+f[DIR_0PM])))+((f[W]+f[DIR_P00])-(f[S]+f[N]))-(vx1*vx1-vx2*vx2));
-   kxxMzz = -3./2.*omega*((((f[DIR_MP0]+f[DIR_PM0])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_MM0]+f[NE])-(f[DIR_0MP]+f[DIR_0PM])))+((f[W]+f[DIR_P00])-(f[B]+f[T]))-(vx1*vx1-vx3*vx3));
+   kxxMyy = -3./2.*omega*((((f[D3Q27System::DIR_M0M]+f[DIR_P0P])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_M0P]+f[DIR_P0M])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_0M0]+f[DIR_0P0]))-(vx1*vx1-vx2*vx2));
+   kxxMzz = -3./2.*omega*((((f[DIR_MP0]+f[DIR_PM0])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_00M]+f[DIR_00P]))-(vx1*vx1-vx3*vx3));
 }
 //////////////////////////////////////////////////////////////////////////
 void RheologyInterpolationProcessor::calcInterpolatedCoefficiets_intern(const D3Q27ICell& icell,
@@ -444,12 +444,12 @@ void RheologyInterpolationProcessor::calcInterpolatedNode(LBMReal* f, /*LBMReal
    D3Q27System::calcIncompFeq(feq,rho,vx1,vx2,vx3);
 
    f[DIR_P00]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[DIR_P00];
-   f[W]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[W];
-   f[N]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[N];
-   f[S]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[S];
-   f[T]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[T];
-   f[B]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[B];
-   f[NE]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[NE];
+   f[DIR_M00]    = f_E    + xs*x_E    + ys*y_E    + zs*z_E    + xs*ys*xy_E    + xs*zs*xz_E    + ys*zs*yz_E    + feq[DIR_M00];
+   f[DIR_0P0]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[DIR_0P0];
+   f[DIR_0M0]    = f_N    + xs*x_N    + ys*y_N    + zs*z_N    + xs*ys*xy_N    + xs*zs*xz_N    + ys*zs*yz_N    + feq[DIR_0M0];
+   f[DIR_00P]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[DIR_00P];
+   f[DIR_00M]    = f_T    + xs*x_T    + ys*y_T    + zs*z_T    + xs*ys*xy_T    + xs*zs*xz_T    + ys*zs*yz_T    + feq[DIR_00M];
+   f[DIR_PP0]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[DIR_PP0];
    f[DIR_MM0]   = f_NE   + xs*x_NE   + ys*y_NE   + zs*z_NE   + xs*ys*xy_NE   + xs*zs*xz_NE   + ys*zs*yz_NE   + feq[DIR_MM0];
    f[DIR_PM0]   = f_SE   + xs*x_SE   + ys*y_SE   + zs*z_SE   + xs*ys*xy_SE   + xs*zs*xz_SE   + ys*zs*yz_SE   + feq[DIR_PM0];
    f[DIR_MP0]   = f_SE   + xs*x_SE   + ys*y_SE   + zs*z_SE   + xs*ys*xy_SE   + xs*zs*xz_SE   + ys*zs*yz_SE   + feq[DIR_MP0];
@@ -633,12 +633,12 @@ void RheologyInterpolationProcessor::calcInterpolatedNodeFC(LBMReal* f, LBMReal
    f_TNW = eps_new*((ay + az + bx - bz + cx - cy+kxyAverage+kxzAverage-kyzAverage)/(72.*o));
 
    f[DIR_P00]    = f_E    + feq[DIR_P00];
-   f[W]    = f_E    + feq[W];
-   f[N]    = f_N    + feq[N];
-   f[S]    = f_N    + feq[S];
-   f[T]    = f_T    + feq[T];
-   f[B]    = f_T    + feq[B];
-   f[NE]   = f_NE   + feq[NE];
+   f[DIR_M00]    = f_E    + feq[DIR_M00];
+   f[DIR_0P0]    = f_N    + feq[DIR_0P0];
+   f[DIR_0M0]    = f_N    + feq[DIR_0M0];
+   f[DIR_00P]    = f_T    + feq[DIR_00P];
+   f[DIR_00M]    = f_T    + feq[DIR_00M];
+   f[DIR_PP0]   = f_NE   + feq[DIR_PP0];
    f[DIR_MM0]   = f_NE   + feq[DIR_MM0];
    f[DIR_PM0]   = f_SE   + feq[DIR_PM0];
    f[DIR_MP0]   = f_SE   + feq[DIR_MP0];
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index 358c61516..0ba49c1a0 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -248,12 +248,12 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
 
 
                f[DIR_P00]    = f_E    + feq[DIR_P00];
-               f[W]    = f_E    + feq[W];
-               f[N]    = f_N    + feq[N];
-               f[S]    = f_N    + feq[S];
-               f[T]    = f_T    + feq[T];
-               f[B]    = f_T    + feq[B];
-               f[NE]   = f_NE   + feq[NE];
+               f[DIR_M00]    = f_E    + feq[DIR_M00];
+               f[DIR_0P0]    = f_N    + feq[DIR_0P0];
+               f[DIR_0M0]    = f_N    + feq[DIR_0M0];
+               f[DIR_00P]    = f_T    + feq[DIR_00P];
+               f[DIR_00M]    = f_T    + feq[DIR_00M];
+               f[DIR_PP0]   = f_NE   + feq[DIR_PP0];
                f[DIR_MM0]   = f_NE   + feq[DIR_MM0];
                f[DIR_PM0]   = f_SE   + feq[DIR_PM0];
                f[DIR_MP0]   = f_SE   + feq[DIR_MP0];
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
index 451ccfc1f..8cd256a9b 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MetisPartitioningGridVisitor.cpp
@@ -256,9 +256,9 @@ void MetisPartitioningGridVisitor::clear()
 int MetisPartitioningGridVisitor::getEdgeWeight(int dir)
 {
     using namespace D3Q27System;
-    if (dir <= B) {
+    if (dir <= DIR_00M) {
         return 100;
-    } else if (dir >= NE && dir <= DIR_0MP) {
+    } else if (dir >= DIR_PP0 && dir <= DIR_0MP) {
         return 10;
     } else if (dir >= DIR_PPP) {
         return 1;
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
index e321725b6..a35fc289b 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseInitDistributionsBlockVisitor.cpp
@@ -228,12 +228,12 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 
 
 					f[DIR_P00]    =  geq[DIR_P00]    ;
-					f[W]    =  geq[W]    ;
-					f[N]    =  geq[N]    ;
-					f[S]    =  geq[S]    ;
-					f[T]    =  geq[T]    ;
-					f[B]    =  geq[B]    ;
-					f[NE]   =  geq[NE]   ;
+					f[DIR_M00]    =  geq[DIR_M00]    ;
+					f[DIR_0P0]    =  geq[DIR_0P0]    ;
+					f[DIR_0M0]    =  geq[DIR_0M0]    ;
+					f[DIR_00P]    =  geq[DIR_00P]    ;
+					f[DIR_00M]    =  geq[DIR_00M]    ;
+					f[DIR_PP0]   =  geq[DIR_PP0]   ;
 					f[DIR_MM0]   =  geq[DIR_MM0]   ;
 					f[DIR_PM0]   =  geq[DIR_PM0]   ;
 					f[DIR_MP0]   =  geq[DIR_MP0]   ;
@@ -259,12 +259,12 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 					distributionsF->setDistributionInv(f, ix1, ix2, ix3);
 
 					f[DIR_P00]    =  phi * feq[DIR_P00]    / rho;
-					f[W]    =  phi * feq[W]    / rho;
-					f[N]    =  phi * feq[N]    / rho;
-					f[S]    =  phi * feq[S]    / rho;
-					f[T]    =  phi * feq[T]    / rho;
-					f[B]    =  phi * feq[B]    / rho;
-					f[NE]   =  phi * feq[NE]   / rho;
+					f[DIR_M00]    =  phi * feq[DIR_M00]    / rho;
+					f[DIR_0P0]    =  phi * feq[DIR_0P0]    / rho;
+					f[DIR_0M0]    =  phi * feq[DIR_0M0]    / rho;
+					f[DIR_00P]    =  phi * feq[DIR_00P]    / rho;
+					f[DIR_00M]    =  phi * feq[DIR_00M]    / rho;
+					f[DIR_PP0]   =  phi * feq[DIR_PP0]   / rho;
 					f[DIR_MM0]   =  phi * feq[DIR_MM0]   / rho;
 					f[DIR_PM0]   =  phi * feq[DIR_PM0]   / rho;
 					f[DIR_MP0]   =  phi * feq[DIR_MP0]   / rho;
@@ -292,12 +292,12 @@ void MultiphaseInitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPt
 					if (distributionsH2) {
 
 						f[DIR_P00]    = (1.-phi) * feq[DIR_P00] / rho;
-						f[W]    = (1.-phi) * feq[W] / rho;
-						f[N]    = (1.-phi) * feq[N] / rho;
-						f[S]    = (1.-phi) * feq[S] / rho;
-						f[T]    = (1.-phi) * feq[T] / rho;
-						f[B]    = (1.-phi) * feq[B] / rho;
-						f[NE]   = (1.-phi) * feq[NE] / rho;
+						f[DIR_M00]    = (1.-phi) * feq[DIR_M00] / rho;
+						f[DIR_0P0]    = (1.-phi) * feq[DIR_0P0] / rho;
+						f[DIR_0M0]    = (1.-phi) * feq[DIR_0M0] / rho;
+						f[DIR_00P]    = (1.-phi) * feq[DIR_00P] / rho;
+						f[DIR_00M]    = (1.-phi) * feq[DIR_00M] / rho;
+						f[DIR_PP0]   = (1.-phi) * feq[DIR_PP0] / rho;
 						f[DIR_MM0]   = (1.-phi) * feq[DIR_MM0] / rho;
 						f[DIR_PM0]   = (1.-phi) * feq[DIR_PM0] / rho;
 						f[DIR_MP0]   = (1.-phi) * feq[DIR_MP0] / rho;
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
index dd9fb72af..623fcade6 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseVelocityFormInitDistributionsBlockVisitor.cpp
@@ -220,12 +220,12 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 
 
 					f[DIR_P00]    =  geq[DIR_P00]    ;
-					f[W]    =  geq[W]    ;
-					f[N]    =  geq[N]    ;
-					f[S]    =  geq[S]    ;
-					f[T]    =  geq[T]    ;
-					f[B]    =  geq[B]    ;
-					f[NE]   =  geq[NE]   ;
+					f[DIR_00M]    =  geq[DIR_00M]    ;
+					f[DIR_0P0]    =  geq[DIR_0P0]    ;
+					f[DIR_0M0]    =  geq[DIR_0M0]    ;
+					f[DIR_00P]    =  geq[DIR_00P]    ;
+					f[DIR_00M]    =  geq[DIR_00M]    ;
+					f[DIR_PP0]   =  geq[DIR_PP0]   ;
 					f[DIR_MM0]   =  geq[DIR_MM0]   ;
 					f[DIR_PM0]   =  geq[DIR_PM0]   ;
 					f[DIR_MP0]   =  geq[DIR_MP0]   ;
@@ -251,12 +251,12 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 					distributionsF->setDistributionInv(f, ix1, ix2, ix3);
 
 					f[DIR_P00]    =  phi * feq[DIR_P00]    ;// / rho;
-					f[W]    =  phi * feq[W]    ;// / rho;
-					f[N]    =  phi * feq[N]    ;// / rho;
-					f[S]    =  phi * feq[S]    ;// / rho;
-					f[T]    =  phi * feq[T]    ;// / rho;
-					f[B]    =  phi * feq[B]    ;// / rho;
-					f[NE]   =  phi * feq[NE]   ;// / rho;
+					f[DIR_M00]    =  phi * feq[DIR_M00]    ;// / rho;
+					f[DIR_0P0]    =  phi * feq[DIR_0P0]    ;// / rho;
+					f[DIR_0M0]    =  phi * feq[DIR_0M0]    ;// / rho;
+					f[DIR_00P]    =  phi * feq[DIR_00P]    ;// / rho;
+					f[DIR_00M]    =  phi * feq[DIR_00M]    ;// / rho;
+					f[DIR_PP0]   =  phi * feq[DIR_PP0]   ;// / rho;
 					f[DIR_MM0]   =  phi * feq[DIR_MM0]   ;// / rho;
 					f[DIR_PM0]   =  phi * feq[DIR_PM0]   ;// / rho;
 					f[DIR_MP0]   =  phi * feq[DIR_MP0]   ;// / rho;
@@ -284,12 +284,12 @@ void MultiphaseVelocityFormInitDistributionsBlockVisitor::visit(const SPtr<Grid3
 					if (distributionsH2) {
 
 						f[DIR_P00]    = (1.-phi) * feq[DIR_P00]   ;// / rho;
-						f[W]    = (1.-phi) * feq[W]   ;// / rho;
-						f[N]    = (1.-phi) * feq[N]   ;// / rho;
-						f[S]    = (1.-phi) * feq[S]   ;// / rho;
-						f[T]    = (1.-phi) * feq[T]   ;// / rho;
-						f[B]    = (1.-phi) * feq[B]   ;// / rho;
-						f[NE]   = (1.-phi) * feq[NE]  ;// / rho;
+						f[DIR_M00]    = (1.-phi) * feq[DIR_M00]   ;// / rho;
+						f[DIR_0P0]    = (1.-phi) * feq[DIR_0P0]   ;// / rho;
+						f[DIR_0M0]    = (1.-phi) * feq[DIR_0M0]   ;// / rho;
+						f[DIR_00P]    = (1.-phi) * feq[DIR_00P]   ;// / rho;
+						f[DIR_00M]    = (1.-phi) * feq[DIR_00M]   ;// / rho;
+						f[DIR_PP0]   = (1.-phi) * feq[DIR_PP0]  ;// / rho;
 						f[DIR_MM0]   = (1.-phi) * feq[DIR_MM0]  ;// / rho;
 						f[DIR_PM0]   = (1.-phi) * feq[DIR_PM0]  ;// / rho;
 						f[DIR_MP0]   = (1.-phi) * feq[DIR_MP0]  ;// / rho;
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
index 07aec29ef..689d84d07 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetInterpolationDirsBlockVisitor.cpp
@@ -32,103 +32,103 @@ void SetInterpolationDirsBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> bl
             if (p_nblock) {
                 bool flagDir;
                 switch (dir) {
-                    case NE:
-                        checkFlagDir(grid, E, N, flagDir, ix1, ix2, ix3, level);
+                    case DIR_PP0:
+                        checkFlagDir(grid, DIR_P00, DIR_0P0, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MM0:
-                        checkFlagDir(grid, W, S, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0M0, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_PM0:
-                        checkFlagDir(grid, E, S, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_0M0, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MP0:
-                        checkFlagDir(grid, W, N, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0P0, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_P0P:
-                        checkFlagDir(grid, E, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_M0M:
-                        checkFlagDir(grid, W, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_P0M:
-                        checkFlagDir(grid, E, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_M0P:
-                        checkFlagDir(grid, W, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_0PP:
-                        checkFlagDir(grid, N, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_0P0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_0MM:
-                        checkFlagDir(grid, S, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_0M0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_0PM:
-                        checkFlagDir(grid, N, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_0P0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_0MP:
-                        checkFlagDir(grid, S, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_0M0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_PPP:
-                        checkFlagDir(grid, E, N, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_0P0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MMP:
-                        checkFlagDir(grid, W, S, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0M0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_PMP:
-                        checkFlagDir(grid, E, S, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_0M0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MPP:
-                        checkFlagDir(grid, W, N, T, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0P0, DIR_00P, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_PPM:
-                        checkFlagDir(grid, E, N, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_0P0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MMM:
-                        checkFlagDir(grid, W, S, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0M0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_PMM:
-                        checkFlagDir(grid, E, S, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_P00, DIR_0M0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
                     case DIR_MPM:
-                        checkFlagDir(grid, W, N, B, flagDir, ix1, ix2, ix3, level);
+                        checkFlagDir(grid, DIR_M00, DIR_0P0, DIR_00M, flagDir, ix1, ix2, ix3, level);
                         if (!flagDir)
                             continue;
                         break;
-- 
GitLab