diff --git a/.gitignore b/.gitignore
index a63828c1c198e6d9685162e21a6ade6128f0a0d2..bb7d38f1349cff7e73d0a8881d1c9dbf3b90f757 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,8 +3,9 @@ build/
 bin/
 cmake-build-debug/
 run/
-buildGCC
-buildWSL
+buildGCC/
+buildWSL/
+buildCygwin/
 pybuild/
 
 # Python
diff --git a/apps/cpu/ConcreteExtrusion/ConcreteExtrusion.cpp b/apps/cpu/ConcreteExtrusion/ConcreteExtrusion.cpp
index 1407d8c4d3f564f6bff351e7450199d983391622..e73ff544b59991bf95874e058e3df3b7113370a0 100644
--- a/apps/cpu/ConcreteExtrusion/ConcreteExtrusion.cpp
+++ b/apps/cpu/ConcreteExtrusion/ConcreteExtrusion.cpp
@@ -144,7 +144,7 @@ void run(string configname)
         real nu_l = nu_h;
         nu_h *= 0.1;
 
-        real rho_h_LB = 1;
+        //real rho_h_LB = 1;
         // surface tension
         real sigma_LB = 0.0; //rho_h_LB *U_LB *U_LB *D_LB / We;
 
@@ -155,12 +155,12 @@ void run(string configname)
         real beta = 12.0 * sigma_LB / interfaceWidth;
         real kappa = 1.5 * interfaceWidth * sigma_LB;
 
-        double tau0 = 715.218181094648*1000.; // Pa
-        double muConcrete = 2.1133054011798826; // [Pa s]
-        real u = Uo; //[m/s]
+        //double tau0 = 715.218181094648*1000.; // Pa
+        //double muConcrete = 2.1133054011798826; // [Pa s]
+        //real u = Uo; //[m/s]
 
 
-        double Bm = (tau0 * D) / (muConcrete * u);
+        //double Bm = (tau0 * D) / (muConcrete * u);
         double tau0_LB = 0.02;
         //Bm *nu_h *U_LB / (D / dx);
 
diff --git a/apps/cpu/JetBreakup/JetBreakup.cpp b/apps/cpu/JetBreakup/JetBreakup.cpp
index f4a3bdc26ac6abc724aded591d691e77b2a50faa..de2f5b95035f8a2b786848941a897ad65cb9cf61 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cpp
+++ b/apps/cpu/JetBreakup/JetBreakup.cpp
@@ -147,7 +147,7 @@ void run(string configname)
         real nu_l = nu_h;
         nu_h *= 0.1;
 
-        real rho_h_LB = 1;
+        //real rho_h_LB = 1;
         // surface tension
         real sigma_LB = 0.0; //rho_h_LB *U_LB *U_LB *D_LB / We;
 
@@ -158,12 +158,12 @@ void run(string configname)
         real beta = 12.0 * sigma_LB / interfaceWidth;
         real kappa = 1.5 * interfaceWidth * sigma_LB;
 
-        double tau0 = 715.218181094648*1000.; // Pa
-        double muConcrete = 2.1133054011798826; // [Pa s]
-        real u = Uo; //[m/s]
+        //double tau0 = 715.218181094648*1000.; // Pa
+        //double muConcrete = 2.1133054011798826; // [Pa s]
+        //real u = Uo; //[m/s]
 
 
-        double Bm = (tau0 * D) / (muConcrete * u);
+        //double Bm = (tau0 * D) / (muConcrete * u);
         double tau0_LB = 0;
         //0.02;
         //Bm *nu_h *U_LB / (D / dx);
diff --git a/apps/cpu/ShotcreteJet/jet.cpp b/apps/cpu/ShotcreteJet/jet.cpp
index 18a018f81c598575313336c4c95a977d2b934b55..baaa572353ecfa91a14c12ab0f7c5ae1bcfd97c1 100644
--- a/apps/cpu/ShotcreteJet/jet.cpp
+++ b/apps/cpu/ShotcreteJet/jet.cpp
@@ -91,17 +91,17 @@ int main(int argc, char *argv[])
         //double g_maxX2 = 0.401582;
         //double g_maxX3 = 0.175;//0.21;
 
-        double dx = 1e-3;
+        double dx = 1;
 
-        double g_maxX3_box = -0.065;        
+        double g_maxX3_box = -0.065e3;        
 
