diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
index 006f45619d46051cea296584e48440acdddd7291..3c96343077f1030f7e6477bb352635c870ece4ce 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
@@ -64,15 +64,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
    switch (direction)
    {
    case E:
-      f[E] = ftemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[E];
-      f[NE] = ftemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[NE];
-      f[SE] = ftemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[SE];
-      f[TE] = ftemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TE];
-      f[BE] = ftemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BE];
-      f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TNE];
-      f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TSE];
-      f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BNE];
-      f[BSE] = ftemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BSE];
+      f[E] = ftemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[E];
+      f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+      f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[SE];
+      f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TE];
+      f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BE];
+      f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TNE];
+      f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TSE];
+      f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BNE];
+      f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BSE];
 
       distributions->setDistributionInvForDirection(f[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
       distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -85,15 +85,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
       break;
    case W:
-      f[W] = ftemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[W];
-      f[NW] = ftemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[NW];
-      f[SW] = ftemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[SW];
-      f[TW] = ftemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TW];
-      f[BW] = ftemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BW];
-      f[TNW] = ftemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TNW];
-      f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TSW];
-      f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BNW];
-      f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BSW];
+      f[W] = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+      f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[NW];
+      f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[SW];
+      f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TW];
+      f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BW];
+      f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TNW];
+      f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TSW];
+      f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BNW];
+      f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BSW];
 
       distributions->setDistributionInvForDirection(f[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
       distributions->setDistributionInvForDirection(f[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -106,15 +106,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
       break;
    case N:
-      f[N] = ftemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[N];
-      f[NE] = ftemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NE];
-      f[NW] = ftemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NW];
-      f[TN] = ftemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TN];
-      f[BN] = ftemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BN];
-      f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNE];
-      f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNW];
-      f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNE];
-      f[BNW] = ftemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNW];
+      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];
+      f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NW];
+      f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TN];
+      f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BN];
+      f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNE];
+      f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNW];
+      f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNE];
+      f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNW];
 
       distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
       distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -127,15 +127,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
       break;
    case S:
-      f[S] = ftemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[S];
-      f[SE] = ftemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SE];
-      f[SW] = ftemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SW];
-      f[TS] = ftemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TS];
-      f[BS] = ftemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BS];
-      f[TSE] = ftemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSE];
-      f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSW];
-      f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSE];
-      f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSW];
+      f[S] = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+      f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SE];
+      f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SW];
+      f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TS];
+      f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BS];
+      f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSE];
+      f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSW];
+      f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSE];
+      f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSW];
 
       distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
       distributions->setDistributionInvForDirection(f[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -148,15 +148,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
       break;
    case T:
-      f[T] = ftemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[T];
-      f[TE] = ftemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TE];
-      f[TW] = ftemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TW];
-      f[TN] = ftemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TN];
-      f[TS] = ftemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TS];
-      f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNE];
-      f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNW];
-      f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSE];
-      f[TSW] = ftemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSW];
+      f[T] = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+      f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TE];
+      f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TW];
+      f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TN];
+      f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TS];
+      f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNE];
+      f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNW];
+      f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSE];
+      f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSW];
 
       distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
       distributions->setDistributionInvForDirection(f[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -169,15 +169,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributions->setDistributionInvForDirection(f[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
       break;
    case B:
-      f[B] = ftemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[B];
-      f[BE] = ftemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BE];
-      f[BW] = ftemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BW];
-      f[BN] = ftemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BN];
-      f[BS] = ftemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BS];
-      f[BNE] = ftemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNE];
-      f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNW];
-      f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSE];
-      f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSW];
+      f[B] = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+      f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BE];
+      f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BW];
+      f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BN];
+      f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BS];
+      f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNE];
+      f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNW];
+      f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSE];
+      f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSW];
 
       distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
       distributions->setDistributionInvForDirection(f[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
@@ -208,15 +208,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
    switch (direction)
    {
    case E:
-      h[E]  = htemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[E];
-      h[NE] = htemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[NE];
-      h[SE] = htemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[SE];
-      h[TE] = htemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TE];
-      h[BE] = htemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BE];
-      h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TNE];
-      h[TSE] = htemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TSE];
-      h[BNE] = htemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BNE];
-      h[BSE] = htemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BSE];
+      h[E]  = htemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[E];
+      h[NE] = htemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[NE];
+      h[SE] = htemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[SE];
+      h[TE] = htemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TE];
+      h[BE] = htemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BE];
+      h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TNE];
+      h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TSE];
+      h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BNE];
+      h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BSE];
 
       distributionsH->setDistributionInvForDirection(h[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
       distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -229,15 +229,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
       break;
    case W:
-      h[W] = htemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[W];
-      h[NW] = htemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[NW];
-      h[SW] = htemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[SW];
-      h[TW] = htemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TW];
-      h[BW] = htemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BW];
-      h[TNW] = htemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TNW];
-      h[TSW] = htemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TSW];
-      h[BNW] = htemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BNW];
-      h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BSW];
+      h[W] = htemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[W];
+      h[NW] = htemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[NW];
+      h[SW] = htemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[SW];
+      h[TW] = htemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TW];
+      h[BW] = htemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BW];
+      h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TNW];
+      h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TSW];
+      h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BNW];
+      h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BSW];
 
       distributionsH->setDistributionInvForDirection(h[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
       distributionsH->setDistributionInvForDirection(h[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -250,15 +250,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
       break;
    case N:
-      h[N] = htemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[N];
-      h[NE] = htemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[NE];
-      h[NW] = htemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[NW];
-      h[TN] = htemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TN];
-      h[BN] = htemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BN];
-      h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TNE];
-      h[TNW] = htemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TNW];
-      h[BNE] = htemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BNE];
-      h[BNW] = htemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BNW];
+      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];
+      h[NW] = htemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[NW];
+      h[TN] = htemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TN];
+      h[BN] = htemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BN];
+      h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TNE];
+      h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TNW];
+      h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BNE];
+      h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BNW];
 
       distributionsH->setDistributionInvForDirection(h[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
       distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -271,15 +271,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
       break;
    case S:
-      h[S] = htemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[S];
-      h[SE] = htemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[SE];
-      h[SW] = htemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[SW];
-      h[TS] = htemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TS];
-      h[BS] = htemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BS];
-      h[TSE] = htemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TSE];
-      h[TSW] = htemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TSW];
-      h[BSE] = htemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BSE];
-      h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BSW];
+      h[S] = htemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[S];
+      h[SE] = htemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[SE];
+      h[SW] = htemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[SW];
+      h[TS] = htemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TS];
+      h[BS] = htemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BS];
+      h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TSE];
+      h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TSW];
+      h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BSE];
+      h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BSW];
 
       distributionsH->setDistributionInvForDirection(h[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
       distributionsH->setDistributionInvForDirection(h[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -292,15 +292,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
       break;
    case T:
-      h[T] = htemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[T];
-      h[TE] = htemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TE];
-      h[TW] = htemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TW];
-      h[TN] = htemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TN];
-      h[TS] = htemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TS];
-      h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TNE];
-      h[TNW] = htemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TNW];
-      h[TSE] = htemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TSE];
-      h[TSW] = htemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TSW];
+      h[T] = htemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[T];
+      h[TE] = htemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TE];
+      h[TW] = htemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TW];
+      h[TN] = htemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TN];
+      h[TS] = htemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TS];
+      h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TNE];
+      h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TNW];
+      h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TSE];
+      h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TSW];
 
       distributionsH->setDistributionInvForDirection(h[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
       distributionsH->setDistributionInvForDirection(h[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -313,15 +313,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
       distributionsH->setDistributionInvForDirection(h[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
       break;
    case B:
-      h[B] = htemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[B];
-      h[BE] = htemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BE];
-      h[BW] = htemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BW];
-      h[BN] = htemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BN];
-      h[BS] = htemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BS];
-      h[BNE] = htemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BNE];
-      h[BNW] = htemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BNW];
-      h[BSE] = htemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BSE];
-      h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BSW];
+      h[B] = htemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[B];
+      h[BE] = htemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BE];
+      h[BW] = htemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BW];
+      h[BN] = htemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BN];
+      h[BS] = htemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BS];
+      h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BNE];
+      h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BNW];
+      h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BSE];
+      h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BSW];
 
       distributionsH->setDistributionInvForDirection(h[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
       distributionsH->setDistributionInvForDirection(h[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
index 0b5884d12229cfccab6647bd0ec52add2e238bbc..0a4280c75832c4eb6c71a62682f346aed4fba64f 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
@@ -66,15 +66,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
 
     switch (direction) {
         case E:
-            f[E]   = ftemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[E];
-            f[NE]  = ftemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[NE];
-            f[SE]  = ftemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[SE];
-            f[TE]  = ftemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TE];
-            f[BE]  = ftemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BE];
-            f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TNE];
-            f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TSE];
-            f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BNE];
-            f[BSE] = ftemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BSE];
+            f[E]   = ftemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[E];
+            f[NE]  = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+            f[SE]  = ftemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[SE];
+            f[TE]  = ftemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TE];
+            f[BE]  = ftemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BE];
+            f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TNE];
+            f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TSE];
+            f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BNE];
+            f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BSE];
 
             distributions->setDistributionInvForDirection(f[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
             distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -87,15 +87,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
             break;
         case W:
-            f[W]   = ftemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[W];
-            f[NW]  = ftemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[NW];
-            f[SW]  = ftemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[SW];
-            f[TW]  = ftemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TW];
-            f[BW]  = ftemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BW];
-            f[TNW] = ftemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TNW];
-            f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TSW];
-            f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BNW];
-            f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BSW];
+            f[W]   = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+            f[NW]  = ftemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[NW];
+            f[SW]  = ftemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[SW];
+            f[TW]  = ftemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TW];
+            f[BW]  = ftemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BW];
+            f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TNW];
+            f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TSW];
+            f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BNW];
+            f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BSW];
 
             distributions->setDistributionInvForDirection(f[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
             distributions->setDistributionInvForDirection(f[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -108,15 +108,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
             break;
         case N:
-            f[N]   = ftemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[N];
-            f[NE]  = ftemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NE];
-            f[NW]  = ftemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NW];
-            f[TN]  = ftemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TN];
-            f[BN]  = ftemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BN];
-            f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNE];
-            f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNW];
-            f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNE];
-            f[BNW] = ftemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNW];
+            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];
+            f[NW]  = ftemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NW];
+            f[TN]  = ftemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TN];
+            f[BN]  = ftemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BN];
+            f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNE];
+            f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNW];
+            f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNE];
+            f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNW];
 
             distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
             distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -129,15 +129,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
             break;
         case S:
-            f[S]   = ftemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[S];
-            f[SE]  = ftemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SE];
-            f[SW]  = ftemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SW];
-            f[TS]  = ftemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TS];
-            f[BS]  = ftemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BS];
-            f[TSE] = ftemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSE];
-            f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSW];
-            f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSE];
-            f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSW];
+            f[S]   = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+            f[SE]  = ftemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SE];
+            f[SW]  = ftemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SW];
+            f[TS]  = ftemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TS];
+            f[BS]  = ftemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BS];
+            f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSE];
+            f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSW];
+            f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSE];
+            f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSW];
 
             distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
             distributions->setDistributionInvForDirection(f[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -150,15 +150,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
             break;
         case T:
-            f[T]   = ftemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[T];
-            f[TE]  = ftemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TE];
-            f[TW]  = ftemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TW];
-            f[TN]  = ftemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TN];
-            f[TS]  = ftemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TS];
-            f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNE];
-            f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNW];
-            f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSE];
-            f[TSW] = ftemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSW];
+            f[T]   = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+            f[TE]  = ftemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TE];
+            f[TW]  = ftemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TW];
+            f[TN]  = ftemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TN];
+            f[TS]  = ftemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TS];
+            f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNE];
+            f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNW];
+            f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSE];
+            f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSW];
 
             distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
             distributions->setDistributionInvForDirection(f[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -171,15 +171,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
             distributions->setDistributionInvForDirection(f[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
             break;
         case B:
-            f[B]   = ftemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[B];
-            f[BE]  = ftemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BE];
-            f[BW]  = ftemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BW];
-            f[BN]  = ftemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BN];
-            f[BS]  = ftemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BS];
-            f[BNE] = ftemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNE];
-            f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNW];
-            f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSE];
-            f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSW];
+            f[B]   = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+            f[BE]  = ftemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BE];
+            f[BW]  = ftemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BW];
+            f[BN]  = ftemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BN];
+            f[BS]  = ftemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BS];
+            f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNE];
+            f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNW];
+            f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSE];
+            f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSW];
 
             distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
             distributions->setDistributionInvForDirection(f[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
index 72e7114e6d13504727ea8d86d7641627df6056f6..af7c6d3930f2d6f5c6ef17ddc046ba5a4f4a2318 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
@@ -217,7 +217,7 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 					//LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
 					//LBMReal omega = Thixotropy::getPowellEyringCollFactor(collFactor, shearRate, rho);
 					LBMReal omega = Thixotropy::getBinghamCollFactor(collFactor, shearRate, rho);
-					LBMReal viscosity = (omega == 0) ? 0 : c1o3 * (c1/omega-c1o2);
+					LBMReal viscosity = (omega == 0) ? 0 : UbMath::c1o3 * (UbMath::c1/omega-UbMath::c1o2);
 
 					
 					data[index++].push_back(viscosity);
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
index 57e26c7e0f4d33ad2cf9983daab85dca37d8ed41..f0955c0b1fd92e4c7fb34ba463c4ddb05d3fa886 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
@@ -8,6 +8,8 @@
 
 #define PROOF_CORRECTNESS
 
+using namespace UbMath;
+
 //////////////////////////////////////////////////////////////////////////
 CumulantLBMKernel::CumulantLBMKernel()
 {
@@ -60,7 +62,7 @@ SPtr<LBMKernel> CumulantLBMKernel::clone()
    }
    else
    {
-      dynamicPointerCast<CumulantLBMKernel>(kernel)->OxxPyyPzz = one;
+      dynamicPointerCast<CumulantLBMKernel>(kernel)->OxxPyyPzz = UbMath::one;
    }
    return kernel;
 }
@@ -1139,7 +1141,7 @@ void CumulantLBMKernel::nodeCollision(int step, int x1, int x2, int x3)
       (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
       ((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) + mfbbb;
 
-   LBMReal rho = one + drho;
+   LBMReal rho = UbMath::one + drho;
    ////////////////////////////////////////////////////////////////////////////////////
    LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
       (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index fc252e1800e62bff6c31f33457898bf72cbd73ed..a0068338c72d4779ee2be95f143a1dea0379c160 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -67,11 +67,6 @@ class for global system-functions
 usage: ...
 */
 
-
-#ifndef SWIG
-   using namespace UbMath;
-#endif
-
    namespace D3Q27System
    {
       //enum COLLISIONMODEL { UNDEFINED, INCOMPLBGKMODEL,   COMPLBGKMODEL,   COMPLBGKWTMODEL,   INCOMPLBGKLESMODEL, INCOMPLBGKNONNEWTONIANMODEL    
@@ -300,7 +295,7 @@ usage: ...
       //ACHTUNG: gilt nicht fuer alle modelle -> praedikat verwenden anstelle static! toDo
       static LBMReal getPressure(const LBMReal* const& f/*[27]*/)
       {
-         return  REAL_CAST(c1o3) * getDensity(f);
+         return  REAL_CAST(UbMath::c1o3) * getDensity(f);
       }
       /*=====================================================================*/
       static LBMReal getIncompVelocityX1(const LBMReal* const& f/*[27]*/)
@@ -411,7 +406,7 @@ usage: ...
          D3Q27System::calcIncompVelocityX1(f, vx1);
          D3Q27System::calcIncompVelocityX2(f, vx2);
          D3Q27System::calcIncompVelocityX3(f, vx3);
-         LBMReal rho = drho + one;
+         LBMReal rho = drho + UbMath::one;
          vx1 /= rho;
          vx2 /= rho;
          vx3 /= rho;
@@ -423,68 +418,68 @@ usage: ...
 
          //switch(direction)    
          //{
-         //   case ZERO : return REAL_CAST( c8o27*rho*(1.0-cu_sq));
-         //   case E : return REAL_CAST(  c2o27*rho*(1.0+3.0*( vx1   )+c9o2*( vx1   )*( vx1   )-cu_sq));
-         //   case W : return REAL_CAST(  c2o27*rho*(1.0+3.0*(-vx1   )+c9o2*(-vx1   )*(-vx1   )-cu_sq));
-         //   case N : return REAL_CAST(  c2o27*rho*(1.0+3.0*(    vx2)+c9o2*(    vx2)*(    vx2)-cu_sq));
-         //   case S : return REAL_CAST(  c2o27*rho*(1.0+3.0*(   -vx2)+c9o2*(   -vx2)*(   -vx2)-cu_sq));
-         //   case T : return REAL_CAST(  c2o27*rho*(1.0+3.0*( vx3   )+c9o2*(    vx3)*(    vx3)-cu_sq));
-         //   case B : return REAL_CAST(  c2o27*rho*(1.0+3.0*(   -vx3)+c9o2*(   -vx3)*(   -vx3)-cu_sq));
-         //   case NE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq));
-         //   case SW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq));
-         //   case SE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1-vx2)+c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq));
-         //   case NW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq));
-         //   case TE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1+vx3)+c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq));
-         //   case BW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1-vx3)+c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq));
-         //   case BE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1-vx3)+c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq));
-         //   case TW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1+vx3)+c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq));
-         //   case TN : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx2+vx3)+c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq));
-         //   case BS : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx2-vx3)+c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq));
-         //   case BN : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx2-vx3)+c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq));
-         //   case TS : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx2+vx3)+c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq));
-         //   case TNE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq));
-         //   case BSW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq));
-         //   case BNE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq));
-         //   case TSW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq));
-         //   case TSE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq));
-         //   case BNW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq));
-         //   case BSE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq));
-         //   case TNW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
+         //   case ZERO : return REAL_CAST( UbMath::c8o27*rho*(1.0-cu_sq));
+         //   case E : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*( vx1   )+UbMath::c9o2*( vx1   )*( vx1   )-cu_sq));
+         //   case W : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*(-vx1   )+UbMath::c9o2*(-vx1   )*(-vx1   )-cu_sq));
+         //   case N : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*(    vx2)+UbMath::c9o2*(    vx2)*(    vx2)-cu_sq));
+         //   case S : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*(   -vx2)+UbMath::c9o2*(   -vx2)*(   -vx2)-cu_sq));
+         //   case T : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*( vx3   )+UbMath::c9o2*(    vx3)*(    vx3)-cu_sq));
+         //   case B : return REAL_CAST(  UbMath::c2o27*rho*(1.0+3.0*(   -vx3)+UbMath::c9o2*(   -vx3)*(   -vx3)-cu_sq));
+         //   case NE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1+vx2)+UbMath::c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq));
+         //   case SW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx2)+UbMath::c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq));
+         //   case SE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1-vx2)+UbMath::c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq));
+         //   case NW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx2)+UbMath::c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq));
+         //   case TE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1+vx3)+UbMath::c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq));
+         //   case BW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx3)+UbMath::c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq));
+         //   case BE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1-vx3)+UbMath::c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq));
+         //   case TW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx3)+UbMath::c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq));
+         //   case TN : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx2+vx3)+UbMath::c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq));
+         //   case BS : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx2-vx3)+UbMath::c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq));
+         //   case BN : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx2-vx3)+UbMath::c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq));
+         //   case TS : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx2+vx3)+UbMath::c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq));
+         //   case TNE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+UbMath::c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq));
+         //   case BSW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+UbMath::c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq));
+         //   case BNE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+UbMath::c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq));
+         //   case TSW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+UbMath::c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq));
+         //   case TSE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+UbMath::c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq));
+         //   case BNW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+UbMath::c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq));
+         //   case BSE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+UbMath::c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq));
+         //   case TNW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+UbMath::c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
          //   default: throw UbException(UB_EXARGS,"unknown dir");
          //}
 
 
          ////-----
-         LBMReal rho = drho + one;
+         LBMReal rho = drho + UbMath::one;
          switch (direction)
          {
-         case ZERO: return REAL_CAST(c8o27 * (drho + rho * (-cu_sq)));
-         case E: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq)));
-         case W: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq)));
-         case N: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq)));
-         case S: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq)));
-         case T: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq)));
-         case B: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq)));
-         case NE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq)));
-         case SW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq)));
-         case SE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq)));
-         case NW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq)));
-         case TE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq)));
-         case BW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq)));
-         case BE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq)));
-         case TW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq)));
-         case TN: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq)));
-         case BS: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq)));
-         case BN: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq)));
-         case TS: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq)));
-         case TNE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq)));
-         case BSW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq)));
-         case BNE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq)));
-         case TSW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq)));
-         case TSE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq)));
-         case BNW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq)));
-         case BSE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq)));
-         case TNW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq)));
+         case ZERO: return REAL_CAST(UbMath::c8o27 * (drho + rho * (-cu_sq)));
+         case E: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq)));
+         case W: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq)));
+         case N: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq)));
+         case S: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq)));
+         case T: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq)));
+         case B: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq)));
+         case NE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq)));
+         case SW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq)));
+         case SE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq)));
+         case NW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq)));
+         case TE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq)));
+         case BW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq)));
+         case BE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq)));
+         case TW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq)));
+         case TN: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq)));
+         case BS: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq)));
+         case BN: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq)));
+         case TS: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq)));
+         case TNE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq)));
+         case BSW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq)));
+         case BNE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq)));
+         case TSW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq)));
+         case TSE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq)));
+         case BNW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq)));
+         case BSE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq)));
+         case TNW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq)));
          default: throw UbException(UB_EXARGS, "unknown dir");
          }
 
