From 29d4d6befac9d77f6be7be288352fb4f31137b31 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Wed, 3 Jan 2018 11:16:16 +0100
Subject: [PATCH] SpongeLayerBlockVisitor is extented

---
 .../Applications/DLR-F16-Solid/f16-solid.cfg  |   21 +-
 source/Applications/DLR-F16-Solid/f16.cpp     |   68 +-
 ...bleOffsetMomentsInterpolationProcessor.cpp | 1140 ++++++++---------
 .../Visitors/SpongeLayerBlockVisitor.cpp      |    8 +-
 .../Visitors/SpongeLayerBlockVisitor.h        |    3 +-
 5 files changed, 639 insertions(+), 601 deletions(-)

diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg
index 5f8088a7a..26854c809 100644
--- a/source/Applications/DLR-F16-Solid/f16-solid.cfg
+++ b/source/Applications/DLR-F16-Solid/f16-solid.cfg
@@ -1,13 +1,14 @@
 pathOut = d:/temp/DLR-F16-Solid
 pathGeo = d:/Projects/SFB880/DLR-F16/A1_Forschungsdaten_Profilgeometrie_STL_CATIA_Rossian
-fngFileWhole = f16-ascii.stl
-ailing-edge-ascii.stl
+#fngFileWhole = f16-ascii.stl
+fngFileWhole = grundgeometrie-direkter-export.stl
+
 zigZagTape = 2zackenbaender0.stl
 
 pathReInit = /work/koskuche/DLR-F16_L7
 stepReInit = 10000
 
-numOfThreads = 1
+numOfThreads = 4
 availMem = 15e9
 
 logToFile = false
@@ -17,7 +18,7 @@ logToFile = false
 #deltaXfine = 0.00075
 #boundingBox = -0.90 1.20 0.035 0.065 -0.65 0.65
 #boundingBox = -0.90 2.1 0.035 0.065 -0.66 0.66
-boundingBox = -0.2 1.1 0.035 0.065 -0.33 0.33
+boundingBox = -0.12 1.2 0.035 0.065 -0.24 0.24
 blockNx = 10 10 10
 
 #deltaXfine = 13393e-9
@@ -30,7 +31,7 @@ blockNx = 10 10 10
 
 refineLevel = 0
 
-deltaXfine = 0.012 #level 0
+deltaXfine = 0.006 #level 0
 #deltaXfine = 0.0015 #level 0
 #deltaXfine = 0.00075 #level 1
 #deltaXfine = 0.000375 #level 2
@@ -47,8 +48,8 @@ refineDistance = 0.6
 
 writeBlocks = true
 
-newStart = true
-restartStep = 1000
+newStart = false
+restartStep = 10000
 
 cpStep = 1000
 cpStart = 1000
@@ -56,13 +57,13 @@ cpStart = 1000
 outTimeStep = 1000
 outTimeStart = 1000
 
-endTime = 100000
+endTime = 20000
 
 #Cp
 pcpStart = 1000000
 pcpStop  = 1000000
 
-timeAvStart = 1000
-timeAvStop  = 67000
+timeAvStart = 10000
+timeAvStop  = 20000
 
 nupsStep = 100 100 10000000
diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp
index b99153d50..7924aa80f 100644
--- a/source/Applications/DLR-F16-Solid/f16.cpp
+++ b/source/Applications/DLR-F16-Solid/f16.cpp
@@ -139,9 +139,11 @@ void run(string configname)
       fct.DefineConst("U", uLB);
       BCAdapterPtr velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
       velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
+      //velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithm()));
 
       BCAdapterPtr outflowBCAdapter(new DensityBCAdapter(rhoLB));
       outflowBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NonReflectingOutflowBCAlgorithm()));
+      //outflowBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NonEqDensityBCAlgorithm()));
 
       BoundaryConditionsBlockVisitor bcVisitor;
       bcVisitor.addBC(noSlipBCAdapter);
@@ -208,15 +210,15 @@ void run(string configname)
 
         
 
-         GbObject3DPtr fngMeshWhole(new GbCylinder3D(15.0, 0.0, 0.0, 15.0, 100.0, 0.0, 25.0));
-         GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut + "/geo/fngMeshWholeCylinder", WbWriterVtkXmlBinary::getInstance());
+         //GbObject3DPtr fngMeshWhole(new GbCylinder3D(15.0, 0.0, 0.0, 15.0, 300.0, 0.0, 25.0));
+         //GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut + "/geo/fngMeshWholeCylinder", WbWriterVtkXmlBinary::getInstance());
 
-         //GbTriFaceMesh3DPtr fngMeshWhole;
-         //if (myid==0) UBLOG(logINFO, "Read fngFileWhole:start");
-         //fngMeshWhole = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(pathGeo+"/"+fngFileWhole, "fngMeshWhole", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
-         //if (myid==0) UBLOG(logINFO, "Read fngFileWhole:end");
-         //fngMeshWhole->rotate(0.0, 0.5, 0.0);
-         //if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut+"/geo/fngMeshWhole2", WbWriterVtkXmlBinary::getInstance());
+         GbTriFaceMesh3DPtr fngMeshWhole;
+         if (myid==0) UBLOG(logINFO, "Read fngFileWhole:start");
+         fngMeshWhole = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(pathGeo+"/"+fngFileWhole, "fngMeshWhole", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+         if (myid==0) UBLOG(logINFO, "Read fngFileWhole:end");
+         fngMeshWhole->rotate(0.0, 0.5, 0.0);
+         if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut+"/geo/fngMeshWhole2", WbWriterVtkXmlBinary::getInstance());
 
          ////////////////////////////////////////////////////////////////////////////
          //// Zackenband
@@ -260,8 +262,8 @@ void run(string configname)
          }
          //////////////////////////////////////////////////////////////////////////
          Interactor3DPtr fngIntrWhole;
-         fngIntrWhole = D3Q27InteractorPtr(new D3Q27Interactor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
-         //fngIntrWhole = D3Q27TriFaceMeshInteractorPtr(new D3Q27TriFaceMeshInteractor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
+         //fngIntrWhole = D3Q27InteractorPtr(new D3Q27Interactor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
+         fngIntrWhole = D3Q27TriFaceMeshInteractorPtr(new D3Q27TriFaceMeshInteractor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
 
          //D3Q27TriFaceMeshInteractorPtr triBand1Interactor(new D3Q27TriFaceMeshInteractor(meshBand1, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
          //D3Q27TriFaceMeshInteractorPtr triBand2Interactor(new D3Q27TriFaceMeshInteractor(meshBand2, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
@@ -422,7 +424,7 @@ void run(string configname)
          D3Q27InteractorPtr inflowIntr = D3Q27InteractorPtr(new D3Q27Interactor(geoInflow, grid, velBCAdapter, Interactor3D::SOLID));
 
          //outflow
-         D3Q27InteractorPtr outflowIntr = D3Q27InteractorPtr(new D3Q27Interactor(geoOutflow, grid, outflowBCAdapter, Interactor3D::SOLID));
+         D3Q27InteractorPtr outflowIntr = D3Q27InteractorPtr(new D3Q27Interactor(geoOutflow, grid, velBCAdapter, Interactor3D::SOLID));
 
          ////////////////////////////////////////////
          //METIS
@@ -552,24 +554,39 @@ void run(string configname)
          ////////////////////////////////////////////////////////////////////////////
          //GbCuboid3DPtr spongeLayerX1min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1+75, g_maxX2+blockLength, g_maxX3+blockLength));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1min.get(), pathOut+"/geo/spongeLayerX1min", WbWriterVtkXmlASCII::getInstance());
-         //SpongeLayerBlockVisitor slVisitorX1min(spongeLayerX1min);
+         //SpongeLayerBlockVisitor slVisitorX1min(spongeLayerX1min, 1.0);
          //grid->accept(slVisitorX1min);
 
-         //GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-75, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+         //GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-3.0*blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance());
-         //SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max);
+         //SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, 1.0);
          //grid->accept(slVisitorX1max);
 
          //GbCuboid3DPtr spongeLayerX3min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+75));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3min.get(), pathOut+"/geo/spongeLayerX3min", WbWriterVtkXmlASCII::getInstance());
-         //SpongeLayerBlockVisitor slVisitorX3min(spongeLayerX3min);
+         //SpongeLayerBlockVisitor slVisitorX3min(spongeLayerX3min, 1.0);
          //grid->accept(slVisitorX3min);
 
          //GbCuboid3DPtr spongeLayerX3max(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_maxX3-75, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3max.get(), pathOut+"/geo/spongeLayerX3max", WbWriterVtkXmlASCII::getInstance());
-         //SpongeLayerBlockVisitor slVisitorX3max(spongeLayerX3max);
+         //SpongeLayerBlockVisitor slVisitorX3max(spongeLayerX3max, 1.0);
          //grid->accept(slVisitorX3max);
 
+         //GbCuboid3DPtr spongeLayerX1max1(new GbCuboid3D(g_maxX1-75*4.0, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-75*2.0, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max1.get(), pathOut+"/geo/spongeLayerX1max1", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max1(spongeLayerX1max1, 1.9);
+         //grid->accept(slVisitorX1max1);
+
+         //GbCuboid3DPtr spongeLayerX1max2(new GbCuboid3D(g_maxX1-75*2.0, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-75.0, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max2.get(), pathOut+"/geo/spongeLayerX1max2", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max2(spongeLayerX1max2, 1.5);
+         //grid->accept(slVisitorX1max2);
+
+         //GbCuboid3DPtr spongeLayerX1max3(new GbCuboid3D(g_maxX1-75, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max3.get(), pathOut+"/geo/spongeLayerX1max3", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max3(spongeLayerX1max3, 1.0);
+         //grid->accept(slVisitorX1max3);
+
          /////////////////////////////////////////////////////////////////////////////
          if (myid==0) UBLOG(logINFO, "Preprozess - end");
       }
@@ -587,6 +604,21 @@ void run(string configname)
          PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
          grid->accept(pqPartVisitor);
 
+         //GbCuboid3DPtr spongeLayerX1max1(new GbCuboid3D(g_maxX1-75*4.0, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-75*2.0, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max1.get(), pathOut+"/geo/spongeLayerX1max1", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max1(spongeLayerX1max1, 1.9);
+         //grid->accept(slVisitorX1max1);
+
+         //GbCuboid3DPtr spongeLayerX1max2(new GbCuboid3D(g_maxX1-75*2.0, g_minX2-blockLength, g_minX3-blockLength, g_maxX1-75.0, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max2.get(), pathOut+"/geo/spongeLayerX1max2", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max2(spongeLayerX1max2, 1.5);
+         //grid->accept(slVisitorX1max2);
+
+         //GbCuboid3DPtr spongeLayerX1max3(new GbCuboid3D(g_maxX1-75, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max3.get(), pathOut+"/geo/spongeLayerX1max3", WbWriterVtkXmlASCII::getInstance());
+         //SpongeLayerBlockVisitor slVisitorX1max3(spongeLayerX1max3, 1.0);
+         //grid->accept(slVisitorX1max3);
+
       }
 
       UbSchedulerPtr nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2]));