-        double g_minX1 = -1.49631;
-        double g_minX2 = 0.193582;
-        double g_minX3 = g_maxX3_box - 0.03;//-0.095; //-0.215; 
+        double g_minX1 = -1.49631e3;
+        double g_minX2 = 0.193582e3;
+        double g_minX3 = g_maxX3_box - 0.03e3;//-0.095; //-0.215; 
 
-        double g_maxX1 = -1.10631;
-        double g_maxX2 = 0.583582;
-        double g_maxX3 = 0.175;
+        double g_maxX1 = -1.10631e3;
+        double g_maxX2 = 0.583582e3;
+        double g_maxX3 = 0.175e3;
 
       
 
@@ -353,9 +353,9 @@ int main(int argc, char *argv[])
         // SPtr<BC> inflowAirBC1(new VelocityBC(true, false, false, fct, 0, BCFunction::INFCONST));
         // inflowAirBC1->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
         // t = U * sin(alpha * _pi / 180) * (1 - (((((x1 - x0) ^ 2 + (x2 - y0) ^ 2 + (x3 - z0) ^ 2) ^ 0.5) / R) ^ NplusOne));
-        cx1 = -1.31416;
-        cx2 = 0.388684;
-        cx3 = 0.138177;
+        cx1 = -1.31416e3;
+        cx2 = 0.388684e3;
+        cx3 = 0.138177e3;
         alpha = 0;
         gamma = 225;
         fctVx1.DefineConst("x0", cx1);
@@ -398,9 +398,9 @@ int main(int argc, char *argv[])
             VF_LOG_INFO("NplusOne = {}", N + 1.0);
             // return 0;
         }
-        cx1 = -1.31303;
-        cx2 = 0.377234;
-        cx3 = 0.138174;
+        cx1 = -1.31303e3;
+        cx2 = 0.377234e3;
+        cx3 = 0.138174e3;
         alpha = 60;
         fctVx1.DefineConst("x0", cx1);
         fctVx1.DefineConst("y0", cx2);
@@ -418,9 +418,9 @@ int main(int argc, char *argv[])
         SPtr<BC> inflowAirBC2(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
         inflowAirBC2->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
 
-        cx1 = -1.2948374155694822;
-        cx2 = 0.37733728717266285;
-        cx3 = 0.13840460401111598;
+        cx1 = -1.2948374155694822e3;
+        cx2 = 0.37733728717266285e3;
+        cx3 = 0.13840460401111598e3;
         alpha = 120;
         fctVx1.DefineConst("x0", cx1);
         fctVx1.DefineConst("y0", cx2);
@@ -438,9 +438,9 @@ int main(int argc, char *argv[])
         SPtr<BC> inflowAirBC3(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
         inflowAirBC3->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
 
-        cx1 = -1.28847;
-        cx2 = 0.3885;
-        cx3 = 0.1385;
+        cx1 = -1.28847e3;
+        cx2 = 0.3885e3;
+        cx3 = 0.1385e3;
         alpha = 180;
         fctVx1.DefineConst("x0", cx1);
         fctVx1.DefineConst("y0", cx2);
@@ -458,9 +458,9 @@ int main(int argc, char *argv[])
         SPtr<BC> inflowAirBC4(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
         inflowAirBC4->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
 
-        cx1 = -1.294771417778694;
-        cx2 = 0.399787947463142;
-        cx3 = 0.1383429692754194;
+        cx1 = -1.294771417778694e3;
+        cx2 = 0.399787947463142e3;
+        cx3 = 0.1383429692754194e3;
         alpha = 240;
         fctVx1.DefineConst("x0", cx1);
         fctVx1.DefineConst("y0", cx2);
@@ -478,9 +478,9 @@ int main(int argc, char *argv[])
         SPtr<BC> inflowAirBC5(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
         inflowAirBC5->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
 
-        cx1 = -1.3077338898450492;
-        cx2 = 0.3998516560596088;
-        cx3 = 0.13843501416896437;
+        cx1 = -1.3077338898450492e3;
+        cx2 = 0.3998516560596088e3;
+        cx3 = 0.13843501416896437e3;
         alpha = 300;
         fctVx1.DefineConst("x0", cx1);
         fctVx1.DefineConst("y0", cx2);
@@ -629,7 +629,7 @@ int main(int argc, char *argv[])
         ///////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshNozzleAirDistributor = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirDistributor:start");
-        meshNozzleAirDistributor->readMeshFromSTLFileASCII(geoPath + "/01_Nozzle_Air_Distributor.stl", false);
+        meshNozzleAirDistributor->readMeshFromSTLFileBinary(geoPath + "/01_Nozzle_Air_Distributor.stl", false);
         if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirDistributor:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleAirDistributor.get(), outputPath + "/geo/meshNozzleAirDistributor", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intrNozzleAirDistributor = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleAirDistributor, grid, noSlipBC, Interactor3D::SOLID, Interactor3D::POINTS);
@@ -682,14 +682,14 @@ int main(int argc, char *argv[])
         ///////////////////////////////////////////////////////////
         // inflow
         //GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.20105, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.23, 0.013));
-        GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3 - 2.0 * dx, -1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3+2.0*dx, 0.013));
+        GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_maxX3 - 2.0 * dx, -1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_maxX3+2.0*dx, 0.013e3));
         if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), outputPath + "/geo/geoInflow", WbWriterVtkXmlBinary::getInstance());
         SPtr<D3Q27Interactor> intrInflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, inflowConcreteBC, Interactor3D::SOLID));
         ///////////////////////////////////////////////////////////
         // outflow
         //GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, -0.22, -1.30181 + 0.0005, 0.390872 - 0.00229, -0.21, 0.013));
         //GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, -0.426, -1.30181 + 0.0005, 0.390872 - 0.00229, -0.415, 0.013));