@@ -494,66 +489,66 @@ usage: ...
       {
          //LBMReal cu_sq=1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
 
-         //feq[ZERO] =  c8o27*rho*(1.0-cu_sq);
-         //feq[E] =   c2o27*rho*(1.0+3.0*( vx1   )+c9o2*( vx1   )*( vx1   )-cu_sq);
-         //feq[W] =   c2o27*rho*(1.0+3.0*(-vx1   )+c9o2*(-vx1   )*(-vx1   )-cu_sq);
-         //feq[N] =   c2o27*rho*(1.0+3.0*(    vx2)+c9o2*(    vx2)*(    vx2)-cu_sq);
-         //feq[S] =   c2o27*rho*(1.0+3.0*(   -vx2)+c9o2*(   -vx2)*(   -vx2)-cu_sq);
-         //feq[T] =   c2o27*rho*(1.0+3.0*( vx3   )+c9o2*(    vx3)*(    vx3)-cu_sq);
-         //feq[B] =   c2o27*rho*(1.0+3.0*(   -vx3)+c9o2*(   -vx3)*(   -vx3)-cu_sq);
-         //feq[NE] =  c1o54*rho*(1.0+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
-         //feq[SW] =  c1o54*rho*(1.0+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
-         //feq[SE] =  c1o54*rho*(1.0+3.0*( vx1-vx2)+c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq);
-         //feq[NW] =  c1o54*rho*(1.0+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
-         //feq[TE] =  c1o54*rho*(1.0+3.0*( vx1+vx3)+c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq);
-         //feq[BW] =  c1o54*rho*(1.0+3.0*(-vx1-vx3)+c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq);
-         //feq[BE] =  c1o54*rho*(1.0+3.0*( vx1-vx3)+c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq);
-         //feq[TW] =  c1o54*rho*(1.0+3.0*(-vx1+vx3)+c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq);
-         //feq[TN] =  c1o54*rho*(1.0+3.0*( vx2+vx3)+c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq);
-         //feq[BS] =  c1o54*rho*(1.0+3.0*(-vx2-vx3)+c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq);
-         //feq[BN] =  c1o54*rho*(1.0+3.0*( vx2-vx3)+c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq);
-         //feq[TS] =  c1o54*rho*(1.0+3.0*(-vx2+vx3)+c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq);
-         //feq[TNE] = c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
-         //feq[BSW] = c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
-         //feq[BNE] = c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
-         //feq[TSW] = c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
-         //feq[TSE] = c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
-         //feq[BNW] = c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
-         //feq[BSE] = c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
-         //feq[TNW] = c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
+         //feq[ZERO] =  UbMath::c8o27*rho*(1.0-cu_sq);
+         //feq[E] =   UbMath::c2o27*rho*(1.0+3.0*( vx1   )+UbMath::c9o2*( vx1   )*( vx1   )-cu_sq);
+         //feq[W] =   UbMath::c2o27*rho*(1.0+3.0*(-vx1   )+UbMath::c9o2*(-vx1   )*(-vx1   )-cu_sq);
+         //feq[N] =   UbMath::c2o27*rho*(1.0+3.0*(    vx2)+UbMath::c9o2*(    vx2)*(    vx2)-cu_sq);
+         //feq[S] =   UbMath::c2o27*rho*(1.0+3.0*(   -vx2)+UbMath::c9o2*(   -vx2)*(   -vx2)-cu_sq);
+         //feq[T] =   UbMath::c2o27*rho*(1.0+3.0*( vx3   )+UbMath::c9o2*(    vx3)*(    vx3)-cu_sq);
+         //feq[B] =   UbMath::c2o27*rho*(1.0+3.0*(   -vx3)+UbMath::c9o2*(   -vx3)*(   -vx3)-cu_sq);
+         //feq[NE] =  UbMath::c1o54*rho*(1.0+3.0*( vx1+vx2)+UbMath::c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
+         //feq[SW] =  UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx2)+UbMath::c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
+         //feq[SE] =  UbMath::c1o54*rho*(1.0+3.0*( vx1-vx2)+UbMath::c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq);
+         //feq[NW] =  UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx2)+UbMath::c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
+         //feq[TE] =  UbMath::c1o54*rho*(1.0+3.0*( vx1+vx3)+UbMath::c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq);
+         //feq[BW] =  UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx3)+UbMath::c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq);
+         //feq[BE] =  UbMath::c1o54*rho*(1.0+3.0*( vx1-vx3)+UbMath::c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq);
+         //feq[TW] =  UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx3)+UbMath::c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq);
+         //feq[TN] =  UbMath::c1o54*rho*(1.0+3.0*( vx2+vx3)+UbMath::c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq);
+         //feq[BS] =  UbMath::c1o54*rho*(1.0+3.0*(-vx2-vx3)+UbMath::c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq);
+         //feq[BN] =  UbMath::c1o54*rho*(1.0+3.0*( vx2-vx3)+UbMath::c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq);
+         //feq[TS] =  UbMath::c1o54*rho*(1.0+3.0*(-vx2+vx3)+UbMath::c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq);
+         //feq[TNE] = UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+UbMath::c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
+         //feq[BSW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+UbMath::c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
+         //feq[BNE] = UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+UbMath::c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
+         //feq[TSW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+UbMath::c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
+         //feq[TSE] = UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+UbMath::c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
+         //feq[BNW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+UbMath::c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
+         //feq[BSE] = UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+UbMath::c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
+         //feq[TNW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+UbMath::c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
 
          //////////////////////////////////////////////////////////////////////////
 
          LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
-         LBMReal rho = drho + one;
+         LBMReal rho = drho + UbMath::one;
 
-         feq[ZERO] = c8o27 * (drho + rho * (-cu_sq));
-         feq[E] = c2o27 * (drho + rho * (3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq));
-         feq[W] = c2o27 * (drho + rho * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
-         feq[N] = c2o27 * (drho + rho * (3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq));
-         feq[S] = c2o27 * (drho + rho * (3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq));
-         feq[T] = c2o27 * (drho + rho * (3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq));
-         feq[B] = c2o27 * (drho + rho * (3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq));
-         feq[NE] = c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
-         feq[SW] = c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
-         feq[SE] = c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
-         feq[NW] = c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
-         feq[TE] = c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
-         feq[BW] = c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
-         feq[BE] = c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
-         feq[TW] = c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
-         feq[TN] = c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
-         feq[BS] = c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
-         feq[BN] = c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
-         feq[TS] = c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
-         feq[TNE] = c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
-         feq[BSW] = c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
-         feq[BNE] = c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
-         feq[TSW] = c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
-         feq[TSE] = c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
-         feq[BNW] = c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
-         feq[BSE] = c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
-         feq[TNW] = c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
+         feq[ZERO] = UbMath::c8o27 * (drho + rho * (-cu_sq));
+         feq[E] = UbMath::c2o27 * (drho + rho * (3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
+         feq[W] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq));
+         feq[N] = UbMath::c2o27 * (drho + rho * (3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq));
+         feq[S] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq));
+         feq[T] = UbMath::c2o27 * (drho + rho * (3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq));
+         feq[B] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq));
+         feq[NE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
+         feq[SW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
+         feq[SE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
+         feq[NW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
+         feq[TE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
+         feq[BW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
+         feq[BE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
+         feq[TW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
+         feq[TN] = UbMath::c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
+         feq[BS] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
+         feq[BN] = UbMath::c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
+         feq[TS] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
+         feq[TNE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
+         feq[BSW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
+         feq[BNE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
+         feq[TSW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
+         feq[TSE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
+         feq[BNW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
+         feq[BSE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
+         feq[TNW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
       }
       //////////////////////////////////////////////////////////////////////////
       static LBMReal getIncompFeqForDirection(const int& direction, const LBMReal& drho, const LBMReal& vx1, const LBMReal& vx2, const LBMReal& vx3)
@@ -562,33 +557,33 @@ usage: ...
 
          switch (direction)
          {
-         case ZERO: return REAL_CAST(c8o27 * (drho - cu_sq));
-         case E: return REAL_CAST(c2o27 * (drho + 3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq));
-         case W: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
-         case N: return REAL_CAST(c2o27 * (drho + 3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq));
-         case S: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq));
-         case T: return REAL_CAST(c2o27 * (drho + 3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq));
-         case B: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq));
-         case NE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
-         case SW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
-         case SE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
-         case NW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
-         case TE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
-         case BW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
-         case BE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
-         case TW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
-         case TN: return REAL_CAST(c1o54 * (drho + 3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
-         case BS: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
-         case BN: return REAL_CAST(c1o54 * (drho + 3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
-         case TS: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
-         case TNE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
-         case BSW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
-         case BNE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
-         case TSW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
-         case TSE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
-         case BNW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
-         case BSE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
-         case TNW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
+         case ZERO: return REAL_CAST(UbMath::c8o27 * (drho - cu_sq));
+         case E: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
+         case W: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq));
+         case N: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq));
+         case S: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq));
+         case T: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq));
+         case B: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq));
+         case NE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
+         case SW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
+         case SE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
+         case NW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
+         case TE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
+         case BW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
+         case BE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
+         case TW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
+         case TN: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
+         case BS: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
+         case BN: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
+         case TS: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
+         case TNE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
+         case BSW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
+         case BNE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
+         case TSW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
+         case TSE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
+         case BNW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
+         case BSE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
+         case TNW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
          default: throw UbException(UB_EXARGS, "unknown dir");
          }
       }
@@ -597,33 +592,33 @@ usage: ...
       {
          LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
 
-         feq[ZERO] = c8o27 * (drho - cu_sq);
-         feq[E] = c2o27 * (drho + 3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq);
-         feq[W] = c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
-         feq[N] = c2o27 * (drho + 3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq);
-         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[SW] = c1o54 * (drho + 3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq);
-         feq[SE] = c1o54 * (drho + 3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq);
-         feq[NW] = c1o54 * (drho + 3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq);
-         feq[TE] = c1o54 * (drho + 3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq);
-         feq[BW] = c1o54 * (drho + 3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq);
-         feq[BE] = c1o54 * (drho + 3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq);
-         feq[TW] = c1o54 * (drho + 3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq);
-         feq[TN] = c1o54 * (drho + 3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq);
-         feq[BS] = c1o54 * (drho + 3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq);
-         feq[BN] = c1o54 * (drho + 3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq);
-         feq[TS] = c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq);
-         feq[TNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq);
-         feq[BSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
-         feq[BNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq);
-         feq[TSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
-         feq[TSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq);
-         feq[BNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
-         feq[BSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
-         feq[TNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
+         feq[ZERO] = UbMath::c8o27 * (drho - cu_sq);
+         feq[E] = UbMath::c2o27 * (drho + 3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq);
+         feq[W] = UbMath::c2o27 * (drho + 3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq);
+         feq[N] = UbMath::c2o27 * (drho + 3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq);
+         feq[S] = UbMath::c2o27 * (drho + 3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq);
+         feq[T] = UbMath::c2o27 * (drho + 3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq);
+         feq[B] = UbMath::c2o27 * (drho + 3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq);
+         feq[NE] = UbMath::c1o54 * (drho + 3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq);
+         feq[SW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq);
+         feq[SE] = UbMath::c1o54 * (drho + 3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq);
+         feq[NW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq);
+         feq[TE] = UbMath::c1o54 * (drho + 3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq);
+         feq[BW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq);
+         feq[BE] = UbMath::c1o54 * (drho + 3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq);
+         feq[TW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq);
+         feq[TN] = UbMath::c1o54 * (drho + 3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq);
+         feq[BS] = UbMath::c1o54 * (drho + 3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq);
+         feq[BN] = UbMath::c1o54 * (drho + 3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq);
+         feq[TS] = UbMath::c1o54 * (drho + 3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq);
+         feq[TNE] = UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq);
+         feq[BSW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
+         feq[BNE] = UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq);
+         feq[TSW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
+         feq[TSE] = UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq);
+         feq[BNW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
+         feq[BSE] = UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
+         feq[TNW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
       }
       //////////////////////////////////////////////////////////////////////////
       static inline float getBoundaryVelocityForDirection(const int& direction, const float& bcVelocityX1, const float& bcVelocityX2, const float& bcVelocityX3)
@@ -772,7 +767,7 @@ usage: ...
       {
          LBMReal op = 1.0;
          return ((f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + 2. * (f[NE] + f[SW] + f[SE] + f[NW] + f[TE] + f[BW] + f[BE] + f[TW] + f[TN] + f[BS] + f[BN] + f[TS]) +
-            3. * (f[TNE] + f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW]) - (vx1 * vx1 + vx2 * vx2 + vx3 * vx3)) * (1 - 0.5 * op) + op * 0.5 * (rho)) * c1o3;
+            3. * (f[TNE] + f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW]) - (vx1 * vx1 + vx2 * vx2 + vx3 * vx3)) * (1 - 0.5 * op) + op * 0.5 * (rho)) * UbMath::c1o3;
 
       }
       //////////////////////////////////////////////////////////////////////////
@@ -870,7 +865,7 @@ usage: ...
          m1 = mfaac - mfaaa;
          m0 = m2 + mfaab;
          mfaaa = m0;
-         m0 += c1o36 * oMdrho;
+         m0 += UbMath::c1o36 * oMdrho;
          mfaab = m1 - m0 * vvz;
          mfaac = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -878,7 +873,7 @@ usage: ...
          m1 = mfabc - mfaba;
          m0 = m2 + mfabb;
          mfaba = m0;
-         m0 += c1o9 * oMdrho;
+         m0 += UbMath::c1o9 * oMdrho;
          mfabb = m1 - m0 * vvz;
          mfabc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -886,7 +881,7 @@ usage: ...
          m1 = mfacc - mfaca;
          m0 = m2 + mfacb;
          mfaca = m0;
-         m0 += c1o36 * oMdrho;
+         m0 += UbMath::c1o36 * oMdrho;
          mfacb = m1 - m0 * vvz;
          mfacc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -895,7 +890,7 @@ usage: ...
          m1 = mfbac - mfbaa;
          m0 = m2 + mfbab;
          mfbaa = m0;
-         m0 += c1o9 * oMdrho;
+         m0 += UbMath::c1o9 * oMdrho;
          mfbab = m1 - m0 * vvz;
          mfbac = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -903,7 +898,7 @@ usage: ...
          m1 = mfbbc - mfbba;
          m0 = m2 + mfbbb;
          mfbba = m0;
-         m0 += c4o9 * oMdrho;
+         m0 += UbMath::c4o9 * oMdrho;
          mfbbb = m1 - m0 * vvz;
          mfbbc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -911,7 +906,7 @@ usage: ...
          m1 = mfbcc - mfbca;
          m0 = m2 + mfbcb;
          mfbca = m0;
-         m0 += c1o9 * oMdrho;
+         m0 += UbMath::c1o9 * oMdrho;
          mfbcb = m1 - m0 * vvz;
          mfbcc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -920,7 +915,7 @@ usage: ...
          m1 = mfcac - mfcaa;
          m0 = m2 + mfcab;
          mfcaa = m0;
-         m0 += c1o36 * oMdrho;
+         m0 += UbMath::c1o36 * oMdrho;
          mfcab = m1 - m0 * vvz;
          mfcac = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -928,7 +923,7 @@ usage: ...
          m1 = mfcbc - mfcba;
          m0 = m2 + mfcbb;
          mfcba = m0;
-         m0 += c1o9 * oMdrho;
+         m0 += UbMath::c1o9 * oMdrho;
          mfcbb = m1 - m0 * vvz;
          mfcbc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -936,7 +931,7 @@ usage: ...
          m1 = mfccc - mfcca;
          m0 = m2 + mfccb;
          mfcca = m0;
-         m0 += c1o36 * oMdrho;
+         m0 += UbMath::c1o36 * oMdrho;
          mfccb = m1 - m0 * vvz;
          mfccc = m2 - 2. * m1 * vvz + vz2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -948,7 +943,7 @@ usage: ...
          m1 = mfaca - mfaaa;
          m0 = m2 + mfaba;
          mfaaa = m0;
-         m0 += c1o6 * oMdrho;
+         m0 += UbMath::c1o6 * oMdrho;
          mfaba = m1 - m0 * vvy;
          mfaca = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -963,7 +958,7 @@ usage: ...
          m1 = mfacc - mfaac;
          m0 = m2 + mfabc;
          mfaac = m0;
-         m0 += c1o18 * oMdrho;
+         m0 += UbMath::c1o18 * oMdrho;
          mfabc = m1 - m0 * vvy;
          mfacc = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -972,7 +967,7 @@ usage: ...
          m1 = mfbca - mfbaa;
          m0 = m2 + mfbba;
          mfbaa = m0;
-         m0 += c2o3 * oMdrho;
+         m0 += UbMath::c2o3 * oMdrho;
          mfbba = m1 - m0 * vvy;
          mfbca = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -987,7 +982,7 @@ usage: ...
          m1 = mfbcc - mfbac;
          m0 = m2 + mfbbc;
          mfbac = m0;
-         m0 += c2o9 * oMdrho;
+         m0 += UbMath::c2o9 * oMdrho;
          mfbbc = m1 - m0 * vvy;
          mfbcc = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -996,7 +991,7 @@ usage: ...
          m1 = mfcca - mfcaa;
          m0 = m2 + mfcba;
          mfcaa = m0;
-         m0 += c1o6 * oMdrho;
+         m0 += UbMath::c1o6 * oMdrho;
          mfcba = m1 - m0 * vvy;
          mfcca = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -1011,7 +1006,7 @@ usage: ...
          m1 = mfccc - mfcac;
          m0 = m2 + mfcbc;
          mfcac = m0;
-         m0 += c1o18 * oMdrho;
+         m0 += UbMath::c1o18 * oMdrho;
          mfcbc = m1 - m0 * vvy;
          mfccc = m2 - 2. * m1 * vvy + vy2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -1038,7 +1033,7 @@ usage: ...
          m1 = mfcca - mfaca;
          m0 = m2 + mfbca;
          mfaca = m0;
-         m0 += c1o3 * oMdrho;
+         m0 += UbMath::c1o3 * oMdrho;
          mfbca = m1 - m0 * vvx;
          mfcca = m2 - 2. * m1 * vvx + vx2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -1069,7 +1064,7 @@ usage: ...
          m1 = mfcac - mfaac;
          m0 = m2 + mfbac;
          mfaac = m0;
-         m0 += c1o3 * oMdrho;
+         m0 += UbMath::c1o3 * oMdrho;
          mfbac = m1 - m0 * vvx;
          mfcac = m2 - 2. * m1 * vvx + vx2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -1084,7 +1079,7 @@ usage: ...
          m1 = mfccc - mfacc;
          m0 = m2 + mfbcc;
          mfacc = m0;
-         m0 += c1o9 * oMdrho;
+         m0 += UbMath::c1o9 * oMdrho;
          mfbcc = m1 - m0 * vvx;
          mfccc = m2 - 2. * m1 * vvx + vx2 * m0;
          ////////////////////////////////////////////////////////////////////////////////////
@@ -1095,22 +1090,22 @@ usage: ...
          // LBMReal OxyyMxzz  = 1.;//2+s9;//
 
          //Cum 4.
-         //LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
-         //LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
-         //LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
+         //LBMReal CUMcbb = mfcbb - ((mfcaa + UbMath::c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+         //LBMReal CUMbcb = mfbcb - ((mfaca + UbMath::c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
+         //LBMReal CUMbbc = mfbbc - ((mfaac + UbMath::c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
 
-         // LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
-         // LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
-         // LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
+         // LBMReal CUMcbb = mfcbb - ((mfcaa + UbMath::c1o3) * mfabb + 2. * mfbba * mfbab);
+         // LBMReal CUMbcb = mfbcb - ((mfaca + UbMath::c1o3) * mfbab + 2. * mfbba * mfabb);
+         // LBMReal CUMbbc = mfbbc - ((mfaac + UbMath::c1o3) * mfbba + 2. * mfbab * mfabb);
 
-         // LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
-         // LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
-         // LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
+         // LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + UbMath::c1o3 * (mfcaa + mfaca) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
+         // LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + UbMath::c1o3 * (mfcaa + mfaac) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
+         // LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) + UbMath::c1o3 * (mfaac + mfaca) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
 
          //Cum 5.
-         // LBMReal CUMbcc = mfbcc - (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) - c1o3 * (mfbca + mfbac) * oMdrho;
-         // LBMReal CUMcbc = mfcbc - (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) - c1o3 * (mfcba + mfabc) * oMdrho;
-         // LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - c1o3 * (mfacb + mfcab) * oMdrho;
+         // LBMReal CUMbcc = mfbcc - (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) - UbMath::c1o3 * (mfbca + mfbac) * oMdrho;
+         // LBMReal CUMcbc = mfcbc - (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) - UbMath::c1o3 * (mfcba + mfabc) * oMdrho;
+         // LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - UbMath::c1o3 * (mfacb + mfcab) * oMdrho;
 
          //Cum 6.
 //         LBMReal CUMccc = mfccc + ((-4. * mfbbb * mfbbb
@@ -1120,8 +1115,8 @@ usage: ...
 //            + (4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
 //               + 2. * (mfcaa * mfaca * mfaac)
 //               + 16. * mfbba * mfbab * mfabb)
-//            - c1o3 * (mfacc + mfcac + mfcca) * oMdrho - c1o9 * oMdrho * oMdrho
-//            - c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
+//            - UbMath::c1o3 * (mfacc + mfcac + mfcca) * oMdrho - UbMath::c1o9 * oMdrho * oMdrho
+//            - UbMath::c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
 //            + (2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
 //               + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) + c1o27 * oMdrho;
 
@@ -1130,17 +1125,17 @@ usage: ...
          LBMReal mxxMyy = mfcaa - mfaca;
          LBMReal mxxMzz = mfcaa - mfaac;
 
-         LBMReal dxux = -c1o2 * collFactorF * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
-         LBMReal dyuy = dxux + collFactorF * c3o2 * mxxMyy;
-         LBMReal dzuz = dxux + collFactorF * c3o2 * mxxMzz;
+         LBMReal dxux = -UbMath::c1o2 * collFactorF * (mxxMyy + mxxMzz) + UbMath::c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+         LBMReal dyuy = dxux + collFactorF * UbMath::c3o2 * mxxMyy;
+         LBMReal dzuz = dxux + collFactorF * UbMath::c3o2 * mxxMzz;
 
-         LBMReal Dxy = -three * collFactorF * mfbba;
-         LBMReal Dxz = -three * collFactorF * mfbab;
-         LBMReal Dyz = -three * collFactorF * mfabb;
+         LBMReal Dxy = -UbMath::three * collFactorF * mfbba;
+         LBMReal Dxz = -UbMath::three * collFactorF * mfbab;
+         LBMReal Dyz = -UbMath::three * collFactorF * mfabb;
 
          
          //TODO: may be factor 2
-         return sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+         return sqrt(UbMath::c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + UbMath::one);
       }
    }
 #endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
index 42c3240077e4b6cc3e3a279608bca0325b304b8c..3a052547388a0e8c23acf15340546bd9aaaadf1c 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
@@ -10,6 +10,8 @@
 
 #define PROOF_CORRECTNESS
 
+using namespace UbMath;
+
 //////////////////////////////////////////////////////////////////////////
 RheologyK17LBMKernel::RheologyK17LBMKernel()
 {
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
index adbca6c6c57eab9fd520fe7e7fe1f8d7c055053c..095f9678a3256957995b43f6cee81805a56bce0f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
@@ -8,6 +8,8 @@
 
 #define PROOF_CORRECTNESS
 
+using namespace UbMath;
+
 //////////////////////////////////////////////////////////////////////////
 ThixotropyExpLBMKernel::ThixotropyExpLBMKernel()
 {
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
index 949e7777d644c5b9051611fbcc4b8277b37d630d..eb313cdd74a6417f3fab6f669b25b08f712667dc 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
@@ -8,6 +8,8 @@
 
 #define PROOF_CORRECTNESS
 
+using namespace UbMath;
+
 //////////////////////////////////////////////////////////////////////////
 ThixotropyLBMKernel::ThixotropyLBMKernel()
 {
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
index d0d0542b63ba68136c0c23ec9cbc86bb37354ea5..1685eb6e1393055f0e04f7f4289fd295e805c33f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
@@ -8,6 +8,8 @@
 
 #define PROOF_CORRECTNESS
 
+using namespace UbMath;
+
 ThixotropyModelLBMKernel::ThixotropyModelLBMKernel() : forcingX1(0), forcingX2(0), forcingX3(0)
 {
    compressible = false;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
index 8839818a70892d30b57de2938360b0199af3e234..5bb4770fd7ddeb0ea599315d8160fcb98ae02b0b 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -7,7 +7,10 @@
 #include "LBMKernel.h"
 #include "Thixotropy.h"
 
-#define PROOF_CORRECTNESS
+#define PROOF_CORRECTNES
+
+using namespace UbMath;
+
 
 ThixotropyModelLBMKernel2::ThixotropyModelLBMKernel2() : forcingX1(0), forcingX2(0), forcingX3(0)
 {