@@ -613,8 +645,8 @@ void run(string configname)
       //TimeseriesCoProcessor tsp1(grid, stepMV, mic1, pathOut+"/mic/mic1", comm);
 
       //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, stepSch));
-      CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, stepSch, CalculationManager::MPI));
-      //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, tavSch, CalculationManager::PrePostBc));
+      //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, stepSch, CalculationManager::MPI));
+      CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, tavSch, CalculationManager::PrePostBc));
 
       if (myid==0) UBLOG(logINFO, "Simulation-start");
       calculation->calculate();
diff --git a/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp b/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
index dc81684e5..83590840f 100644
--- a/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
+++ b/source/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
@@ -489,32 +489,32 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeCF(LBM
    LBMReal vx2  = b0 + 0.25*( xs*bx + ys*by + zs*bz) + 0.0625*(bxx + xs*ys*bxy + xs*zs*bxz + byy + ys*zs*byz + bzz) + 0.015625*(xs*ys*zs*bxyz);
    LBMReal vx3  = c0 + 0.25*( xs*cx + ys*cy + zs*cz) + 0.0625*(cxx + xs*ys*cxy + xs*zs*cxz + cyy + ys*zs*cyz + czz) + 0.015625*(xs*ys*zs*cxyz);
 
-   LBMReal mfcbb = zeroReal;
-   LBMReal mfabb = zeroReal;
-   LBMReal mfbcb = zeroReal;
-   LBMReal mfbab = zeroReal;
-   LBMReal mfbbc = zeroReal;
-   LBMReal mfbba = zeroReal;
-   LBMReal mfccb = zeroReal;
-   LBMReal mfaab = zeroReal;
-   LBMReal mfcab = zeroReal;
-   LBMReal mfacb = zeroReal;
-   LBMReal mfcbc = zeroReal;
-   LBMReal mfaba = zeroReal;
-   LBMReal mfcba = zeroReal;
-   LBMReal mfabc = zeroReal;
-   LBMReal mfbcc = zeroReal;
-   LBMReal mfbaa = zeroReal;
-   LBMReal mfbca = zeroReal;
-   LBMReal mfbac = zeroReal;
-   LBMReal mfbbb = zeroReal;
-   LBMReal mfccc = zeroReal;
-   LBMReal mfaac = zeroReal;
-   LBMReal mfcac = zeroReal;
-   LBMReal mfacc = zeroReal;
-   LBMReal mfcca = zeroReal;
-   LBMReal mfaaa = zeroReal;
-   LBMReal mfcaa = zeroReal;
+   LBMReal mfcbb = zeroReal;
+   LBMReal mfabb = zeroReal;
+   LBMReal mfbcb = zeroReal;
+   LBMReal mfbab = zeroReal;
+   LBMReal mfbbc = zeroReal;
+   LBMReal mfbba = zeroReal;
+   LBMReal mfccb = zeroReal;
+   LBMReal mfaab = zeroReal;
+   LBMReal mfcab = zeroReal;
+   LBMReal mfacb = zeroReal;
+   LBMReal mfcbc = zeroReal;
+   LBMReal mfaba = zeroReal;
+   LBMReal mfcba = zeroReal;
+   LBMReal mfabc = zeroReal;
+   LBMReal mfbcc = zeroReal;
+   LBMReal mfbaa = zeroReal;
+   LBMReal mfbca = zeroReal;
+   LBMReal mfbac = zeroReal;
+   LBMReal mfbbb = zeroReal;
+   LBMReal mfccc = zeroReal;
+   LBMReal mfaac = zeroReal;
+   LBMReal mfcac = zeroReal;
+   LBMReal mfacc = zeroReal;
+   LBMReal mfcca = zeroReal;
+   LBMReal mfaaa = zeroReal;
+   LBMReal mfcaa = zeroReal;
    LBMReal mfaca = zeroReal;
 
    mfaaa = press; // if drho is interpolated directly
@@ -572,237 +572,237 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeCF(LBM
    ////////////////////////////////////////////////////////////////////////////////////
    //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
    ////////////////////////////////////////////////////////////////////////////////////
-   // Z - Dir
-   LBMReal m0 =  mfaac * c1o2 +      mfaab * (vx3 - c1o2) + (mfaaa + one * oMdrho) * (vx3Sq - vx3) * c1o2;
-   LBMReal m1 = -mfaac        - two * mfaab *  vx3         +  mfaaa                * (one - vx3Sq)              - one * oMdrho * vx3Sq;
-   LBMReal m2 =  mfaac * c1o2 +      mfaab * (vx3 + c1o2) + (mfaaa + one * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfaaa = m0;
-   mfaab = m1;
-   mfaac = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfabc * c1o2 +      mfabb * (vx3 - c1o2) + mfaba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfabc        - two * mfabb *  vx3         + mfaba * (one - vx3Sq);
-   m2 =  mfabc * c1o2 +      mfabb * (vx3 + c1o2) + mfaba * (vx3Sq + vx3) * c1o2;
-   mfaba = m0;
-   mfabb = m1;
-   mfabc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacc * c1o2 +      mfacb * (vx3 - c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfacc        - two * mfacb *  vx3         +  mfaca                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
-   m2 =  mfacc * c1o2 +      mfacb * (vx3 + c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfaca = m0;
-   mfacb = m1;
-   mfacc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfbac * c1o2 +      mfbab * (vx3 - c1o2) + mfbaa * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbac        - two * mfbab *  vx3         + mfbaa * (one - vx3Sq);
-   m2 =  mfbac * c1o2 +      mfbab * (vx3 + c1o2) + mfbaa * (vx3Sq + vx3) * c1o2;
-   mfbaa = m0;
-   mfbab = m1;
-   mfbac = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbbc * c1o2 +      mfbbb * (vx3 - c1o2) + mfbba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbbc        - two * mfbbb *  vx3         + mfbba * (one - vx3Sq);
-   m2 =  mfbbc * c1o2 +      mfbbb * (vx3 + c1o2) + mfbba * (vx3Sq + vx3) * c1o2;
-   mfbba = m0;
-   mfbbb = m1;
-   mfbbc = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcc * c1o2 +      mfbcb * (vx3 - c1o2) + mfbca * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbcc        - two * mfbcb *  vx3         + mfbca * (one - vx3Sq);
-   m2 =  mfbcc * c1o2 +      mfbcb * (vx3 + c1o2) + mfbca * (vx3Sq + vx3) * c1o2;
-   mfbca = m0;
-   mfbcb = m1;
-   mfbcc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcac * c1o2 +      mfcab * (vx3 - c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfcac        - two * mfcab *  vx3         +  mfcaa                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
-   m2 =  mfcac * c1o2 +      mfcab * (vx3 + c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfcaa = m0;
-   mfcab = m1;
-   mfcac = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbc * c1o2 +      mfcbb * (vx3 - c1o2) + mfcba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfcbc        - two * mfcbb *  vx3         + mfcba * (one - vx3Sq);
-   m2 =  mfcbc * c1o2 +      mfcbb * (vx3 + c1o2) + mfcba * (vx3Sq + vx3) * c1o2;
-   mfcba = m0;
-   mfcbb = m1;
-   mfcbc = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfccb * (vx3 - c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfccc        - two * mfccb *  vx3         +  mfcca                  * (one - vx3Sq)              - c1o9 * oMdrho * vx3Sq;
-   m2 =  mfccc * c1o2 +      mfccb * (vx3 + c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfcca = m0;
-   mfccb = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18   Konditionieren
-   ////////////////////////////////////////////////////////////////////////////////////
-   // Y - Dir
-   m0 =  mfaca * c1o2 +      mfaba * (vx2 - c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfaca        - two * mfaba *  vx2         +  mfaaa                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
-   m2 =  mfaca * c1o2 +      mfaba * (vx2 + c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaaa = m0;
-   mfaba = m1;
-   mfaca = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacb * c1o2 +      mfabb * (vx2 - c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfacb        - two * mfabb *  vx2         +  mfaab                  * (one - vx2Sq)              - c2o3 * oMdrho * vx2Sq;
-   m2 =  mfacb * c1o2 +      mfabb * (vx2 + c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaab = m0;
-   mfabb = m1;
-   mfacb = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacc * c1o2 +      mfabc * (vx2 - c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfacc        - two * mfabc *  vx2         +  mfaac                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
-   m2 =  mfacc * c1o2 +      mfabc * (vx2 + c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaac = m0;
-   mfabc = m1;
-   mfacc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfbca * c1o2 +      mfbba * (vx2 - c1o2) + mfbaa * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbca        - two * mfbba *  vx2         + mfbaa * (one - vx2Sq);
-   m2 =  mfbca * c1o2 +      mfbba * (vx2 + c1o2) + mfbaa * (vx2Sq + vx2) * c1o2;
-   mfbaa = m0;
-   mfbba = m1;
-   mfbca = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcb * c1o2 +      mfbbb * (vx2 - c1o2) + mfbab * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbcb        - two * mfbbb *  vx2         + mfbab * (one - vx2Sq);
-   m2 =  mfbcb * c1o2 +      mfbbb * (vx2 + c1o2) + mfbab * (vx2Sq + vx2) * c1o2;
-   mfbab = m0;
-   mfbbb = m1;
-   mfbcb = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcc * c1o2 +      mfbbc * (vx2 - c1o2) + mfbac * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbcc        - two * mfbbc *  vx2         + mfbac * (one - vx2Sq);
-   m2 =  mfbcc * c1o2 +      mfbbc * (vx2 + c1o2) + mfbac * (vx2Sq + vx2) * c1o2;
-   mfbac = m0;
-   mfbbc = m1;
-   mfbcc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcca * c1o2 +      mfcba * (vx2 - c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfcca        - two * mfcba *  vx2         +  mfcaa                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
-   m2 =  mfcca * c1o2 +      mfcba * (vx2 + c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcaa = m0;
-   mfcba = m1;
-   mfcca = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccb * c1o2 +      mfcbb * (vx2 - c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfccb        - two * mfcbb *  vx2         +  mfcab                  * (one - vx2Sq)              - c2o9 * oMdrho * vx2Sq;
-   m2 =  mfccb * c1o2 +      mfcbb * (vx2 + c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcab = m0;
-   mfcbb = m1;
-   mfccb = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfcbc * (vx2 - c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfccc        - two * mfcbc *  vx2         +  mfcac                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
-   m2 =  mfccc * c1o2 +      mfcbc * (vx2 + c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcac = m0;
-   mfcbc = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
-   ////////////////////////////////////////////////////////////////////////////////////
-   // X - Dir
-   m0 =  mfcaa * c1o2 +      mfbaa * (vx1 - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcaa        - two * mfbaa *  vx1         +  mfaaa                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcaa * c1o2 +      mfbaa * (vx1 + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaaa = m0;
-   mfbaa = m1;
-   mfcaa = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcba * c1o2 +      mfbba * (vx1 - c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcba        - two * mfbba *  vx1         +  mfaba                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcba * c1o2 +      mfbba * (vx1 + c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaba = m0;
-   mfbba = m1;
-   mfcba = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcca * c1o2 +      mfbca * (vx1 - c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcca        - two * mfbca *  vx1         +  mfaca                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcca * c1o2 +      mfbca * (vx1 + c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaca = m0;
-   mfbca = m1;
-   mfcca = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcab * c1o2 +      mfbab * (vx1 - c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcab        - two * mfbab *  vx1         +  mfaab                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcab * c1o2 +      mfbab * (vx1 + c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaab = m0;
-   mfbab = m1;
-   mfcab = m2;
-   ///////////b////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbb * c1o2 +      mfbbb * (vx1 - c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcbb        - two * mfbbb *  vx1         +  mfabb                  * (one - vx1Sq)              - c4o9 * oMdrho * vx1Sq;
-   m2 =  mfcbb * c1o2 +      mfbbb * (vx1 + c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfabb = m0;
-   mfbbb = m1;
-   mfcbb = m2;
-   ///////////b////////////////////////////////////////////////////////////////////////
-   m0 =  mfccb * c1o2 +      mfbcb * (vx1 - c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfccb        - two * mfbcb *  vx1         +  mfacb                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfccb * c1o2 +      mfbcb * (vx1 + c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfacb = m0;
-   mfbcb = m1;
-   mfccb = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcac * c1o2 +      mfbac * (vx1 - c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcac        - two * mfbac *  vx1         +  mfaac                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcac * c1o2 +      mfbac * (vx1 + c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaac = m0;
-   mfbac = m1;
-   mfcac = m2;
-   ///////////c////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbc * c1o2 +      mfbbc * (vx1 - c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcbc        - two * mfbbc *  vx1         +  mfabc                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcbc * c1o2 +      mfbbc * (vx1 + c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfabc = m0;
-   mfbbc = m1;
-   mfcbc = m2;
-   ///////////c////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfbcc * (vx1 - c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfccc        - two * mfbcc *  vx1         +  mfacc                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfccc * c1o2 +      mfbcc * (vx1 + c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfacc = m0;
-   mfbcc = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-
-   f[E]    = mfcbb;
-   f[W]    = mfabb;
-   f[N]    = mfbcb;
-   f[S]    = mfbab;
-   f[T]    = mfbbc;
-   f[B]    = mfbba;
-   f[NE]   = mfccb;
-   f[SW]   = mfaab;
-   f[SE]   = mfcab;
-   f[NW]   = mfacb;
-   f[TE]   = mfcbc;
-   f[BW]   = mfaba;
-   f[BE]   = mfcba;
-   f[TW]   = mfabc;
-   f[TN]   = mfbcc;
-   f[BS]   = mfbaa;
-   f[BN]   = mfbca;
-   f[TS]   = mfbac;
-   f[ZERO] = mfbbb;
-   f[TNE]  = mfccc;
-   f[TSE]  = mfcac;
-   f[BNE]  = mfcca;
-   f[BSE]  = mfcaa;
-   f[TNW]  = mfacc;
-   f[TSW]  = mfaac;
-   f[BNW]  = mfaca;
+   // Z - Dir
+   LBMReal m0 =  mfaac * c1o2 +      mfaab * (vx3 - c1o2) + (mfaaa + one * oMdrho) * (vx3Sq - vx3) * c1o2;
+   LBMReal m1 = -mfaac        - two * mfaab *  vx3         +  mfaaa                * (one - vx3Sq)              - one * oMdrho * vx3Sq;
+   LBMReal m2 =  mfaac * c1o2 +      mfaab * (vx3 + c1o2) + (mfaaa + one * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfaaa = m0;
+   mfaab = m1;
+   mfaac = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfabc * c1o2 +      mfabb * (vx3 - c1o2) + mfaba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfabc        - two * mfabb *  vx3         + mfaba * (one - vx3Sq);
+   m2 =  mfabc * c1o2 +      mfabb * (vx3 + c1o2) + mfaba * (vx3Sq + vx3) * c1o2;
+   mfaba = m0;
+   mfabb = m1;
+   mfabc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacc * c1o2 +      mfacb * (vx3 - c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfacc        - two * mfacb *  vx3         +  mfaca                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
+   m2 =  mfacc * c1o2 +      mfacb * (vx3 + c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfaca = m0;
+   mfacb = m1;
+   mfacc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfbac * c1o2 +      mfbab * (vx3 - c1o2) + mfbaa * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbac        - two * mfbab *  vx3         + mfbaa * (one - vx3Sq);
+   m2 =  mfbac * c1o2 +      mfbab * (vx3 + c1o2) + mfbaa * (vx3Sq + vx3) * c1o2;
+   mfbaa = m0;
+   mfbab = m1;
+   mfbac = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbbc * c1o2 +      mfbbb * (vx3 - c1o2) + mfbba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbbc        - two * mfbbb *  vx3         + mfbba * (one - vx3Sq);
+   m2 =  mfbbc * c1o2 +      mfbbb * (vx3 + c1o2) + mfbba * (vx3Sq + vx3) * c1o2;
+   mfbba = m0;
+   mfbbb = m1;
+   mfbbc = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcc * c1o2 +      mfbcb * (vx3 - c1o2) + mfbca * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbcc        - two * mfbcb *  vx3         + mfbca * (one - vx3Sq);
+   m2 =  mfbcc * c1o2 +      mfbcb * (vx3 + c1o2) + mfbca * (vx3Sq + vx3) * c1o2;
+   mfbca = m0;
+   mfbcb = m1;
+   mfbcc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcac * c1o2 +      mfcab * (vx3 - c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfcac        - two * mfcab *  vx3         +  mfcaa                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
+   m2 =  mfcac * c1o2 +      mfcab * (vx3 + c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfcaa = m0;
+   mfcab = m1;
+   mfcac = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbc * c1o2 +      mfcbb * (vx3 - c1o2) + mfcba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfcbc        - two * mfcbb *  vx3         + mfcba * (one - vx3Sq);
+   m2 =  mfcbc * c1o2 +      mfcbb * (vx3 + c1o2) + mfcba * (vx3Sq + vx3) * c1o2;
+   mfcba = m0;
+   mfcbb = m1;
+   mfcbc = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfccb * (vx3 - c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfccc        - two * mfccb *  vx3         +  mfcca                  * (one - vx3Sq)              - c1o9 * oMdrho * vx3Sq;
+   m2 =  mfccc * c1o2 +      mfccb * (vx3 + c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfcca = m0;
+   mfccb = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18   Konditionieren
+   ////////////////////////////////////////////////////////////////////////////////////
+   // Y - Dir
+   m0 =  mfaca * c1o2 +      mfaba * (vx2 - c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfaca        - two * mfaba *  vx2         +  mfaaa                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
+   m2 =  mfaca * c1o2 +      mfaba * (vx2 + c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaaa = m0;
+   mfaba = m1;
+   mfaca = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacb * c1o2 +      mfabb * (vx2 - c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfacb        - two * mfabb *  vx2         +  mfaab                  * (one - vx2Sq)              - c2o3 * oMdrho * vx2Sq;
+   m2 =  mfacb * c1o2 +      mfabb * (vx2 + c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaab = m0;
+   mfabb = m1;
+   mfacb = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacc * c1o2 +      mfabc * (vx2 - c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfacc        - two * mfabc *  vx2         +  mfaac                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
+   m2 =  mfacc * c1o2 +      mfabc * (vx2 + c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaac = m0;
+   mfabc = m1;
+   mfacc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfbca * c1o2 +      mfbba * (vx2 - c1o2) + mfbaa * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbca        - two * mfbba *  vx2         + mfbaa * (one - vx2Sq);
+   m2 =  mfbca * c1o2 +      mfbba * (vx2 + c1o2) + mfbaa * (vx2Sq + vx2) * c1o2;
+   mfbaa = m0;
+   mfbba = m1;
+   mfbca = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcb * c1o2 +      mfbbb * (vx2 - c1o2) + mfbab * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbcb        - two * mfbbb *  vx2         + mfbab * (one - vx2Sq);
+   m2 =  mfbcb * c1o2 +      mfbbb * (vx2 + c1o2) + mfbab * (vx2Sq + vx2) * c1o2;
+   mfbab = m0;
+   mfbbb = m1;
+   mfbcb = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcc * c1o2 +      mfbbc * (vx2 - c1o2) + mfbac * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbcc        - two * mfbbc *  vx2         + mfbac * (one - vx2Sq);
+   m2 =  mfbcc * c1o2 +      mfbbc * (vx2 + c1o2) + mfbac * (vx2Sq + vx2) * c1o2;
+   mfbac = m0;
+   mfbbc = m1;
+   mfbcc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcca * c1o2 +      mfcba * (vx2 - c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfcca        - two * mfcba *  vx2         +  mfcaa                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
+   m2 =  mfcca * c1o2 +      mfcba * (vx2 + c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcaa = m0;
+   mfcba = m1;
+   mfcca = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccb * c1o2 +      mfcbb * (vx2 - c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfccb        - two * mfcbb *  vx2         +  mfcab                  * (one - vx2Sq)              - c2o9 * oMdrho * vx2Sq;
+   m2 =  mfccb * c1o2 +      mfcbb * (vx2 + c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcab = m0;
+   mfcbb = m1;
+   mfccb = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfcbc * (vx2 - c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfccc        - two * mfcbc *  vx2         +  mfcac                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
+   m2 =  mfccc * c1o2 +      mfcbc * (vx2 + c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcac = m0;
+   mfcbc = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+   ////////////////////////////////////////////////////////////////////////////////////
+   // X - Dir
+   m0 =  mfcaa * c1o2 +      mfbaa * (vx1 - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcaa        - two * mfbaa *  vx1         +  mfaaa                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcaa * c1o2 +      mfbaa * (vx1 + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaaa = m0;
+   mfbaa = m1;
+   mfcaa = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcba * c1o2 +      mfbba * (vx1 - c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcba        - two * mfbba *  vx1         +  mfaba                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcba * c1o2 +      mfbba * (vx1 + c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaba = m0;
+   mfbba = m1;
+   mfcba = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcca * c1o2 +      mfbca * (vx1 - c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcca        - two * mfbca *  vx1         +  mfaca                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcca * c1o2 +      mfbca * (vx1 + c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaca = m0;
+   mfbca = m1;
+   mfcca = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcab * c1o2 +      mfbab * (vx1 - c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcab        - two * mfbab *  vx1         +  mfaab                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcab * c1o2 +      mfbab * (vx1 + c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaab = m0;
+   mfbab = m1;
+   mfcab = m2;
+   ///////////b////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbb * c1o2 +      mfbbb * (vx1 - c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcbb        - two * mfbbb *  vx1         +  mfabb                  * (one - vx1Sq)              - c4o9 * oMdrho * vx1Sq;
+   m2 =  mfcbb * c1o2 +      mfbbb * (vx1 + c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfabb = m0;
+   mfbbb = m1;
+   mfcbb = m2;
+   ///////////b////////////////////////////////////////////////////////////////////////
+   m0 =  mfccb * c1o2 +      mfbcb * (vx1 - c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfccb        - two * mfbcb *  vx1         +  mfacb                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfccb * c1o2 +      mfbcb * (vx1 + c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfacb = m0;
+   mfbcb = m1;
+   mfccb = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcac * c1o2 +      mfbac * (vx1 - c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcac        - two * mfbac *  vx1         +  mfaac                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcac * c1o2 +      mfbac * (vx1 + c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaac = m0;
+   mfbac = m1;
+   mfcac = m2;
+   ///////////c////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbc * c1o2 +      mfbbc * (vx1 - c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcbc        - two * mfbbc *  vx1         +  mfabc                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcbc * c1o2 +      mfbbc * (vx1 + c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfabc = m0;
+   mfbbc = m1;
+   mfcbc = m2;
+   ///////////c////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfbcc * (vx1 - c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfccc        - two * mfbcc *  vx1         +  mfacc                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfccc * c1o2 +      mfbcc * (vx1 + c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfacc = m0;
+   mfbcc = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+
+   f[E]    = mfcbb;
+   f[W]    = mfabb;
+   f[N]    = mfbcb;
+   f[S]    = mfbab;
+   f[T]    = mfbbc;
+   f[B]    = mfbba;
+   f[NE]   = mfccb;
+   f[SW]   = mfaab;
+   f[SE]   = mfcab;
+   f[NW]   = mfacb;
+   f[TE]   = mfcbc;
+   f[BW]   = mfaba;
+   f[BE]   = mfcba;
+   f[TW]   = mfabc;
+   f[TN]   = mfbcc;
+   f[BS]   = mfbaa;
+   f[BN]   = mfbca;
+   f[TS]   = mfbac;
+   f[ZERO] = mfbbb;
+   f[TNE]  = mfccc;
+   f[TSE]  = mfcac;
+   f[BNE]  = mfcca;
+   f[BSE]  = mfcaa;
+   f[TNW]  = mfacc;
+   f[TSW]  = mfaac;
+   f[BNW]  = mfaca;
    f[BSW]  = mfaaa;
 }
 //////////////////////////////////////////////////////////////////////////
@@ -934,323 +934,323 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeFC(LBM
    //bulk viscosity
    LBMReal oP = OxxPyyPzzC;
 
-   LBMReal mfcbb = zeroReal;
-   LBMReal mfabb = zeroReal;
-   LBMReal mfbcb = zeroReal;
-   LBMReal mfbab = zeroReal;
-   LBMReal mfbbc = zeroReal;
-   LBMReal mfbba = zeroReal;
-   LBMReal mfccb = zeroReal;
-   LBMReal mfaab = zeroReal;
-   LBMReal mfcab = zeroReal;
-   LBMReal mfacb = zeroReal;
-   LBMReal mfcbc = zeroReal;
-   LBMReal mfaba = zeroReal;
-   LBMReal mfcba = zeroReal;
-   LBMReal mfabc = zeroReal;
-   LBMReal mfbcc = zeroReal;
-   LBMReal mfbaa = zeroReal;
-   LBMReal mfbca = zeroReal;
-   LBMReal mfbac = zeroReal;
-   LBMReal mfbbb = zeroReal;
-   LBMReal mfccc = zeroReal;
-   LBMReal mfaac = zeroReal;
-   LBMReal mfcac = zeroReal;
-   LBMReal mfacc = zeroReal;
-   LBMReal mfcca = zeroReal;
-   LBMReal mfaaa = zeroReal;
-   LBMReal mfcaa = zeroReal;
-   LBMReal mfaca = zeroReal;
-
-   mfaaa = press; // if drho is interpolated directly
-
+   LBMReal mfcbb = zeroReal;
+   LBMReal mfabb = zeroReal;
+   LBMReal mfbcb = zeroReal;
+   LBMReal mfbab = zeroReal;
+   LBMReal mfbbc = zeroReal;
+   LBMReal mfbba = zeroReal;
+   LBMReal mfccb = zeroReal;
+   LBMReal mfaab = zeroReal;
+   LBMReal mfcab = zeroReal;
+   LBMReal mfacb = zeroReal;
+   LBMReal mfcbc = zeroReal;
+   LBMReal mfaba = zeroReal;
+   LBMReal mfcba = zeroReal;
+   LBMReal mfabc = zeroReal;
+   LBMReal mfbcc = zeroReal;
+   LBMReal mfbaa = zeroReal;
+   LBMReal mfbca = zeroReal;
+   LBMReal mfbac = zeroReal;
+   LBMReal mfbbb = zeroReal;
+   LBMReal mfccc = zeroReal;
+   LBMReal mfaac = zeroReal;
+   LBMReal mfcac = zeroReal;
+   LBMReal mfacc = zeroReal;
+   LBMReal mfcca = zeroReal;
+   LBMReal mfaaa = zeroReal;
+   LBMReal mfcaa = zeroReal;
+   LBMReal mfaca = zeroReal;
+
+   mfaaa = press; // if drho is interpolated directly
+
    LBMReal vx1Sq = vx1*vx1;
    LBMReal vx2Sq = vx2*vx2;
    LBMReal vx3Sq = vx3*vx3;
-   LBMReal oMdrho = one;
-   //oMdrho = one - mfaaa;
-
-   //2.f
-   // linear combinations
-
+   LBMReal oMdrho = one;
+   //oMdrho = one - mfaaa;
+
+   //2.f
+   // linear combinations
+
 /////////////////////////
-   LBMReal mxxPyyPzz = mfaaa    -c2o3*(ax+by+cz)*eps_new/oP*(one+press);
-
-   LBMReal mxxMyy    = -c2o3*((ax - by)+kxxMyyAverage)*eps_new/o * (one + press);
-   LBMReal mxxMzz    = -c2o3*((ax - cz)+kxxMzzAverage)*eps_new/o * (one + press);
-
-   mfabb     = -c1o3 * ((bz + cy)+kyzAverage)*eps_new/o * (one + press);
-   mfbab     = -c1o3 * ((az + cx)+kxzAverage)*eps_new/o * (one + press);
-   mfbba     = -c1o3 * ((ay + bx)+kxyAverage)*eps_new/o * (one + press);
-
-   ////////////////////////
-   // linear combinations back
-   mfcaa = c1o3 * (mxxMyy +       mxxMzz + mxxPyyPzz);
-   mfaca = c1o3 * (-two * mxxMyy +       mxxMzz + mxxPyyPzz);
-   mfaac = c1o3 * (mxxMyy - two * mxxMzz + mxxPyyPzz);
-
-   //three
-   mfbbb = zeroReal;
-
-   LBMReal mxxyPyzz = zeroReal;
-   LBMReal mxxyMyzz = zeroReal;
-   LBMReal mxxzPyyz = zeroReal;
-   LBMReal mxxzMyyz = zeroReal;
-   LBMReal mxyyPxzz =  zeroReal;
-   LBMReal mxyyMxzz = zeroReal;
-
-   // linear combinations back
-   mfcba = (mxxyMyzz + mxxyPyzz) * c1o2;
-   mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
-   mfcab = (mxxzMyyz + mxxzPyyz) * c1o2;
-   mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
-   mfbca = (mxyyMxzz + mxyyPxzz) * c1o2;
-   mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
-
-   //4.f
-   mfacc = mfaaa*c1o9;
-   mfcac = mfacc;
-   mfcca = mfacc;
-   //5.
-
-   //6.
-   mfccc = mfaaa*c1o27;
-   ////////////////////////////////////////////////////////////////////////////////////
-   //back
-   ////////////////////////////////////////////////////////////////////////////////////
-   //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
-   ////////////////////////////////////////////////////////////////////////////////////
-   // Z - Dir
-   LBMReal m0 =  mfaac * c1o2 +      mfaab * (vx3 - c1o2) + (mfaaa + one * oMdrho) * (vx3Sq - vx3) * c1o2;
-   LBMReal m1 = -mfaac        - two * mfaab *  vx3         +  mfaaa                * (one - vx3Sq)              - one * oMdrho * vx3Sq;
-   LBMReal m2 =  mfaac * c1o2 +      mfaab * (vx3 + c1o2) + (mfaaa + one * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfaaa = m0;
-   mfaab = m1;
-   mfaac = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfabc * c1o2 +      mfabb * (vx3 - c1o2) + mfaba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfabc        - two * mfabb *  vx3         + mfaba * (one - vx3Sq);
-   m2 =  mfabc * c1o2 +      mfabb * (vx3 + c1o2) + mfaba * (vx3Sq + vx3) * c1o2;
-   mfaba = m0;
-   mfabb = m1;
-   mfabc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacc * c1o2 +      mfacb * (vx3 - c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfacc        - two * mfacb *  vx3         +  mfaca                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
-   m2 =  mfacc * c1o2 +      mfacb * (vx3 + c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfaca = m0;
-   mfacb = m1;
-   mfacc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfbac * c1o2 +      mfbab * (vx3 - c1o2) + mfbaa * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbac        - two * mfbab *  vx3         + mfbaa * (one - vx3Sq);
-   m2 =  mfbac * c1o2 +      mfbab * (vx3 + c1o2) + mfbaa * (vx3Sq + vx3) * c1o2;
-   mfbaa = m0;
-   mfbab = m1;
-   mfbac = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbbc * c1o2 +      mfbbb * (vx3 - c1o2) + mfbba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbbc        - two * mfbbb *  vx3         + mfbba * (one - vx3Sq);
-   m2 =  mfbbc * c1o2 +      mfbbb * (vx3 + c1o2) + mfbba * (vx3Sq + vx3) * c1o2;
-   mfbba = m0;
-   mfbbb = m1;
-   mfbbc = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcc * c1o2 +      mfbcb * (vx3 - c1o2) + mfbca * (vx3Sq - vx3) * c1o2;
-   m1 = -mfbcc        - two * mfbcb *  vx3         + mfbca * (one - vx3Sq);
-   m2 =  mfbcc * c1o2 +      mfbcb * (vx3 + c1o2) + mfbca * (vx3Sq + vx3) * c1o2;
-   mfbca = m0;
-   mfbcb = m1;
-   mfbcc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcac * c1o2 +      mfcab * (vx3 - c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfcac        - two * mfcab *  vx3         +  mfcaa                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
-   m2 =  mfcac * c1o2 +      mfcab * (vx3 + c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfcaa = m0;
-   mfcab = m1;
-   mfcac = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbc * c1o2 +      mfcbb * (vx3 - c1o2) + mfcba * (vx3Sq - vx3) * c1o2;
-   m1 = -mfcbc        - two * mfcbb *  vx3         + mfcba * (one - vx3Sq);
-   m2 =  mfcbc * c1o2 +      mfcbb * (vx3 + c1o2) + mfcba * (vx3Sq + vx3) * c1o2;
-   mfcba = m0;
-   mfcbb = m1;
-   mfcbc = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfccb * (vx3 - c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq - vx3) * c1o2;
-   m1 = -mfccc        - two * mfccb *  vx3         +  mfcca                  * (one - vx3Sq)              - c1o9 * oMdrho * vx3Sq;
-   m2 =  mfccc * c1o2 +      mfccb * (vx3 + c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq + vx3) * c1o2;
-   mfcca = m0;
-   mfccb = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18   Konditionieren
-   ////////////////////////////////////////////////////////////////////////////////////
-   // Y - Dir
-   m0 =  mfaca * c1o2 +      mfaba * (vx2 - c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfaca        - two * mfaba *  vx2         +  mfaaa                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
-   m2 =  mfaca * c1o2 +      mfaba * (vx2 + c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaaa = m0;
-   mfaba = m1;
-   mfaca = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacb * c1o2 +      mfabb * (vx2 - c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfacb        - two * mfabb *  vx2         +  mfaab                  * (one - vx2Sq)              - c2o3 * oMdrho * vx2Sq;
-   m2 =  mfacb * c1o2 +      mfabb * (vx2 + c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaab = m0;
-   mfabb = m1;
-   mfacb = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfacc * c1o2 +      mfabc * (vx2 - c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfacc        - two * mfabc *  vx2         +  mfaac                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
-   m2 =  mfacc * c1o2 +      mfabc * (vx2 + c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfaac = m0;
-   mfabc = m1;
-   mfacc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfbca * c1o2 +      mfbba * (vx2 - c1o2) + mfbaa * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbca        - two * mfbba *  vx2         + mfbaa * (one - vx2Sq);
-   m2 =  mfbca * c1o2 +      mfbba * (vx2 + c1o2) + mfbaa * (vx2Sq + vx2) * c1o2;
-   mfbaa = m0;
-   mfbba = m1;
-   mfbca = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcb * c1o2 +      mfbbb * (vx2 - c1o2) + mfbab * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbcb        - two * mfbbb *  vx2         + mfbab * (one - vx2Sq);
-   m2 =  mfbcb * c1o2 +      mfbbb * (vx2 + c1o2) + mfbab * (vx2Sq + vx2) * c1o2;
-   mfbab = m0;
-   mfbbb = m1;
-   mfbcb = m2;
-   /////////b//////////////////////////////////////////////////////////////////////////
-   m0 =  mfbcc * c1o2 +      mfbbc * (vx2 - c1o2) + mfbac * (vx2Sq - vx2) * c1o2;
-   m1 = -mfbcc        - two * mfbbc *  vx2         + mfbac * (one - vx2Sq);
-   m2 =  mfbcc * c1o2 +      mfbbc * (vx2 + c1o2) + mfbac * (vx2Sq + vx2) * c1o2;
-   mfbac = m0;
-   mfbbc = m1;
-   mfbcc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcca * c1o2 +      mfcba * (vx2 - c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfcca        - two * mfcba *  vx2         +  mfcaa                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
-   m2 =  mfcca * c1o2 +      mfcba * (vx2 + c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcaa = m0;
-   mfcba = m1;
-   mfcca = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccb * c1o2 +      mfcbb * (vx2 - c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfccb        - two * mfcbb *  vx2         +  mfcab                  * (one - vx2Sq)              - c2o9 * oMdrho * vx2Sq;
-   m2 =  mfccb * c1o2 +      mfcbb * (vx2 + c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcab = m0;
-   mfcbb = m1;
-   mfccb = m2;
-   /////////c//////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfcbc * (vx2 - c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
-   m1 = -mfccc        - two * mfcbc *  vx2         +  mfcac                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
-   m2 =  mfccc * c1o2 +      mfcbc * (vx2 + c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
-   mfcac = m0;
-   mfcbc = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
-   ////////////////////////////////////////////////////////////////////////////////////
-   // X - Dir
-   m0 =  mfcaa * c1o2 +      mfbaa * (vx1 - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcaa        - two * mfbaa *  vx1         +  mfaaa                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcaa * c1o2 +      mfbaa * (vx1 + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaaa = m0;
-   mfbaa = m1;
-   mfcaa = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcba * c1o2 +      mfbba * (vx1 - c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcba        - two * mfbba *  vx1         +  mfaba                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcba * c1o2 +      mfbba * (vx1 + c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaba = m0;
-   mfbba = m1;
-   mfcba = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcca * c1o2 +      mfbca * (vx1 - c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcca        - two * mfbca *  vx1         +  mfaca                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcca * c1o2 +      mfbca * (vx1 + c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaca = m0;
-   mfbca = m1;
-   mfcca = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcab * c1o2 +      mfbab * (vx1 - c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcab        - two * mfbab *  vx1         +  mfaab                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcab * c1o2 +      mfbab * (vx1 + c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaab = m0;
-   mfbab = m1;
-   mfcab = m2;
-   ///////////b////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbb * c1o2 +      mfbbb * (vx1 - c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcbb        - two * mfbbb *  vx1         +  mfabb                  * (one - vx1Sq)              - c4o9 * oMdrho * vx1Sq;
-   m2 =  mfcbb * c1o2 +      mfbbb * (vx1 + c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfabb = m0;
-   mfbbb = m1;
-   mfcbb = m2;
-   ///////////b////////////////////////////////////////////////////////////////////////
-   m0 =  mfccb * c1o2 +      mfbcb * (vx1 - c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfccb        - two * mfbcb *  vx1         +  mfacb                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfccb * c1o2 +      mfbcb * (vx1 + c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfacb = m0;
-   mfbcb = m1;
-   mfccb = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-   ////////////////////////////////////////////////////////////////////////////////////
-   m0 =  mfcac * c1o2 +      mfbac * (vx1 - c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcac        - two * mfbac *  vx1         +  mfaac                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfcac * c1o2 +      mfbac * (vx1 + c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfaac = m0;
-   mfbac = m1;
-   mfcac = m2;
-   ///////////c////////////////////////////////////////////////////////////////////////
-   m0 =  mfcbc * c1o2 +      mfbbc * (vx1 - c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfcbc        - two * mfbbc *  vx1         +  mfabc                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
-   m2 =  mfcbc * c1o2 +      mfbbc * (vx1 + c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfabc = m0;
-   mfbbc = m1;
-   mfcbc = m2;
-   ///////////c////////////////////////////////////////////////////////////////////////
-   m0 =  mfccc * c1o2 +      mfbcc * (vx1 - c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
-   m1 = -mfccc        - two * mfbcc *  vx1         +  mfacc                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
-   m2 =  mfccc * c1o2 +      mfbcc * (vx1 + c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
-   mfacc = m0;
-   mfbcc = m1;
-   mfccc = m2;
-   ////////////////////////////////////////////////////////////////////////////////////
-
-   f[E]    = mfcbb;
-   f[W]    = mfabb;
-   f[N]    = mfbcb;
-   f[S]    = mfbab;
-   f[T]    = mfbbc;
-   f[B]    = mfbba;
-   f[NE]   = mfccb;
-   f[SW]   = mfaab;
-   f[SE]   = mfcab;
-   f[NW]   = mfacb;
-   f[TE]   = mfcbc;
-   f[BW]   = mfaba;
-   f[BE]   = mfcba;
-   f[TW]   = mfabc;
-   f[TN]   = mfbcc;
-   f[BS]   = mfbaa;
-   f[BN]   = mfbca;
-   f[TS]   = mfbac;
-   f[ZERO] = mfbbb;
-   f[TNE]  = mfccc;
-   f[TSE]  = mfcac;
-   f[BNE]  = mfcca;
-   f[BSE]  = mfcaa;
-   f[TNW]  = mfacc;
-   f[TSW]  = mfaac;
-   f[BNW]  = mfaca;
+   LBMReal mxxPyyPzz = mfaaa    -c2o3*(ax+by+cz)*eps_new/oP*(one+press);
+
+   LBMReal mxxMyy    = -c2o3*((ax - by)+kxxMyyAverage)*eps_new/o * (one + press);
+   LBMReal mxxMzz    = -c2o3*((ax - cz)+kxxMzzAverage)*eps_new/o * (one + press);
+
+   mfabb     = -c1o3 * ((bz + cy)+kyzAverage)*eps_new/o * (one + press);
+   mfbab     = -c1o3 * ((az + cx)+kxzAverage)*eps_new/o * (one + press);
+   mfbba     = -c1o3 * ((ay + bx)+kxyAverage)*eps_new/o * (one + press);
+
+   ////////////////////////
+   // linear combinations back
+   mfcaa = c1o3 * (mxxMyy +       mxxMzz + mxxPyyPzz);
+   mfaca = c1o3 * (-two * mxxMyy +       mxxMzz + mxxPyyPzz);
+   mfaac = c1o3 * (mxxMyy - two * mxxMzz + mxxPyyPzz);
+
+   //three
+   mfbbb = zeroReal;
+
+   LBMReal mxxyPyzz = zeroReal;
+   LBMReal mxxyMyzz = zeroReal;
+   LBMReal mxxzPyyz = zeroReal;
+   LBMReal mxxzMyyz = zeroReal;
+   LBMReal mxyyPxzz =  zeroReal;
+   LBMReal mxyyMxzz = zeroReal;
+
+   // linear combinations back
+   mfcba = (mxxyMyzz + mxxyPyzz) * c1o2;
+   mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
+   mfcab = (mxxzMyyz + mxxzPyyz) * c1o2;
+   mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
+   mfbca = (mxyyMxzz + mxyyPxzz) * c1o2;
+   mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
+
+   //4.f
+   mfacc = mfaaa*c1o9;
+   mfcac = mfacc;
+   mfcca = mfacc;
+   //5.
+
+   //6.
+   mfccc = mfaaa*c1o27;
+   ////////////////////////////////////////////////////////////////////////////////////
+   //back
+   ////////////////////////////////////////////////////////////////////////////////////
+   //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
+   ////////////////////////////////////////////////////////////////////////////////////
+   // Z - Dir
+   LBMReal m0 =  mfaac * c1o2 +      mfaab * (vx3 - c1o2) + (mfaaa + one * oMdrho) * (vx3Sq - vx3) * c1o2;
+   LBMReal m1 = -mfaac        - two * mfaab *  vx3         +  mfaaa                * (one - vx3Sq)              - one * oMdrho * vx3Sq;
+   LBMReal m2 =  mfaac * c1o2 +      mfaab * (vx3 + c1o2) + (mfaaa + one * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfaaa = m0;
+   mfaab = m1;
+   mfaac = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfabc * c1o2 +      mfabb * (vx3 - c1o2) + mfaba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfabc        - two * mfabb *  vx3         + mfaba * (one - vx3Sq);
+   m2 =  mfabc * c1o2 +      mfabb * (vx3 + c1o2) + mfaba * (vx3Sq + vx3) * c1o2;
+   mfaba = m0;
+   mfabb = m1;
+   mfabc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacc * c1o2 +      mfacb * (vx3 - c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfacc        - two * mfacb *  vx3         +  mfaca                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
+   m2 =  mfacc * c1o2 +      mfacb * (vx3 + c1o2) + (mfaca + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfaca = m0;
+   mfacb = m1;
+   mfacc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfbac * c1o2 +      mfbab * (vx3 - c1o2) + mfbaa * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbac        - two * mfbab *  vx3         + mfbaa * (one - vx3Sq);
+   m2 =  mfbac * c1o2 +      mfbab * (vx3 + c1o2) + mfbaa * (vx3Sq + vx3) * c1o2;
+   mfbaa = m0;
+   mfbab = m1;
+   mfbac = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbbc * c1o2 +      mfbbb * (vx3 - c1o2) + mfbba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbbc        - two * mfbbb *  vx3         + mfbba * (one - vx3Sq);
+   m2 =  mfbbc * c1o2 +      mfbbb * (vx3 + c1o2) + mfbba * (vx3Sq + vx3) * c1o2;
+   mfbba = m0;
+   mfbbb = m1;
+   mfbbc = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcc * c1o2 +      mfbcb * (vx3 - c1o2) + mfbca * (vx3Sq - vx3) * c1o2;
+   m1 = -mfbcc        - two * mfbcb *  vx3         + mfbca * (one - vx3Sq);
+   m2 =  mfbcc * c1o2 +      mfbcb * (vx3 + c1o2) + mfbca * (vx3Sq + vx3) * c1o2;
+   mfbca = m0;
+   mfbcb = m1;
+   mfbcc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcac * c1o2 +      mfcab * (vx3 - c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfcac        - two * mfcab *  vx3         +  mfcaa                  * (one - vx3Sq)              - c1o3 * oMdrho * vx3Sq;
+   m2 =  mfcac * c1o2 +      mfcab * (vx3 + c1o2) + (mfcaa + c1o3 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfcaa = m0;
+   mfcab = m1;
+   mfcac = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbc * c1o2 +      mfcbb * (vx3 - c1o2) + mfcba * (vx3Sq - vx3) * c1o2;
+   m1 = -mfcbc        - two * mfcbb *  vx3         + mfcba * (one - vx3Sq);
+   m2 =  mfcbc * c1o2 +      mfcbb * (vx3 + c1o2) + mfcba * (vx3Sq + vx3) * c1o2;
+   mfcba = m0;
+   mfcbb = m1;
+   mfcbc = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfccb * (vx3 - c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq - vx3) * c1o2;
+   m1 = -mfccc        - two * mfccb *  vx3         +  mfcca                  * (one - vx3Sq)              - c1o9 * oMdrho * vx3Sq;
+   m2 =  mfccc * c1o2 +      mfccb * (vx3 + c1o2) + (mfcca + c1o9 * oMdrho) * (vx3Sq + vx3) * c1o2;
+   mfcca = m0;
+   mfccb = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18   Konditionieren
+   ////////////////////////////////////////////////////////////////////////////////////
+   // Y - Dir
+   m0 =  mfaca * c1o2 +      mfaba * (vx2 - c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfaca        - two * mfaba *  vx2         +  mfaaa                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
+   m2 =  mfaca * c1o2 +      mfaba * (vx2 + c1o2) + (mfaaa + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaaa = m0;
+   mfaba = m1;
+   mfaca = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacb * c1o2 +      mfabb * (vx2 - c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfacb        - two * mfabb *  vx2         +  mfaab                  * (one - vx2Sq)              - c2o3 * oMdrho * vx2Sq;
+   m2 =  mfacb * c1o2 +      mfabb * (vx2 + c1o2) + (mfaab + c2o3 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaab = m0;
+   mfabb = m1;
+   mfacb = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfacc * c1o2 +      mfabc * (vx2 - c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfacc        - two * mfabc *  vx2         +  mfaac                  * (one - vx2Sq)              - c1o6 * oMdrho * vx2Sq;
+   m2 =  mfacc * c1o2 +      mfabc * (vx2 + c1o2) + (mfaac + c1o6 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfaac = m0;
+   mfabc = m1;
+   mfacc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfbca * c1o2 +      mfbba * (vx2 - c1o2) + mfbaa * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbca        - two * mfbba *  vx2         + mfbaa * (one - vx2Sq);
+   m2 =  mfbca * c1o2 +      mfbba * (vx2 + c1o2) + mfbaa * (vx2Sq + vx2) * c1o2;
+   mfbaa = m0;
+   mfbba = m1;
+   mfbca = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcb * c1o2 +      mfbbb * (vx2 - c1o2) + mfbab * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbcb        - two * mfbbb *  vx2         + mfbab * (one - vx2Sq);
+   m2 =  mfbcb * c1o2 +      mfbbb * (vx2 + c1o2) + mfbab * (vx2Sq + vx2) * c1o2;
+   mfbab = m0;
+   mfbbb = m1;
+   mfbcb = m2;
+   /////////b//////////////////////////////////////////////////////////////////////////
+   m0 =  mfbcc * c1o2 +      mfbbc * (vx2 - c1o2) + mfbac * (vx2Sq - vx2) * c1o2;
+   m1 = -mfbcc        - two * mfbbc *  vx2         + mfbac * (one - vx2Sq);
+   m2 =  mfbcc * c1o2 +      mfbbc * (vx2 + c1o2) + mfbac * (vx2Sq + vx2) * c1o2;
+   mfbac = m0;
+   mfbbc = m1;
+   mfbcc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcca * c1o2 +      mfcba * (vx2 - c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfcca        - two * mfcba *  vx2         +  mfcaa                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
+   m2 =  mfcca * c1o2 +      mfcba * (vx2 + c1o2) + (mfcaa + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcaa = m0;
+   mfcba = m1;
+   mfcca = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccb * c1o2 +      mfcbb * (vx2 - c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfccb        - two * mfcbb *  vx2         +  mfcab                  * (one - vx2Sq)              - c2o9 * oMdrho * vx2Sq;
+   m2 =  mfccb * c1o2 +      mfcbb * (vx2 + c1o2) + (mfcab + c2o9 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcab = m0;
+   mfcbb = m1;
+   mfccb = m2;
+   /////////c//////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfcbc * (vx2 - c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq - vx2) * c1o2;
+   m1 = -mfccc        - two * mfcbc *  vx2         +  mfcac                   * (one - vx2Sq)              - c1o18 * oMdrho * vx2Sq;
+   m2 =  mfccc * c1o2 +      mfcbc * (vx2 + c1o2) + (mfcac + c1o18 * oMdrho) * (vx2Sq + vx2) * c1o2;
+   mfcac = m0;
+   mfcbc = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+   ////////////////////////////////////////////////////////////////////////////////////
+   // X - Dir
+   m0 =  mfcaa * c1o2 +      mfbaa * (vx1 - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcaa        - two * mfbaa *  vx1         +  mfaaa                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcaa * c1o2 +      mfbaa * (vx1 + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaaa = m0;
+   mfbaa = m1;
+   mfcaa = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcba * c1o2 +      mfbba * (vx1 - c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcba        - two * mfbba *  vx1         +  mfaba                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcba * c1o2 +      mfbba * (vx1 + c1o2) + (mfaba + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaba = m0;
+   mfbba = m1;
+   mfcba = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcca * c1o2 +      mfbca * (vx1 - c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcca        - two * mfbca *  vx1         +  mfaca                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcca * c1o2 +      mfbca * (vx1 + c1o2) + (mfaca + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaca = m0;
+   mfbca = m1;
+   mfcca = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcab * c1o2 +      mfbab * (vx1 - c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcab        - two * mfbab *  vx1         +  mfaab                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcab * c1o2 +      mfbab * (vx1 + c1o2) + (mfaab + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaab = m0;
+   mfbab = m1;
+   mfcab = m2;
+   ///////////b////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbb * c1o2 +      mfbbb * (vx1 - c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcbb        - two * mfbbb *  vx1         +  mfabb                  * (one - vx1Sq)              - c4o9 * oMdrho * vx1Sq;
+   m2 =  mfcbb * c1o2 +      mfbbb * (vx1 + c1o2) + (mfabb + c4o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfabb = m0;
+   mfbbb = m1;
+   mfcbb = m2;
+   ///////////b////////////////////////////////////////////////////////////////////////
+   m0 =  mfccb * c1o2 +      mfbcb * (vx1 - c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfccb        - two * mfbcb *  vx1         +  mfacb                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfccb * c1o2 +      mfbcb * (vx1 + c1o2) + (mfacb + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfacb = m0;
+   mfbcb = m1;
+   mfccb = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+   ////////////////////////////////////////////////////////////////////////////////////
+   m0 =  mfcac * c1o2 +      mfbac * (vx1 - c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcac        - two * mfbac *  vx1         +  mfaac                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfcac * c1o2 +      mfbac * (vx1 + c1o2) + (mfaac + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfaac = m0;
+   mfbac = m1;
+   mfcac = m2;
+   ///////////c////////////////////////////////////////////////////////////////////////
+   m0 =  mfcbc * c1o2 +      mfbbc * (vx1 - c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfcbc        - two * mfbbc *  vx1         +  mfabc                  * (one - vx1Sq)              - c1o9 * oMdrho * vx1Sq;
+   m2 =  mfcbc * c1o2 +      mfbbc * (vx1 + c1o2) + (mfabc + c1o9 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfabc = m0;
+   mfbbc = m1;
+   mfcbc = m2;
+   ///////////c////////////////////////////////////////////////////////////////////////
+   m0 =  mfccc * c1o2 +      mfbcc * (vx1 - c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq - vx1) * c1o2;
+   m1 = -mfccc        - two * mfbcc *  vx1         +  mfacc                   * (one - vx1Sq)              - c1o36 * oMdrho * vx1Sq;
+   m2 =  mfccc * c1o2 +      mfbcc * (vx1 + c1o2) + (mfacc + c1o36 * oMdrho) * (vx1Sq + vx1) * c1o2;
+   mfacc = m0;
+   mfbcc = m1;
+   mfccc = m2;
+   ////////////////////////////////////////////////////////////////////////////////////
+
+   f[E]    = mfcbb;
+   f[W]    = mfabb;
+   f[N]    = mfbcb;
+   f[S]    = mfbab;
+   f[T]    = mfbbc;
+   f[B]    = mfbba;
+   f[NE]   = mfccb;
+   f[SW]   = mfaab;
+   f[SE]   = mfcab;
+   f[NW]   = mfacb;
+   f[TE]   = mfcbc;
+   f[BW]   = mfaba;
+   f[BE]   = mfcba;
+   f[TW]   = mfabc;
+   f[TN]   = mfbcc;
+   f[BS]   = mfbaa;
+   f[BN]   = mfbca;
+   f[TS]   = mfbac;
+   f[ZERO] = mfbbb;
+   f[TNE]  = mfccc;
+   f[TSE]  = mfcac;
+   f[BNE]  = mfcca;
+   f[BSE]  = mfcaa;
+   f[TNW]  = mfacc;
+   f[TSW]  = mfaac;
+   f[BNW]  = mfaca;
    f[BSW]  = mfaaa;
 }
 //////////////////////////////////////////////////////////////////////////
diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
index 0249f9f64..77cfa4d8d 100644
--- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
@@ -4,7 +4,10 @@
 
 using namespace std;
 
-SpongeLayerBlockVisitor::SpongeLayerBlockVisitor(GbCuboid3DPtr boundingBox) : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), boundingBox(boundingBox)
+SpongeLayerBlockVisitor::SpongeLayerBlockVisitor(GbCuboid3DPtr boundingBox, LBMReal collFactor) : 
+   Block3DVisitor(0, Grid3DSystem::MAXLEVEL), 
+   boundingBox(boundingBox),
+   collFactor(collFactor)
 {
 
 }
@@ -31,7 +34,8 @@ void SpongeLayerBlockVisitor::visit(Grid3DPtr grid, Block3DPtr block)
       if (boundingBox->isCellInsideGbObject3D(minX1, minX2, minX3, maxX1, maxX2, maxX3))
       {
          LBMKernelPtr kernel = block->getKernel();
-         kernel->setCollisionFactor(LBMSystem::calcCollisionFactor(0.01, block->getLevel()));
+         //kernel->setCollisionFactor(LBMSystem::calcCollisionFactor(0.01, block->getLevel()));
+         kernel->setCollisionFactor(collFactor);
       }
    }
 }
diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h
index fc7c6d328..621c6ffc3 100644
--- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h
+++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.h
@@ -12,12 +12,13 @@
 class SpongeLayerBlockVisitor : public Block3DVisitor
 {
 public:
-   SpongeLayerBlockVisitor(GbCuboid3DPtr boundingBox);
+   SpongeLayerBlockVisitor(GbCuboid3DPtr boundingBox, LBMReal collFactor);
    virtual ~SpongeLayerBlockVisitor();
    virtual void visit(Grid3DPtr grid, Block3DPtr block);
 protected:
 private:
    GbCuboid3DPtr boundingBox;
+   LBMReal collFactor;
 };
 
 #endif // SetSpongeLayerBlockVisitor_h__
-- 
GitLab