-        GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3, -1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3+2.*dx, 0.013));
+        GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_minX3, -1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_minX3+2.*dx, 0.013e3));
         if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), outputPath + "/geo/geoOutflow", WbWriterVtkXmlBinary::getInstance());
         SPtr<D3Q27Interactor> intrOutflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowBC, Interactor3D::SOLID));
         ///////////////////////////////////////////////////////////
@@ -702,7 +702,7 @@ int main(int argc, char *argv[])
         /////////////////////////////////////////////////////////////
         // Fluid area
         //GbCylinder3DPtr geoFluidArea(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3, -1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3, 0.013));
-        GbCylinder3DPtr geoFluidArea(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3_box, -1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3, 0.013));
+        GbCylinder3DPtr geoFluidArea(new GbCylinder3D(-1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_maxX3_box, -1.30181e3 + 0.0005e3, 0.390872e3 - 0.00229e3, g_maxX3, 0.013e3));
         if (myid == 0) GbSystem3D::writeGeoObject(geoFluidArea.get(), outputPath + "/geo/geoFluidArea", WbWriterVtkXmlBinary::getInstance());
         SPtr<D3Q27Interactor> intrFluidArea = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoFluidArea, grid, noSlipBC, Interactor3D::INVERSESOLID));
 
@@ -772,42 +772,42 @@ int main(int argc, char *argv[])
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet1 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet1:start");
-        meshAirInlet1->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_1.stl", true);
+        meshAirInlet1->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_1.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet1:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet1.get(), outputPath + "/geo/meshAirInlet1", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet1 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet1, grid, inflowAirBC1, Interactor3D::SOLID, Interactor3D::POINTS);
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet2 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet2:start");
-        meshAirInlet2->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_2.stl", true);
+        meshAirInlet2->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_2.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet2:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet2.get(), outputPath + "/geo/meshAirInlet2", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet2 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet2, grid, inflowAirBC2, Interactor3D::SOLID, Interactor3D::POINTS);
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet3 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet3:start");
-        meshAirInlet3->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_3.stl", true);
+        meshAirInlet3->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_3.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet3:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet3.get(), outputPath + "/geo/meshAirInlet3", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet3 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet3, grid, inflowAirBC3, Interactor3D::SOLID, Interactor3D::POINTS);
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet4 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet4:start");
-        meshAirInlet4->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_4.stl", true);
+        meshAirInlet4->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_4.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet4:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet4.get(), outputPath + "/geo/meshAirInlet4", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet4 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet4, grid, inflowAirBC4, Interactor3D::SOLID, Interactor3D::POINTS);
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet5 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet5:start");
-        meshAirInlet5->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_5.stl", true);
+        meshAirInlet5->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_5.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet5:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet5.get(), outputPath + "/geo/meshAirInlet5", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet5 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet5, grid, inflowAirBC5, Interactor3D::SOLID, Interactor3D::POINTS);
         ///////////////////////////////////////////////////////////
         SPtr<GbTriFaceMesh3D> meshAirInlet6 = std::make_shared<GbTriFaceMesh3D>();
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet6:start");
-        meshAirInlet6->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_6.stl", true);
+        meshAirInlet6->readMeshFromSTLFileBinary(geoPath + "/Air_Inlet_6.stl", true);
         if (myid == 0) UBLOG(logINFO, "Read meshAirInlet6:end");
         if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet6.get(), outputPath + "/geo/meshAirInlet6", WbWriterVtkXmlBinary::getInstance());
         SPtr<Interactor3D> intAirInlet6 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet6, grid, inflowAirBC6, Interactor3D::SOLID, Interactor3D::POINTS);
@@ -993,7 +993,7 @@ int main(int argc, char *argv[])
             fct1.DefineConst("interfaceThickness", interfaceThickness * dx);
 
             MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
-            initVisitor.setPhi(fct1);
+            //initVisitor.setPhi(fct1);
             grid->accept(initVisitor);
         }
         //else
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
index 36d13de99dc0f5b351c8ebda9d2fd063186f4237..64a08aa86f29104ad725169adf6bce5887ceb105 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
@@ -483,43 +483,15 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 					//surfacetension
  
 
-					if ((phi2[DIR_000] <= phiLim) || phi[DIR_000] <= phiLim &&
-                        (
-						(phi[DIR_P00] > phiLim) ||
-						(phi[DIR_M00] > phiLim) ||
-						(phi[DIR_00P] > phiLim) ||
-						(phi[DIR_00M] > phiLim) ||
-						(phi[DIR_0M0] > phiLim) ||
-						(phi[DIR_0P0] > phiLim) ||
-						(phi[DIR_PP0] > phiLim) ||
-						(phi[DIR_PM0] > phiLim) ||
-						(phi[DIR_P0P] > phiLim) ||
-						(phi[DIR_P0M] > phiLim) ||
-						(phi[DIR_MP0] > phiLim) ||
-						(phi[DIR_MM0] > phiLim) ||
-						(phi[DIR_M0P] > phiLim) ||
-						(phi[DIR_M0M] > phiLim) ||
-						(phi[DIR_0PM] > phiLim) ||
-						(phi[DIR_0MM] > phiLim) ||
-						(phi[DIR_0PP] > phiLim) ||
-						(phi[DIR_0MP] > phiLim) ||
-						(phi[DIR_PPP] > phiLim) ||
-						(phi[DIR_PMP] > phiLim) ||
-						(phi[DIR_MPP] > phiLim) ||
-						(phi[DIR_MMP] > phiLim) ||
-						(phi[DIR_PPM] > phiLim) ||
-						(phi[DIR_PMM] > phiLim) ||
-						(phi[DIR_MPM] > phiLim) ||
-						(phi[DIR_MMM] > phiLim)
-						)) {
+					if (this->isGas(phiLim, phi, phi2)) {
 						real vx = (*vxNode)(x1, x2, x3);
 						real vy = (*vyNode)(x1, x2, x3);
 						real vz = (*vzNode)(x1, x2, x3);
                         //real rho = (*rhoNode)(x1, x2, x3);
 						findNeighbors(phaseField, x1, x2, x3);
-                        real dX1_phi = gradX1_phi();
-                        real dX2_phi = gradX2_phi();
-                        real dX3_phi = gradX3_phi();
+                        //real dX1_phi = gradX1_phi();
+                        //real dX2_phi = gradX2_phi();
+                        //real dX3_phi = gradX3_phi();
 						//real curv = computeCurvature_phi();
                         real laplacePressure = c12o1 * sigma * computeCurvature_phi();
 						findNeighbors(phaseFieldOld, x1, x2, x3);
@@ -551,25 +523,25 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 							rhoG = sumRho / sumWeight;// uncheck excpetion: what if there is no adequate neighbor?
 							for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
                                 if ((phi[fdir] > phiLim)) {
-									real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-									real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-									real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+									// real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+									// real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+									// real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                     //vx = vxBC;
                                     //vy = vyBC;
                                     //vz = vzBC;
-									real fPEQNeighbor = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir],c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+									//real fPEQNeighbor = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir],c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
 									//real fPEQNeighborInv = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                     //real vBC = (fPEQNeighborInv - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
 									//real fPEQHere = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , vx,vy,vz);
                                     //real vBC = (fPEQHere - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
-									real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
-									real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
+									//real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
+									//real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
 									//vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
                                    // real dvDir = vBC - vDir;
                                     // 27.04.23
-                                    real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                    real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                    real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                    //real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                    //real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                    //real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
                                     //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                    // real dvDir = (vBC - vIDir) * c1o2;
                                     //real dvDir = (vBC - vDir) ;
@@ -652,25 +624,25 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
                             if (phi2[DIR_000] <= phiLim) { // no refill liquid
 								for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
                                     if ((phi[fdir] > phiLim)) {
-										real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-										real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-										real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										// real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										// real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										// real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                         //vx = vxBC;
                                         //vy = vyBC;
                                         //vz = vzBC;
-										real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
+										//real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
 									//real fPEQNeighbor = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir],c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
 									//real fPEQNeighborInv = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                     //real vBC = (fPEQNeighborInv - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
                                     //									real fPEQHere = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , vx,vy,vz);
 //                                  real vBC = (fPEQHere - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
-										real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
+										//real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
 										//real dvDir = vBC - vDir;
                                         // real dvDir = vBC - vDir;
                                         // 27.04.23
-                                        real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                        real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                        real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        //real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        //real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        //real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
                                         //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                        // real dvDir = (vBC - vIDir) * c1o2;
                                        // real dvDir = (vBC - vDir) ;
@@ -844,22 +816,22 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 
 								for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
                                     if ((phi[fdir] > phiLim)) {
-										real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-										real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-										real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-										real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
+										// real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										// real vyBC = ((*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										// real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+										//real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
 									//real fPEQNeighbor = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir],c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
 									//real fPEQNeighborInv = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                     //real vBC = (fPEQNeighborInv - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
 									//real fPEQHere = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , vx,vy,vz);
                                     //real vBC = (fPEQHere - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
-										real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
+										//real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
 										//real dvDir = vBC - vDir;
 										//27.04.23
-                                        real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                        real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                        real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                        real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+                                        // real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        // real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        // real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                        //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                       //  real dvDir = (vBC - vIDir)*c1o2;
                                         //real dvDir = (vBC - vDir) ;
 
@@ -1041,16 +1013,16 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 											real vzBC = ((*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
 											real vBC = (D3Q27System::DX1[fdir] * vxBC + D3Q27System::DX2[fdir] * vyBC + D3Q27System::DX3[fdir] * vzBC);
 									//real fPEQNeighbor = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir],c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-									real fPEQNeighborInv = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+									//real fPEQNeighborInv = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
                                     //real vBC = (fPEQNeighborInv - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
 									//real fPEQHere = D3Q27System::getIncompFeqForDirection(fdir,c0o1 , vx,vy,vz);
                                     //real vBC = (fPEQHere - fPEQNeighbor) / WEIGTH[fdir] * c1o6;
 											real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
 											//real dvDir = vBC - vDir;
                                             // 27.04.23
-                                            real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                            real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-                                            real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                            //real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                            //real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+                                            //real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
                                             //real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                             //real dvDir = (vBC - vIDir) * c1o2;
                                             //real dvDir = (vBC - vDir) ;
@@ -3396,7 +3368,7 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
 					//collFactorM=(((*phaseField)(x1, x2, x3) > c1o2) && ((*phaseFieldOld)(x1, x2, x3) <= c1o2)) ? 1.8 : collFactorM;
 					real collFactorMInv = phi[DIR_000] > phiLim ? collFactorG : collFactorL;
 
-					real mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
+					//real mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
 
 					//----------- Calculating Macroscopic Values -------------
 					real rho = phi[DIR_000] > phiLim ? rhoH : rhoL;//rhoH + rhoToPhi * (phi[DIR_000] - phiH); //Incompressible
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.h b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.h
index 67a6680bb7d11980e357c0f797c185899618b324..85ada0fc290658fae8d3d92c1b2d8ae8dd0df471 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.h
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.h
@@ -115,6 +115,7 @@ protected:
     void computePhasefield();
     void findNeighbors(CbArray3D<real,IndexerX3X2X1>::CbArray3DPtr ph /*Phase-Field*/, int x1, int x2, int x3);
     void findNeighbors2(CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr ph, int x1, int x2, int x3);
+    bool isGas(real phiLim, real* phi, real* phi2);
 
     real nabla2_phi();
 
@@ -129,4 +130,43 @@ protected:
     real forcingX3;
 };
 
+/// @brief The function computes a fancy expression
+/// @param phiLim 
+/// @param phi 
+/// @param phi2 
+/// @return 
+inline bool MultiphaseScaleDistributionLBMKernel::isGas(real phiLim, real* phi, real* phi2)
+{
+    using namespace vf::lbm::dir;
+    return (phi2[DIR_000] <= phiLim) || ((phi[DIR_000] <= phiLim) &&
+                        (
+						(phi[DIR_P00] > phiLim) ||
+						(phi[DIR_M00] > phiLim) ||
+						(phi[DIR_00P] > phiLim) ||
+						(phi[DIR_00M] > phiLim) ||
+						(phi[DIR_0M0] > phiLim) ||
+						(phi[DIR_0P0] > phiLim) ||
+						(phi[DIR_PP0] > phiLim) ||
+						(phi[DIR_PM0] > phiLim) ||
+						(phi[DIR_P0P] > phiLim) ||
+						(phi[DIR_P0M] > phiLim) ||
+						(phi[DIR_MP0] > phiLim) ||
+						(phi[DIR_MM0] > phiLim) ||
+						(phi[DIR_M0P] > phiLim) ||
+						(phi[DIR_M0M] > phiLim) ||
+						(phi[DIR_0PM] > phiLim) ||
+						(phi[DIR_0MM] > phiLim) ||
+						(phi[DIR_0PP] > phiLim) ||
+						(phi[DIR_0MP] > phiLim) ||
+						(phi[DIR_PPP] > phiLim) ||
+						(phi[DIR_PMP] > phiLim) ||
+						(phi[DIR_MPP] > phiLim) ||
+						(phi[DIR_MMP] > phiLim) ||
+						(phi[DIR_PPM] > phiLim) ||
+						(phi[DIR_PMM] > phiLim) ||
+						(phi[DIR_MPM] > phiLim) ||
+						(phi[DIR_MMM] > phiLim)
+						));
+}
+
 #endif
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
index 915fddca7578a292e51cbd39bca35156117e5aee..5d32955a948f9c997832cfd9757a1887d0c58466 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseSharpInterfaceLBMKernel.cpp
@@ -361,9 +361,9 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
 									vBC = (vBC + vDir) / (c2o1 + vBC - vDir);	
 
 									//27.04.23
-									real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-									real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-									real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+									// real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+									// real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+									// real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
 									//real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
 									//real dvDir = (vBC - vIDir)*c1o2;
 									real dvDir = (vBC - vDir);
@@ -410,9 +410,9 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
 										real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
                                         //real dvDir = vBC - vDir;
 										//27.04.23
-										real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
 										//real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                         //real dvDir = (vBC - vIDir)*c1o2;
 										real dvDir = (vBC - vDir);
@@ -499,10 +499,10 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
 										real vDir = (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy + D3Q27System::DX3[fdir] * vz);
                                         //real dvDir = vBC - vDir;
 										//27.04.23
-										real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+										// real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										//real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
                                         //real dvDir = (vBC - vIDir)*c1o2;
 										real dvDir = (vBC - vDir);
 
@@ -609,10 +609,10 @@ void MultiphaseSharpInterfaceLBMKernel::calculate(int step)
 										vBC = (vBC + vDir) / (c2o1 + vBC - vDir);
 
 										//27.04.23
-										real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
-										real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
+										// real vxI = ((*vxNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vyI = ((*vyNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										// real vzI = ((*vzNode)(x1 - D3Q27System::DX1[fdir], x2 - D3Q27System::DX2[fdir], x3 - D3Q27System::DX3[fdir]));
+										//real vIDir = (D3Q27System::DX1[fdir] * vxI + D3Q27System::DX2[fdir] * vyI + D3Q27System::DX3[fdir] * vzI);
 										//real dvDir = (vBC - vIDir)*c1o2;
 										real dvDir = (vBC - vDir);