diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
index 6c1dd5df6301af0f1f3b62499fc3d007b858567b..7ed86bf8f4fff30cac707c567cde9f8526f4ae6a 100644
--- a/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.cpp
@@ -27,10 +27,12 @@ PorousMedia::PorousMedia(double porosity, unsigned int geoID, double darcySI, do
 	this->forchheimerSI = forchheimerSI;
 	this->dxLBM = dxLBM;
 	this->dtLBM = dtLBM;
-	this->forchheimerLBM = this->forchheimerSI * this->dxLBM;
-	this->darcyLBM = this->darcySI * this->dtLBM;
+	//this->lengthOfPorousMedia = 1.0;
+	//this->forchheimerLBM = this->forchheimerSI * this->dxLBM / this->lengthOfPorousMedia;
+	//this->darcyLBM = this->darcySI * this->dtLBM / this->lengthOfPorousMedia;
 	this->level = level;
 	setCoordinatesToZero();
+	setResistanceLBM();
 }
 
 PorousMedia::~PorousMedia()
@@ -72,6 +74,20 @@ void PorousMedia::setSizePM(unsigned int sizePM)
 	this->sizePM = sizePM;
 }
 
+void PorousMedia::setResistanceLBM() 
+{
+	if ((this->endCoordX - this->startCoordX) > 1.0)
+	{
+		this->lengthOfPorousMedia = (this->endCoordX - this->startCoordX);
+	}
+	else
+	{
+		this->lengthOfPorousMedia = 1.0;
+	}
+	this->forchheimerLBM = this->forchheimerSI * this->dxLBM / this->lengthOfPorousMedia;
+	this->darcyLBM = this->darcySI * this->dtLBM / this->lengthOfPorousMedia;
+}
+
 //void PorousMedia::definePMarea(Parameter* para, unsigned int level)
 //{
 //	unsigned int counter = 0;
diff --git a/src/VirtualFluids_GPU/Calculation/PorousMedia.h b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
index a6d49037d43f5872211b6da124fc3a8f0b60c7ec..961a873ce329b5e657c64db4e5e852de515385ee 100644
--- a/src/VirtualFluids_GPU/Calculation/PorousMedia.h
+++ b/src/VirtualFluids_GPU/Calculation/PorousMedia.h
@@ -27,6 +27,7 @@ public:
 	void setHostNodeIDsPM(unsigned int* hostNodeIDsPM);
 	void setDeviceNodeIDsPM(unsigned int* deviceNodeIDsPM);
 	void setSizePM(unsigned int sizePM);
+	void setResistanceLBM();
 
 	//getter
 	double getPorosity();
@@ -60,6 +61,7 @@ private:
 	double endCoordX;
 	double endCoordY;
 	double endCoordZ;
+	double lengthOfPorousMedia;
 	unsigned int geoID;
 	//std::vector< unsigned int > nodeIDsPorousMedia;
 	unsigned int sizePM;
diff --git a/src/VirtualFluids_GPU/GPU/CalcMac27.cu b/src/VirtualFluids_GPU/GPU/CalcMac27.cu
index d70cc9724f70afa4549589f7ed28953fb72a9a74..c6e46d7913f0d9a346596fb23ecfa6cecc625ee1 100644
--- a/src/VirtualFluids_GPU/GPU/CalcMac27.cu
+++ b/src/VirtualFluids_GPU/GPU/CalcMac27.cu
@@ -553,7 +553,7 @@ extern "C" __global__ void LBCalcMacCompSP27( doubflo* vxD,
 	  vyD[k]    = zero;
 	  vzD[k]    = zero;
 
-      if(geoD[k] == GEO_FLUID)
+      if(geoD[k] == GEO_FLUID || geoD[k] == GEO_PM_0)
       {
          rhoD[k]    =   (D.f[dirE   ])[ke  ]+ (D.f[dirW   ])[kw  ]+ 
                         (D.f[dirN   ])[kn  ]+ (D.f[dirS   ])[ks  ]+
diff --git a/src/VirtualFluids_GPU/GPU/Cumulant27.cu b/src/VirtualFluids_GPU/GPU/Cumulant27.cu
index 9eae9857791d926912242f84abf343eaa84e3c61..323572ab3a9dff4a937c7a8febe7800cfa204943 100644
--- a/src/VirtualFluids_GPU/GPU/Cumulant27.cu
+++ b/src/VirtualFluids_GPU/GPU/Cumulant27.cu
@@ -9773,7 +9773,7 @@ extern "C" __global__ void LB_Kernel_Kum_New_Comp_SP_27(doubflo omega,
 		unsigned int BC;
 		BC = bcMatD[k];
 
-		if( (BC != GEO_SOLID) && (BC != GEO_VOID) )
+		if( BC >= GEO_FLUID/*(BC != GEO_SOLID) && (BC != GEO_VOID)*/ )
 		{
 			Distributions27 D;
 			if (EvenOrOdd==true)
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
index e3cd75a7587b2d70ddbfbca77a653947d339f9bf..306a8ce7c286fb0a864d38b528ca7bfbeeb62239 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Interface.h
+++ b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
@@ -313,7 +313,6 @@ extern "C" void KernelWaleCumAA2016CompSP27( unsigned int numberOfThreads,
 
 extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads, 
 									   doubflo omega,
-									   unsigned int* bcMatD,
 									   unsigned int* neighborX,
 									   unsigned int* neighborY,
 									   unsigned int* neighborZ,
@@ -324,7 +323,8 @@ extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads,
 									   doubflo porosity,
 									   doubflo darcy,
 									   doubflo forchheimer,
-									   unsigned int porousMedia,
+									   unsigned int sizeOfPorousMedia,
+									   unsigned int* nodeIdsPorousMedia, 
 									   bool EvenOrOdd);
 
 extern "C" void KernelThS7(unsigned int numberOfThreads, 
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index 4eadea9766a24daa934fc5949d06fc6482cab778..f17a904c08724cdea0d9dcaf79507e82f7f60f2e 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -279,7 +279,6 @@ extern "C" __global__ void LB_Kernel_Wale_Cum_AA2016_Comp_SP_27( doubflo omega,
 																 bool EvenOrOdd);
 
 extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
-															unsigned int* bcMatD,
 															unsigned int* neighborX,
 															unsigned int* neighborY,
 															unsigned int* neighborZ,
@@ -290,7 +289,8 @@ extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
 															doubflo porosity,
 															doubflo darcy,
 															doubflo forchheimer,
-															unsigned int porousMedia,
+															unsigned int sizeOfPorousMedia,
+															unsigned int* nodeIdsPorousMedia,
 															bool EvenOrOdd);
 
 extern "C" __global__ void LB_Kernel_ThS7(doubflo diffusivity,
diff --git a/src/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
index 6c4d91e169b37d2705ea3e321c29ca1101bccc20..6c432d5d0ceb1ed8d6a7b14fe38b6a76e8d9f056 100644
--- a/src/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -959,7 +959,6 @@ extern "C" void KernelWaleCumOneCompSP27(unsigned int numberOfThreads,
 //////////////////////////////////////////////////////////////////////////
 extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads, 
 									   doubflo omega,
-									   unsigned int* bcMatD,
 									   unsigned int* neighborX,
 									   unsigned int* neighborY,
 									   unsigned int* neighborZ,
@@ -970,7 +969,8 @@ extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads,
 									   doubflo porosity,
 									   doubflo darcy,
 									   doubflo forchheimer,
-									   unsigned int porousMedia,
+									   unsigned int sizeOfPorousMedia,
+									   unsigned int* nodeIdsPorousMedia, 
 									   bool EvenOrOdd)
 {
 	int Grid = (size_Mat / numberOfThreads) + 1;
@@ -989,7 +989,6 @@ extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads,
 	dim3 threads(numberOfThreads, 1, 1);
 
 	LB_Kernel_PM_Cum_One_Comp_SP_27 <<< grid, threads >>>(omega,
-														  bcMatD,
 														  neighborX,
 														  neighborY,
 														  neighborZ,
@@ -1000,7 +999,8 @@ extern "C" void KernelPMCumOneCompSP27(unsigned int numberOfThreads,
 														  porosity,
 														  darcy,
 														  forchheimer,
-														  porousMedia,
+														  sizeOfPorousMedia,
+														  nodeIdsPorousMedia,
 														  EvenOrOdd); 
 	getLastCudaError("LB_Kernel_PM_Cum_One_Comp_SP_27 execution failed"); 
 }
diff --git a/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu b/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu
index fdbc535bbad930a4d2f403e8e725e08eeec42c36..124eab96c37b33b6bcdc33daa7ef035c95b27125 100644
--- a/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu
+++ b/src/VirtualFluids_GPU/GPU/PorousMediaCumulant27.cu
@@ -5,7 +5,6 @@
 
 ////////////////////////////////////////////////////////////////////////////////
 extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
-															unsigned int* bcMatD,
 															unsigned int* neighborX,
 															unsigned int* neighborY,
 															unsigned int* neighborZ,
@@ -16,9 +15,72 @@ extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
 															doubflo porosity,
 															doubflo darcy,
 															doubflo forchheimer,
-															unsigned int porousMedia,
+															unsigned int sizeOfPorousMedia,
+															unsigned int* nodeIdsPorousMedia,
 															bool EvenOrOdd)
 {
+	Distributions27 D;
+	if (EvenOrOdd==true)
+	{
+		D.f[dirE   ] = &DDStart[dirE   *size_Mat];
+		D.f[dirW   ] = &DDStart[dirW   *size_Mat];
+		D.f[dirN   ] = &DDStart[dirN   *size_Mat];
+		D.f[dirS   ] = &DDStart[dirS   *size_Mat];
+		D.f[dirT   ] = &DDStart[dirT   *size_Mat];
+		D.f[dirB   ] = &DDStart[dirB   *size_Mat];
+		D.f[dirNE  ] = &DDStart[dirNE  *size_Mat];
+		D.f[dirSW  ] = &DDStart[dirSW  *size_Mat];
+		D.f[dirSE  ] = &DDStart[dirSE  *size_Mat];
+		D.f[dirNW  ] = &DDStart[dirNW  *size_Mat];
+		D.f[dirTE  ] = &DDStart[dirTE  *size_Mat];
+		D.f[dirBW  ] = &DDStart[dirBW  *size_Mat];
+		D.f[dirBE  ] = &DDStart[dirBE  *size_Mat];
+		D.f[dirTW  ] = &DDStart[dirTW  *size_Mat];
+		D.f[dirTN  ] = &DDStart[dirTN  *size_Mat];
+		D.f[dirBS  ] = &DDStart[dirBS  *size_Mat];
+		D.f[dirBN  ] = &DDStart[dirBN  *size_Mat];
+		D.f[dirTS  ] = &DDStart[dirTS  *size_Mat];
+		D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+		D.f[dirTNE ] = &DDStart[dirTNE *size_Mat];
+		D.f[dirTSW ] = &DDStart[dirTSW *size_Mat];
+		D.f[dirTSE ] = &DDStart[dirTSE *size_Mat];
+		D.f[dirTNW ] = &DDStart[dirTNW *size_Mat];
+		D.f[dirBNE ] = &DDStart[dirBNE *size_Mat];
+		D.f[dirBSW ] = &DDStart[dirBSW *size_Mat];
+		D.f[dirBSE ] = &DDStart[dirBSE *size_Mat];
+		D.f[dirBNW ] = &DDStart[dirBNW *size_Mat];
+	}
+	else
+	{
+		D.f[dirW   ] = &DDStart[dirE   *size_Mat];
+		D.f[dirE   ] = &DDStart[dirW   *size_Mat];
+		D.f[dirS   ] = &DDStart[dirN   *size_Mat];
+		D.f[dirN   ] = &DDStart[dirS   *size_Mat];
+		D.f[dirB   ] = &DDStart[dirT   *size_Mat];
+		D.f[dirT   ] = &DDStart[dirB   *size_Mat];
+		D.f[dirSW  ] = &DDStart[dirNE  *size_Mat];
+		D.f[dirNE  ] = &DDStart[dirSW  *size_Mat];
+		D.f[dirNW  ] = &DDStart[dirSE  *size_Mat];
+		D.f[dirSE  ] = &DDStart[dirNW  *size_Mat];
+		D.f[dirBW  ] = &DDStart[dirTE  *size_Mat];
+		D.f[dirTE  ] = &DDStart[dirBW  *size_Mat];
+		D.f[dirTW  ] = &DDStart[dirBE  *size_Mat];
+		D.f[dirBE  ] = &DDStart[dirTW  *size_Mat];
+		D.f[dirBS  ] = &DDStart[dirTN  *size_Mat];
+		D.f[dirTN  ] = &DDStart[dirBS  *size_Mat];
+		D.f[dirTS  ] = &DDStart[dirBN  *size_Mat];
+		D.f[dirBN  ] = &DDStart[dirTS  *size_Mat];
+		D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
+		D.f[dirBSW ] = &DDStart[dirTNE *size_Mat];
+		D.f[dirBNE ] = &DDStart[dirTSW *size_Mat];
+		D.f[dirBNW ] = &DDStart[dirTSE *size_Mat];
+		D.f[dirBSE ] = &DDStart[dirTNW *size_Mat];
+		D.f[dirTSW ] = &DDStart[dirBNE *size_Mat];
+		D.f[dirTNE ] = &DDStart[dirBSW *size_Mat];
+		D.f[dirTNW ] = &DDStart[dirBSE *size_Mat];
+		D.f[dirTSE ] = &DDStart[dirBNW *size_Mat];
+	}
+
 	////////////////////////////////////////////////////////////////////////////////
 	const unsigned  x = threadIdx.x;  // Globaler x-Index 
 	const unsigned  y = blockIdx.x;   // Globaler y-Index 
@@ -27,915 +89,847 @@ extern "C" __global__ void LB_Kernel_PM_Cum_One_Comp_SP_27( doubflo omega,
 	const unsigned nx = blockDim.x;
 	const unsigned ny = gridDim.x;
 
-	const unsigned k = nx*(ny*z + y) + x;
+	const unsigned kt = nx*(ny*z + y) + x;
 	//////////////////////////////////////////////////////////////////////////
 
-	if(k<size_Mat)
+	if(kt<sizeOfPorousMedia)
 	{
 		////////////////////////////////////////////////////////////////////////////////
-		unsigned int BC;
-		BC = bcMatD[k];
-
-		if(BC == porousMedia)
+		//index
+		unsigned int k    = nodeIdsPorousMedia[kt];
+		unsigned int kw   = neighborX[k];
+		unsigned int ks   = neighborY[k];
+		unsigned int kb   = neighborZ[k];
+		unsigned int ksw  = neighborY[kw];
+		unsigned int kbw  = neighborZ[kw];
+		unsigned int kbs  = neighborZ[ks];
+		unsigned int kbsw = neighborZ[ksw];
+		//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		doubflo mfcbb = (D.f[dirE   ])[k  ];
+		doubflo mfabb = (D.f[dirW   ])[kw ];
+		doubflo mfbcb = (D.f[dirN   ])[k  ];
+		doubflo mfbab = (D.f[dirS   ])[ks ];
+		doubflo mfbbc = (D.f[dirT   ])[k  ];
+		doubflo mfbba = (D.f[dirB   ])[kb ];
+		doubflo mfccb = (D.f[dirNE  ])[k  ];
+		doubflo mfaab = (D.f[dirSW  ])[ksw];
+		doubflo mfcab = (D.f[dirSE  ])[ks ];
+		doubflo mfacb = (D.f[dirNW  ])[kw ];
+		doubflo mfcbc = (D.f[dirTE  ])[k  ];
+		doubflo mfaba = (D.f[dirBW  ])[kbw];
+		doubflo mfcba = (D.f[dirBE  ])[kb ];
+		doubflo mfabc = (D.f[dirTW  ])[kw ];
+		doubflo mfbcc = (D.f[dirTN  ])[k  ];
+		doubflo mfbaa = (D.f[dirBS  ])[kbs];
+		doubflo mfbca = (D.f[dirBN  ])[kb ];
+		doubflo mfbac = (D.f[dirTS  ])[ks ];
+		doubflo mfbbb = (D.f[dirZERO])[k  ];
+		doubflo mfccc = (D.f[dirTNE ])[k  ];
+		doubflo mfaac = (D.f[dirTSW ])[ksw];
+		doubflo mfcac = (D.f[dirTSE ])[ks ];
+		doubflo mfacc = (D.f[dirTNW ])[kw ];
+		doubflo mfcca = (D.f[dirBNE ])[kb ];
+		doubflo mfaaa = (D.f[dirBSW])[kbsw];
+		doubflo mfcaa = (D.f[dirBSE ])[kbs];
+		doubflo mfaca = (D.f[dirBNW ])[kbw];
+		////////////////////////////////////////////////////////////////////////////////////
+		doubflo drho = ((((mfccc+mfaaa) + (mfaca+mfcac)) + ((mfacc+mfcaa) + (mfaac+mfcca))) + 
+						(((mfbac+mfbca) + (mfbaa+mfbcc)) + ((mfabc+mfcba) + (mfaba+mfcbc)) + ((mfacb+mfcab) + (mfaab+mfccb))) +
+						((mfabb+mfcbb) + (mfbab+mfbcb)) + (mfbba+mfbbc)) + mfbbb;
+
+		doubflo rho = one+drho;
+		////////////////////////////////////////////////////////////////////////////////////
+		doubflo vvx    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfcaa-mfacc) + (mfcca-mfaac))) + 
+					     (((mfcba-mfabc) + (mfcbc-mfaba)) + ((mfcab-mfacb) + (mfccb-mfaab))) +
+					       (mfcbb-mfabb)) / rho;
+		doubflo vvy    =((((mfccc-mfaaa) + (mfaca-mfcac)) + ((mfacc-mfcaa) + (mfcca-mfaac))) + 
+			             (((mfbca-mfbac) + (mfbcc-mfbaa)) + ((mfacb-mfcab) + (mfccb-mfaab))) +
+			               (mfbcb-mfbab)) / rho;
+		doubflo vvz    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfacc-mfcaa) + (mfaac-mfcca))) + 
+			             (((mfbac-mfbca) + (mfbcc-mfbaa)) + ((mfabc-mfcba) + (mfcbc-mfaba))) +
+			               (mfbbc-mfbba)) / rho;
+		doubflo vx2    = vvx*vvx;
+		doubflo vy2    = vvy*vvy;
+		doubflo vz2    = vvz*vvz;
+		////////////////////////////////////////////////////////////////////////////////////
+		//porous media
+		vvx = -(two * vvx) / (-two - darcy - forchheimer * sqrtf(vx2 + vy2 + vz2));
+		vvy = -(two * vvy) / (-two - darcy - forchheimer * sqrtf(vx2 + vy2 + vz2));
+		vvz = -(two * vvz) / (-two - darcy - forchheimer * sqrtf(vx2 + vy2 + vz2));
+		//vvx = (two * vvx) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+		//vvy = (two * vvy) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+		//vvz = (two * vvz) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
+		////////////////////////////////////////////////////////////////////////////////////
+		//the force be with you
+		doubflo fx = forces[0]/(pow(two,level)); //zero;//0.0032653/(pow(two,level)); //0.000000005;//(two/1600000.0) / 120.0; //
+		doubflo fy = forces[1]/(pow(two,level)); //zero;
+		doubflo fz = forces[2]/(pow(two,level)); //zero;
+		vvx += fx*c1o2;
+		vvy += fy*c1o2;
+		vvz += fz*c1o2;
+		////////////////////////////////////////////////////////////////////////////////////
+		//doubflo omega = omega_in;
+		////////////////////////////////////////////////////////////////////////////////////
+		doubflo oMdrho = one; // comp special
+		////////////////////////////////////////////////////////////////////////////////////
+		doubflo m0, m1, m2;	
+		//////////////////////////////////////////////////////////////////////////////////////
+		doubflo wadjust;
+		doubflo qudricLimitP = 0.01f;// * 0.0001f;
+		doubflo qudricLimitM = 0.01f;// * 0.0001f;
+		doubflo qudricLimitD = 0.01f;// * 0.001f;
+		////////////////////////////////////////////////////////////////////////////////////
+		//Hin
+		////////////////////////////////////////////////////////////////////////////////////
+		// mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36  Konditionieren
+		////////////////////////////////////////////////////////////////////////////////////
+		// Z - Dir
+		m2    = mfaaa	+ mfaac;
+		m1    = mfaac	- mfaaa;
+		m0    = m2		+ mfaab;
+		mfaaa = m0;
+		m0   += c1o36 * oMdrho;	
+		mfaab = m1 -		m0 * vvz;
+		mfaac = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaba  + mfabc;
+		m1    = mfabc  - mfaba;
+		m0    = m2		+ mfabb;
+		mfaba = m0;
+		m0   += c1o9 * oMdrho;
+		mfabb = m1 -		m0 * vvz;
+		mfabc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaca  + mfacc;
+		m1    = mfacc  - mfaca;
+		m0    = m2		+ mfacb;
+		mfaca = m0;
+		m0   += c1o36 * oMdrho;
+		mfacb = m1 -		m0 * vvz;
+		mfacc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbaa	+ mfbac;
+		m1    = mfbac	- mfbaa;
+		m0    = m2		+ mfbab;
+		mfbaa = m0;
+		m0   += c1o9 * oMdrho;
+		mfbab = m1 -		m0 * vvz;
+		mfbac = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbba  + mfbbc;
+		m1    = mfbbc  - mfbba;
+		m0    = m2		+ mfbbb;
+		mfbba = m0;
+		m0   += c4o9 * oMdrho;
+		mfbbb = m1 -		m0 * vvz;
+		mfbbc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbca  + mfbcc;
+		m1    = mfbcc  - mfbca;
+		m0    = m2		+ mfbcb;
+		mfbca = m0;
+		m0   += c1o9 * oMdrho;
+		mfbcb = m1 -		m0 * vvz;
+		mfbcc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcaa	+ mfcac;
+		m1    = mfcac	- mfcaa;
+		m0    = m2		+ mfcab;
+		mfcaa = m0;
+		m0   += c1o36 * oMdrho;
+		mfcab = m1 -		m0 * vvz;
+		mfcac = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcba  + mfcbc;
+		m1    = mfcbc  - mfcba;
+		m0    = m2		+ mfcbb;
+		mfcba = m0;
+		m0   += c1o9 * oMdrho;
+		mfcbb = m1 -		m0 * vvz;
+		mfcbc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcca  + mfccc;
+		m1    = mfccc  - mfcca;
+		m0    = m2		+ mfccb;
+		mfcca = m0;
+		m0   += c1o36 * oMdrho;
+		mfccb = m1 -		m0 * vvz;
+		mfccc = m2 - two*	m1 * vvz + vz2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		// mit  1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
+		////////////////////////////////////////////////////////////////////////////////////
+		// Y - Dir
+		m2    = mfaaa	+ mfaca;
+		m1    = mfaca	- mfaaa;
+		m0    = m2		+ mfaba;
+		mfaaa = m0;
+		m0   += c1o6 * oMdrho;
+		mfaba = m1 -		m0 * vvy;
+		mfaca = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaab  + mfacb;
+		m1    = mfacb  - mfaab;
+		m0    = m2		+ mfabb;
+		mfaab = m0;
+		mfabb = m1 -		m0 * vvy;
+		mfacb = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaac  + mfacc;
+		m1    = mfacc  - mfaac;
+		m0    = m2		+ mfabc;
+		mfaac = m0;
+		m0   += c1o18 * oMdrho;
+		mfabc = m1 -		m0 * vvy;
+		mfacc = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbaa	+ mfbca;
+		m1    = mfbca	- mfbaa;
+		m0    = m2		+ mfbba;
+		mfbaa = m0;
+		m0   += c2o3 * oMdrho;
+		mfbba = m1 -		m0 * vvy;
+		mfbca = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbab  + mfbcb;
+		m1    = mfbcb  - mfbab;
+		m0    = m2		+ mfbbb;
+		mfbab = m0;
+		mfbbb = m1 -		m0 * vvy;
+		mfbcb = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfbac  + mfbcc;
+		m1    = mfbcc  - mfbac;
+		m0    = m2		+ mfbbc;
+		mfbac = m0;
+		m0   += c2o9 * oMdrho;
+		mfbbc = m1 -		m0 * vvy;
+		mfbcc = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcaa	+ mfcca;
+		m1    = mfcca	- mfcaa;
+		m0    = m2		+ mfcba;
+		mfcaa = m0;
+		m0   += c1o6 * oMdrho;
+		mfcba = m1 -		m0 * vvy;
+		mfcca = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcab  + mfccb;
+		m1    = mfccb  - mfcab;
+		m0    = m2		+ mfcbb;
+		mfcab = m0;
+		mfcbb = m1 -		m0 * vvy;
+		mfccb = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfcac  + mfccc;
+		m1    = mfccc  - mfcac;
+		m0    = m2		+ mfcbc;
+		mfcac = m0;
+		m0   += c1o18 * oMdrho;
+		mfcbc = m1 -		m0 * vvy;
+		mfccc = m2 - two*	m1 * vvy + vy2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		// mit     1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9		Konditionieren
+		////////////////////////////////////////////////////////////////////////////////////
+		// X - Dir
+		m2    = mfaaa	+ mfcaa;
+		m1    = mfcaa	- mfaaa;
+		m0    = m2		+ mfbaa;
+		mfaaa = m0;
+		m0   += one* oMdrho;
+		mfbaa = m1 -		m0 * vvx;
+		mfcaa = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaba  + mfcba;
+		m1    = mfcba  - mfaba;
+		m0    = m2		+ mfbba;
+		mfaba = m0;
+		mfbba = m1 -		m0 * vvx;
+		mfcba = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaca  + mfcca;
+		m1    = mfcca  - mfaca;
+		m0    = m2		+ mfbca;
+		mfaca = m0;
+		m0   += c1o3 * oMdrho;
+		mfbca = m1 -		m0 * vvx;
+		mfcca = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaab	+ mfcab;
+		m1    = mfcab	- mfaab;
+		m0    = m2		+ mfbab;
+		mfaab = m0;
+		mfbab = m1 -		m0 * vvx;
+		mfcab = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfabb  + mfcbb;
+		m1    = mfcbb  - mfabb;
+		m0    = m2		+ mfbbb;
+		mfabb = m0;
+		mfbbb = m1 -		m0 * vvx;
+		mfcbb = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfacb  + mfccb;
+		m1    = mfccb  - mfacb;
+		m0    = m2		+ mfbcb;
+		mfacb = m0;
+		mfbcb = m1 -		m0 * vvx;
+		mfccb = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfaac	+ mfcac;
+		m1    = mfcac	- mfaac;
+		m0    = m2		+ mfbac;
+		mfaac = m0;
+		m0   += c1o3 * oMdrho;
+		mfbac = m1 -		m0 * vvx;
+		mfcac = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfabc  + mfcbc;
+		m1    = mfcbc  - mfabc;
+		m0    = m2		+ mfbbc;
+		mfabc = m0;
+		mfbbc = m1 -		m0 * vvx;
+		mfcbc = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		m2    = mfacc  + mfccc;
+		m1    = mfccc  - mfacc;
+		m0    = m2		+ mfbcc;
+		mfacc = m0;
+		m0   += c1o9 * oMdrho;
+		mfbcc = m1 -		m0 * vvx;
+		mfccc = m2 - two*	m1 * vvx + vx2 * m0;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+
+
+		////////////////////////////////////////////////////////////////////////////////////
+		// Cumulants
+		////////////////////////////////////////////////////////////////////////////////////
+		doubflo OxxPyyPzz = one;	//set the bulk viscosity one is high / two is very low and zero is (too) high
+
+		////////////////////////////////////////////////////////////
+		//3.
+		//////////////////////////////
+		doubflo OxyyPxzz  = one;
+		doubflo OxyyMxzz  = one;
+		doubflo Oxyz      = one;
+		////////////////////////////////////////////////////////////
+		//4.
+		//////////////////////////////
+		doubflo O4        = one;
+		////////////////////////////////////////////////////////////
+		//5.
+		//////////////////////////////
+		doubflo O5        = one;
+		////////////////////////////////////////////////////////////
+		//6.
+		//////////////////////////////
+		doubflo O6        = one;
+		////////////////////////////////////////////////////////////
+
+
+		//central moments to cumulants
+		//4.
+		doubflo CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+		doubflo CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+		doubflo CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+			  	 		
+		doubflo CUMcca = mfcca - (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
+		doubflo CUMcac = mfcac - (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
+		doubflo CUMacc = mfacc - (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
+
+		//5.
+		doubflo CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
+		doubflo CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
+		doubflo CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
+		
+		//6.
+
+		doubflo CUMccc = mfccc + ((-four *  mfbbb * mfbbb  
+						-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+						-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+						-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+						+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+						+     two * (mfcaa * mfaca * mfaac)
+						+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
+						-    c1o3 * (mfacc + mfcac + mfcca) /rho 
+						-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
+						+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
+						+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
+						+ c1o27*((drho * drho - drho)/(rho*rho)));
+						//+ c1o27*(one -three/rho +two/(rho*rho)));
+
+
+
+
+		//2.
+		// linear combinations
+		doubflo mxxPyyPzz = mfcaa + mfaca + mfaac;
+		doubflo mxxMyy    = mfcaa - mfaca;
+		doubflo mxxMzz	  = mfcaa - mfaac;
+		
+		//////////////////////////////////////////////////////////////////////////
+		doubflo magicBulk=(CUMacc+CUMcac+CUMcca)*(one/OxxPyyPzz-c1o2)*c3o2*8.;
+
+		//////////////////////////////////////////////////////////////////////////
+		//limiter-Scheise Teil 1
+		//doubflo oxxyy,oxxzz,oxy,oxz,oyz;
+		//doubflo smag=0.001;
+		//oxxyy    = omega+(one-omega)*abs(mxxMyy)/(abs(mxxMyy)+smag);
+		//oxxzz    = omega+(one-omega)*abs(mxxMzz)/(abs(mxxMzz)+smag);
+		//oxy      = omega+(one-omega)*abs(mfbba)/(abs(mfbba)+smag);
+		//oxz      = omega+(one-omega)*abs(mfbab)/(abs(mfbab)+smag);
+		//oyz      = omega+(one-omega)*abs(mfabb)/(abs(mfabb)+smag);
+
+		////////////////////////////////////////////////////////////////////////////
+		////Teil 1b
+		//doubflo constante = 1000.0;
+		//doubflo nuEddi = constante * abs(mxxPyyPzz);
+		//doubflo omegaLimit = one / (one / omega + three * nuEddi);
+
+		//{
+		//	doubflo dxux = c1o2 * (-omegaLimit) *(mxxMyy + mxxMzz) +  OxxPyyPzz * (mfaaa - mxxPyyPzz);
+		//	doubflo dyuy = dxux + omegaLimit * c3o2 * mxxMyy;
+		//	doubflo dzuz = dxux + omegaLimit * c3o2 * mxxMzz;
+
+			////relax
+			//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+			//mxxMyy    += omegaLimit * (-mxxMyy) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vy2 * dyuy);
+			//mxxMzz    += omegaLimit * (-mxxMzz) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vz2 * dzuz);
+
+		//}
+		//mfabb     += omegaLimit * (-mfabb);
+		//mfbab     += omegaLimit * (-mfbab);
+		//mfbba     += omegaLimit * (-mfbba);
+		////////////////////////////////////////////////////////////////////////////
+
+		///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		//incl. correction		(hat noch nicht so gut funktioniert...Optimierungsbedarf??)
 		{
-			Distributions27 D;
-			if (EvenOrOdd==true)
-			{
-				D.f[dirE   ] = &DDStart[dirE   *size_Mat];
-				D.f[dirW   ] = &DDStart[dirW   *size_Mat];
-				D.f[dirN   ] = &DDStart[dirN   *size_Mat];
-				D.f[dirS   ] = &DDStart[dirS   *size_Mat];
-				D.f[dirT   ] = &DDStart[dirT   *size_Mat];
-				D.f[dirB   ] = &DDStart[dirB   *size_Mat];
-				D.f[dirNE  ] = &DDStart[dirNE  *size_Mat];
-				D.f[dirSW  ] = &DDStart[dirSW  *size_Mat];
-				D.f[dirSE  ] = &DDStart[dirSE  *size_Mat];
-				D.f[dirNW  ] = &DDStart[dirNW  *size_Mat];
-				D.f[dirTE  ] = &DDStart[dirTE  *size_Mat];
-				D.f[dirBW  ] = &DDStart[dirBW  *size_Mat];
-				D.f[dirBE  ] = &DDStart[dirBE  *size_Mat];
-				D.f[dirTW  ] = &DDStart[dirTW  *size_Mat];
-				D.f[dirTN  ] = &DDStart[dirTN  *size_Mat];
-				D.f[dirBS  ] = &DDStart[dirBS  *size_Mat];
-				D.f[dirBN  ] = &DDStart[dirBN  *size_Mat];
-				D.f[dirTS  ] = &DDStart[dirTS  *size_Mat];
-				D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
-				D.f[dirTNE ] = &DDStart[dirTNE *size_Mat];
-				D.f[dirTSW ] = &DDStart[dirTSW *size_Mat];
-				D.f[dirTSE ] = &DDStart[dirTSE *size_Mat];
-				D.f[dirTNW ] = &DDStart[dirTNW *size_Mat];
-				D.f[dirBNE ] = &DDStart[dirBNE *size_Mat];
-				D.f[dirBSW ] = &DDStart[dirBSW *size_Mat];
-				D.f[dirBSE ] = &DDStart[dirBSE *size_Mat];
-				D.f[dirBNW ] = &DDStart[dirBNW *size_Mat];
-			}
-			else
-			{
-				D.f[dirW   ] = &DDStart[dirE   *size_Mat];
-				D.f[dirE   ] = &DDStart[dirW   *size_Mat];
-				D.f[dirS   ] = &DDStart[dirN   *size_Mat];
-				D.f[dirN   ] = &DDStart[dirS   *size_Mat];
-				D.f[dirB   ] = &DDStart[dirT   *size_Mat];
-				D.f[dirT   ] = &DDStart[dirB   *size_Mat];
-				D.f[dirSW  ] = &DDStart[dirNE  *size_Mat];
-				D.f[dirNE  ] = &DDStart[dirSW  *size_Mat];
-				D.f[dirNW  ] = &DDStart[dirSE  *size_Mat];
-				D.f[dirSE  ] = &DDStart[dirNW  *size_Mat];
-				D.f[dirBW  ] = &DDStart[dirTE  *size_Mat];
-				D.f[dirTE  ] = &DDStart[dirBW  *size_Mat];
-				D.f[dirTW  ] = &DDStart[dirBE  *size_Mat];
-				D.f[dirBE  ] = &DDStart[dirTW  *size_Mat];
-				D.f[dirBS  ] = &DDStart[dirTN  *size_Mat];
-				D.f[dirTN  ] = &DDStart[dirBS  *size_Mat];
-				D.f[dirTS  ] = &DDStart[dirBN  *size_Mat];
-				D.f[dirBN  ] = &DDStart[dirTS  *size_Mat];
-				D.f[dirZERO] = &DDStart[dirZERO*size_Mat];
-				D.f[dirBSW ] = &DDStart[dirTNE *size_Mat];
-				D.f[dirBNE ] = &DDStart[dirTSW *size_Mat];
-				D.f[dirBNW ] = &DDStart[dirTSE *size_Mat];
-				D.f[dirBSE ] = &DDStart[dirTNW *size_Mat];
-				D.f[dirTSW ] = &DDStart[dirBNE *size_Mat];
-				D.f[dirTNE ] = &DDStart[dirBSW *size_Mat];
-				D.f[dirTNW ] = &DDStart[dirBSE *size_Mat];
-				D.f[dirTSE ] = &DDStart[dirBNW *size_Mat];
-			}
-
-			////////////////////////////////////////////////////////////////////////////////
-			//index
-			unsigned int kw   = neighborX[k];
-			unsigned int ks   = neighborY[k];
-			unsigned int kb   = neighborZ[k];
-			unsigned int ksw  = neighborY[kw];
-			unsigned int kbw  = neighborZ[kw];
-			unsigned int kbs  = neighborZ[ks];
-			unsigned int kbsw = neighborZ[ksw];
-			//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			doubflo mfcbb = (D.f[dirE   ])[k  ];
-			doubflo mfabb = (D.f[dirW   ])[kw ];
-			doubflo mfbcb = (D.f[dirN   ])[k  ];
-			doubflo mfbab = (D.f[dirS   ])[ks ];
-			doubflo mfbbc = (D.f[dirT   ])[k  ];
-			doubflo mfbba = (D.f[dirB   ])[kb ];
-			doubflo mfccb = (D.f[dirNE  ])[k  ];
-			doubflo mfaab = (D.f[dirSW  ])[ksw];
-			doubflo mfcab = (D.f[dirSE  ])[ks ];
-			doubflo mfacb = (D.f[dirNW  ])[kw ];
-			doubflo mfcbc = (D.f[dirTE  ])[k  ];
-			doubflo mfaba = (D.f[dirBW  ])[kbw];
-			doubflo mfcba = (D.f[dirBE  ])[kb ];
-			doubflo mfabc = (D.f[dirTW  ])[kw ];
-			doubflo mfbcc = (D.f[dirTN  ])[k  ];
-			doubflo mfbaa = (D.f[dirBS  ])[kbs];
-			doubflo mfbca = (D.f[dirBN  ])[kb ];
-			doubflo mfbac = (D.f[dirTS  ])[ks ];
-			doubflo mfbbb = (D.f[dirZERO])[k  ];
-			doubflo mfccc = (D.f[dirTNE ])[k  ];
-			doubflo mfaac = (D.f[dirTSW ])[ksw];
-			doubflo mfcac = (D.f[dirTSE ])[ks ];
-			doubflo mfacc = (D.f[dirTNW ])[kw ];
-			doubflo mfcca = (D.f[dirBNE ])[kb ];
-			doubflo mfaaa = (D.f[dirBSW])[kbsw];
-			doubflo mfcaa = (D.f[dirBSE ])[kbs];
-			doubflo mfaca = (D.f[dirBNW ])[kbw];
-			////////////////////////////////////////////////////////////////////////////////////
-			doubflo drho = ((((mfccc+mfaaa) + (mfaca+mfcac)) + ((mfacc+mfcaa) + (mfaac+mfcca))) + 
-							(((mfbac+mfbca) + (mfbaa+mfbcc)) + ((mfabc+mfcba) + (mfaba+mfcbc)) + ((mfacb+mfcab) + (mfaab+mfccb))) +
-							((mfabb+mfcbb) + (mfbab+mfbcb)) + (mfbba+mfbbc)) + mfbbb;
-
-			doubflo rho = one+drho;
-			////////////////////////////////////////////////////////////////////////////////////
-			doubflo vvx    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfcaa-mfacc) + (mfcca-mfaac))) + 
-						     (((mfcba-mfabc) + (mfcbc-mfaba)) + ((mfcab-mfacb) + (mfccb-mfaab))) +
-						       (mfcbb-mfabb)) / rho;
-			doubflo vvy    =((((mfccc-mfaaa) + (mfaca-mfcac)) + ((mfacc-mfcaa) + (mfcca-mfaac))) + 
-				             (((mfbca-mfbac) + (mfbcc-mfbaa)) + ((mfacb-mfcab) + (mfccb-mfaab))) +
-				               (mfbcb-mfbab)) / rho;
-			doubflo vvz    =((((mfccc-mfaaa) + (mfcac-mfaca)) + ((mfacc-mfcaa) + (mfaac-mfcca))) + 
-				             (((mfbac-mfbca) + (mfbcc-mfbaa)) + ((mfabc-mfcba) + (mfcbc-mfaba))) +
-				               (mfbbc-mfbba)) / rho;
-			doubflo vx2    = vvx*vvx;
-			doubflo vy2    = vvy*vvy;
-			doubflo vz2    = vvz*vvz;
-			////////////////////////////////////////////////////////////////////////////////////
-			//porous media
-			vvx = -(two * vvx) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
-			vvy = -(two * vvy) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
-			vvz = -(two * vvz) / (-two + darcy + forchheimer * sqrtf(vx2 + vy2 + vz2));
-			//vvx = (two * vvx) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
-			//vvy = (two * vvy) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
-			//vvz = (two * vvz) / (two + 134.4 + 0.0068287 * sqrtf(vx2 + vy2 + vz2));
-			////////////////////////////////////////////////////////////////////////////////////
-			//the force be with you
-			doubflo fx = forces[0]/(pow(two,level)); //zero;//0.0032653/(pow(two,level)); //0.000000005;//(two/1600000.0) / 120.0; //
-			doubflo fy = forces[1]/(pow(two,level)); //zero;
-			doubflo fz = forces[2]/(pow(two,level)); //zero;
-			vvx += fx*c1o2;
-			vvy += fy*c1o2;
-			vvz += fz*c1o2;
-			////////////////////////////////////////////////////////////////////////////////////
-			//doubflo omega = omega_in;
-			////////////////////////////////////////////////////////////////////////////////////
-			doubflo oMdrho = one; // comp special
-			////////////////////////////////////////////////////////////////////////////////////
-			doubflo m0, m1, m2;	
-			//////////////////////////////////////////////////////////////////////////////////////
-			doubflo wadjust;
-			doubflo qudricLimitP = 0.01f;// * 0.0001f;
-			doubflo qudricLimitM = 0.01f;// * 0.0001f;
-			doubflo qudricLimitD = 0.01f;// * 0.001f;
-			////////////////////////////////////////////////////////////////////////////////////
-			//Hin
-			////////////////////////////////////////////////////////////////////////////////////
-			// mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36  Konditionieren
-			////////////////////////////////////////////////////////////////////////////////////
-			// Z - Dir
-			m2    = mfaaa	+ mfaac;
-			m1    = mfaac	- mfaaa;
-			m0    = m2		+ mfaab;
-			mfaaa = m0;
-			m0   += c1o36 * oMdrho;	
-			mfaab = m1 -		m0 * vvz;
-			mfaac = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaba  + mfabc;
-			m1    = mfabc  - mfaba;
-			m0    = m2		+ mfabb;
-			mfaba = m0;
-			m0   += c1o9 * oMdrho;
-			mfabb = m1 -		m0 * vvz;
-			mfabc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaca  + mfacc;
-			m1    = mfacc  - mfaca;
-			m0    = m2		+ mfacb;
-			mfaca = m0;
-			m0   += c1o36 * oMdrho;
-			mfacb = m1 -		m0 * vvz;
-			mfacc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbaa	+ mfbac;
-			m1    = mfbac	- mfbaa;
-			m0    = m2		+ mfbab;
-			mfbaa = m0;
-			m0   += c1o9 * oMdrho;
-			mfbab = m1 -		m0 * vvz;
-			mfbac = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbba  + mfbbc;
-			m1    = mfbbc  - mfbba;
-			m0    = m2		+ mfbbb;
-			mfbba = m0;
-			m0   += c4o9 * oMdrho;
-			mfbbb = m1 -		m0 * vvz;
-			mfbbc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbca  + mfbcc;
-			m1    = mfbcc  - mfbca;
-			m0    = m2		+ mfbcb;
-			mfbca = m0;
-			m0   += c1o9 * oMdrho;
-			mfbcb = m1 -		m0 * vvz;
-			mfbcc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcaa	+ mfcac;
-			m1    = mfcac	- mfcaa;
-			m0    = m2		+ mfcab;
-			mfcaa = m0;
-			m0   += c1o36 * oMdrho;
-			mfcab = m1 -		m0 * vvz;
-			mfcac = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcba  + mfcbc;
-			m1    = mfcbc  - mfcba;
-			m0    = m2		+ mfcbb;
-			mfcba = m0;
-			m0   += c1o9 * oMdrho;
-			mfcbb = m1 -		m0 * vvz;
-			mfcbc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcca  + mfccc;
-			m1    = mfccc  - mfcca;
-			m0    = m2		+ mfccb;
-			mfcca = m0;
-			m0   += c1o36 * oMdrho;
-			mfccb = m1 -		m0 * vvz;
-			mfccc = m2 - two*	m1 * vvz + vz2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			// mit  1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
-			////////////////////////////////////////////////////////////////////////////////////
-			// Y - Dir
-			m2    = mfaaa	+ mfaca;
-			m1    = mfaca	- mfaaa;
-			m0    = m2		+ mfaba;
-			mfaaa = m0;
-			m0   += c1o6 * oMdrho;
-			mfaba = m1 -		m0 * vvy;
-			mfaca = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaab  + mfacb;
-			m1    = mfacb  - mfaab;
-			m0    = m2		+ mfabb;
-			mfaab = m0;
-			mfabb = m1 -		m0 * vvy;
-			mfacb = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaac  + mfacc;
-			m1    = mfacc  - mfaac;
-			m0    = m2		+ mfabc;
-			mfaac = m0;
-			m0   += c1o18 * oMdrho;
-			mfabc = m1 -		m0 * vvy;
-			mfacc = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbaa	+ mfbca;
-			m1    = mfbca	- mfbaa;
-			m0    = m2		+ mfbba;
-			mfbaa = m0;
-			m0   += c2o3 * oMdrho;
-			mfbba = m1 -		m0 * vvy;
-			mfbca = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbab  + mfbcb;
-			m1    = mfbcb  - mfbab;
-			m0    = m2		+ mfbbb;
-			mfbab = m0;
-			mfbbb = m1 -		m0 * vvy;
-			mfbcb = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfbac  + mfbcc;
-			m1    = mfbcc  - mfbac;
-			m0    = m2		+ mfbbc;
-			mfbac = m0;
-			m0   += c2o9 * oMdrho;
-			mfbbc = m1 -		m0 * vvy;
-			mfbcc = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcaa	+ mfcca;
-			m1    = mfcca	- mfcaa;
-			m0    = m2		+ mfcba;
-			mfcaa = m0;
-			m0   += c1o6 * oMdrho;
-			mfcba = m1 -		m0 * vvy;
-			mfcca = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcab  + mfccb;
-			m1    = mfccb  - mfcab;
-			m0    = m2		+ mfcbb;
-			mfcab = m0;
-			mfcbb = m1 -		m0 * vvy;
-			mfccb = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfcac  + mfccc;
-			m1    = mfccc  - mfcac;
-			m0    = m2		+ mfcbc;
-			mfcac = m0;
-			m0   += c1o18 * oMdrho;
-			mfcbc = m1 -		m0 * vvy;
-			mfccc = m2 - two*	m1 * vvy + vy2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			// mit     1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9		Konditionieren
-			////////////////////////////////////////////////////////////////////////////////////
-			// X - Dir
-			m2    = mfaaa	+ mfcaa;
-			m1    = mfcaa	- mfaaa;
-			m0    = m2		+ mfbaa;
-			mfaaa = m0;
-			m0   += one* oMdrho;
-			mfbaa = m1 -		m0 * vvx;
-			mfcaa = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaba  + mfcba;
-			m1    = mfcba  - mfaba;
-			m0    = m2		+ mfbba;
-			mfaba = m0;
-			mfbba = m1 -		m0 * vvx;
-			mfcba = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaca  + mfcca;
-			m1    = mfcca  - mfaca;
-			m0    = m2		+ mfbca;
-			mfaca = m0;
-			m0   += c1o3 * oMdrho;
-			mfbca = m1 -		m0 * vvx;
-			mfcca = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaab	+ mfcab;
-			m1    = mfcab	- mfaab;
-			m0    = m2		+ mfbab;
-			mfaab = m0;
-			mfbab = m1 -		m0 * vvx;
-			mfcab = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfabb  + mfcbb;
-			m1    = mfcbb  - mfabb;
-			m0    = m2		+ mfbbb;
-			mfabb = m0;
-			mfbbb = m1 -		m0 * vvx;
-			mfcbb = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfacb  + mfccb;
-			m1    = mfccb  - mfacb;
-			m0    = m2		+ mfbcb;
-			mfacb = m0;
-			mfbcb = m1 -		m0 * vvx;
-			mfccb = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfaac	+ mfcac;
-			m1    = mfcac	- mfaac;
-			m0    = m2		+ mfbac;
-			mfaac = m0;
-			m0   += c1o3 * oMdrho;
-			mfbac = m1 -		m0 * vvx;
-			mfcac = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfabc  + mfcbc;
-			m1    = mfcbc  - mfabc;
-			m0    = m2		+ mfbbc;
-			mfabc = m0;
-			mfbbc = m1 -		m0 * vvx;
-			mfcbc = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			m2    = mfacc  + mfccc;
-			m1    = mfccc  - mfacc;
-			m0    = m2		+ mfbcc;
-			mfacc = m0;
-			m0   += c1o9 * oMdrho;
-			mfbcc = m1 -		m0 * vvx;
-			mfccc = m2 - two*	m1 * vvx + vx2 * m0;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-
-
-			////////////////////////////////////////////////////////////////////////////////////
-			// Cumulants
-			////////////////////////////////////////////////////////////////////////////////////
-			doubflo OxxPyyPzz = one;	//set the bulk viscosity one is high / two is very low and zero is (too) high
-
-			////////////////////////////////////////////////////////////
-			//3.
-			//////////////////////////////
-			doubflo OxyyPxzz  = one;
-			doubflo OxyyMxzz  = one;
-			doubflo Oxyz      = one;
-			////////////////////////////////////////////////////////////
-			//4.
-			//////////////////////////////
-			doubflo O4        = one;
-			////////////////////////////////////////////////////////////
-			//5.
-			//////////////////////////////
-			doubflo O5        = one;
-			////////////////////////////////////////////////////////////
-			//6.
-			//////////////////////////////
-			doubflo O6        = one;
-			////////////////////////////////////////////////////////////
-
-
-			//central moments to cumulants
-			//4.
-			doubflo CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
-			doubflo CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
-			doubflo CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
-				  	 		
-			doubflo CUMcca = mfcca - (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
-			doubflo CUMcac = mfcac - (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
-			doubflo CUMacc = mfacc - (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
-
-			//5.
-			doubflo CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
-			doubflo CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
-			doubflo CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
-			
-			//6.
-
-			doubflo CUMccc = mfccc + ((-four *  mfbbb * mfbbb  
-							-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
-							-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
-							-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
-							+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
-							+     two * (mfcaa * mfaca * mfaac)
-							+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
-							-    c1o3 * (mfacc + mfcac + mfcca) /rho 
-							-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
-							+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
-							+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
-							+ c1o27*((drho * drho - drho)/(rho*rho)));
-							//+ c1o27*(one -three/rho +two/(rho*rho)));
-
-
-
-
-			//2.
-			// linear combinations
-			doubflo mxxPyyPzz = mfcaa + mfaca + mfaac;
-			doubflo mxxMyy    = mfcaa - mfaca;
-			doubflo mxxMzz	  = mfcaa - mfaac;
-			
-			//////////////////////////////////////////////////////////////////////////
-// 			doubflo magicBulk=(CUMacc+CUMcac+CUMcca)*(one/OxxPyyPzz-c1o2)*c3o2*8.;
-
-			//////////////////////////////////////////////////////////////////////////
-			//limiter-Scheise Teil 1
-			//doubflo oxxyy,oxxzz,oxy,oxz,oyz;
-			//doubflo smag=0.001;
-			//oxxyy    = omega+(one-omega)*abs(mxxMyy)/(abs(mxxMyy)+smag);
-			//oxxzz    = omega+(one-omega)*abs(mxxMzz)/(abs(mxxMzz)+smag);
-			//oxy      = omega+(one-omega)*abs(mfbba)/(abs(mfbba)+smag);
-			//oxz      = omega+(one-omega)*abs(mfbab)/(abs(mfbab)+smag);
-			//oyz      = omega+(one-omega)*abs(mfabb)/(abs(mfabb)+smag);
-
-			////////////////////////////////////////////////////////////////////////////
-			////Teil 1b
-			//doubflo constante = 1000.0;
-			//doubflo nuEddi = constante * abs(mxxPyyPzz);
-			//doubflo omegaLimit = one / (one / omega + three * nuEddi);
-
-			//{
-			//	doubflo dxux = c1o2 * (-omegaLimit) *(mxxMyy + mxxMzz) +  OxxPyyPzz * (mfaaa - mxxPyyPzz);
-			//	doubflo dyuy = dxux + omegaLimit * c3o2 * mxxMyy;
-			//	doubflo dzuz = dxux + omegaLimit * c3o2 * mxxMzz;
-
-				////relax
-				//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
-				//mxxMyy    += omegaLimit * (-mxxMyy) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vy2 * dyuy);
-				//mxxMzz    += omegaLimit * (-mxxMzz) - three * (one + c1o2 * (-omegaLimit)) * (vx2 * dxux + vz2 * dzuz);
-
-			//}
-			//mfabb     += omegaLimit * (-mfabb);
-			//mfbab     += omegaLimit * (-mfbab);
-			//mfbba     += omegaLimit * (-mfbba);
-			////////////////////////////////////////////////////////////////////////////
-
- 			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- 			//incl. correction		(hat noch nicht so gut funktioniert...Optimierungsbedarf??)
- 			{
- 				doubflo dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 *  OxxPyyPzz * (mfaaa - mxxPyyPzz);
- 				doubflo dyuy = dxux + omega * c3o2 * mxxMyy;
- 				doubflo dzuz = dxux + omega * c3o2 * mxxMzz;
- 
- 				////relax original
- 				//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
- 				//mxxMyy    += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
- 				//mxxMzz    += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
-				//relax porous media Mp
-				mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz + rho * ( vx2 + vy2 + vz2) * (one / porosity - one)) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
-			    //relax porous media Ms
-				mxxMyy += omega * (rho * (vx2 - vy2) * (one / porosity - one) - mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
-				mxxMzz += omega * (rho * (vx2 - vz2) * (one / porosity - one) - mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
-
- 				//////////////////////////////////////////////////////////////////////////
- 				//limiter-Scheise Teil 2
- 				//mxxMyy    += oxxyy * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vy2 * dyuy);
- 				//mxxMzz    += oxxzz * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vz2 * dzuz);
- 				//////////////////////////////////////////////////////////////////////////
- 
- 			}
- 			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- 			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- 			////no correction
- 			//mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz);//-magicBulk*OxxPyyPzz;
- 			//mxxMyy    += -(-omega) * (-mxxMyy);
- 			//mxxMzz    += -(-omega) * (-mxxMzz);
- 			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			////original
-			//mfabb     += omega * (-mfabb);
-			//mfbab     += omega * (-mfbab);
-			//mfbba     += omega * (-mfbba);
-			/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-			//porous media M11
-			mfabb += omega * (rho * vvy * vvz * (one / porosity - one) - mfabb);
-			mfbab += omega * (rho * vvx * vvz * (one / porosity - one) - mfbab);
-			mfbba += omega * (rho * vvx * vvy * (one / porosity - one) - mfbba);
+			doubflo dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 *  OxxPyyPzz * (mfaaa - mxxPyyPzz);
+			doubflo dyuy = dxux + omega * c3o2 * mxxMyy;
+			doubflo dzuz = dxux + omega * c3o2 * mxxMzz;
+
+			////relax original
+			//mxxPyyPzz += OxxPyyPzz*(mfaaa  - mxxPyyPzz)- three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
+			//mxxMyy    += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+			//mxxMzz    += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+			//relax porous media Mp
+			mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz + rho * ( vx2 + vy2 + vz2) * (one / porosity - one)) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);//-magicBulk*OxxPyyPzz;
+		    //relax porous media Ms
+			mxxMyy += omega * (rho * (vx2 - vy2) * (one / porosity - one) - mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+			mxxMzz += omega * (rho * (vx2 - vz2) * (one / porosity - one) - mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
 
 			//////////////////////////////////////////////////////////////////////////
-			//limiter-Scheise Teil 3
-			//mfabb     += oyz * (-mfabb);
-			//mfbab     += oxz * (-mfbab);
-			//mfbba     += oxy * (-mfbba);
+			//limiter-Scheise Teil 2
+			//mxxMyy    += oxxyy * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vy2 * dyuy);
+			//mxxMzz    += oxxzz * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux + vz2 * dzuz);
 			//////////////////////////////////////////////////////////////////////////
 
-			// linear combinations back
-			mfcaa = c1o3 * (       mxxMyy +      mxxMzz + mxxPyyPzz);
-			mfaca = c1o3 * (-two*  mxxMyy +      mxxMzz + mxxPyyPzz);
-			mfaac = c1o3 * (       mxxMyy - two* mxxMzz + mxxPyyPzz);
-
-			//3.
-			// linear combinations
-
-			doubflo mxxyPyzz = mfcba + mfabc;
-			doubflo mxxyMyzz = mfcba - mfabc;
-
-			doubflo mxxzPyyz = mfcab + mfacb;
-			doubflo mxxzMyyz = mfcab - mfacb;
-
-			doubflo mxyyPxzz = mfbca + mfbac;
-			doubflo mxyyMxzz = mfbca - mfbac;
-
-			//relax
-			//////////////////////////////////////////////////////////////////////////
-			//das ist der limiter
- 			wadjust    = Oxyz+(one-Oxyz)*abs(mfbbb)/(abs(mfbbb)+qudricLimitD);
- 			mfbbb     += wadjust * (-mfbbb);
- 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxyPyzz)/(abs(mxxyPyzz)+qudricLimitP);
- 			mxxyPyzz  += wadjust * (-mxxyPyzz);
- 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxyMyzz)/(abs(mxxyMyzz)+qudricLimitM);
- 			mxxyMyzz  += wadjust * (-mxxyMyzz);
- 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxzPyyz)/(abs(mxxzPyyz)+qudricLimitP);
- 			mxxzPyyz  += wadjust * (-mxxzPyyz);
- 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxzMyyz)/(abs(mxxzMyyz)+qudricLimitM);
- 			mxxzMyyz  += wadjust * (-mxxzMyyz);
- 			wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxyyPxzz)/(abs(mxyyPxzz)+qudricLimitP);
- 			mxyyPxzz  += wadjust * (-mxyyPxzz);
- 			wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxyyMxzz)/(abs(mxyyMxzz)+qudricLimitM);
- 			mxyyMxzz  += wadjust * (-mxyyMxzz);
-			//////////////////////////////////////////////////////////////////////////
-			//ohne limiter
-			//mfbbb     += OxyyMxzz * (-mfbbb);
-			//mxxyPyzz  += OxyyPxzz * (-mxxyPyzz);
-			//mxxyMyzz  += OxyyMxzz * (-mxxyMyzz);
-			//mxxzPyyz  += OxyyPxzz * (-mxxzPyyz);
-			//mxxzMyyz  += OxyyMxzz * (-mxxzMyyz);
-			//mxyyPxzz  += OxyyPxzz * (-mxyyPxzz);
-			//mxyyMxzz  += OxyyMxzz * (-mxyyMxzz);
-			//////////////////////////////////////////////////////////////////////////
-
-			//// 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.
-			//////////////////////////////////////////////////////////////////////////
-			//mit limiter
- 		//	wadjust    = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit);
-			//CUMacc    += wadjust * (-CUMacc);
- 		//	wadjust    = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit);
-			//CUMcac    += wadjust * (-CUMcac); 
- 		//	wadjust    = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit);
-			//CUMcca    += wadjust * (-CUMcca); 
-
- 		//	wadjust    = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit);
-			//CUMbbc    += wadjust * (-CUMbbc); 
- 		//	wadjust    = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit);
-			//CUMbcb    += wadjust * (-CUMbcb); 
- 		//	wadjust    = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit);
-			//CUMcbb    += wadjust * (-CUMcbb); 
-			//////////////////////////////////////////////////////////////////////////
-			//ohne limiter
-			CUMacc += O4 * (-CUMacc); 
-			CUMcac += O4 * (-CUMcac); 
-			CUMcca += O4 * (-CUMcca); 
-
-			CUMbbc += O4 * (-CUMbbc); 
-			CUMbcb += O4 * (-CUMbcb); 
-			CUMcbb += O4 * (-CUMcbb); 
-			//////////////////////////////////////////////////////////////////////////
-			
-					
-			//5.
-			CUMbcc += O5 * (-CUMbcc);
-			CUMcbc += O5 * (-CUMcbc);
-			CUMccb += O5 * (-CUMccb);
-
-			//6.
-			CUMccc += O6 * (-CUMccc);
-			
-
-
-			//back cumulants to central moments
-			//4.
-			mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
-			mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
-			mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
-						   
-			mfcca = CUMcca + (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
-			mfcac = CUMcac + (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
-			mfacc = CUMacc + (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
-
-			//5.
-			mfbcc = CUMbcc + ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
-			mfcbc = CUMcbc + ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
-			mfccb = CUMccb + ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
-			
-			//6.
-
-			mfccc = CUMccc - ((-four *  mfbbb * mfbbb  
-							-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
-							-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
-							-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
-							+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
-							+     two * (mfcaa * mfaca * mfaac)
-							+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
-							-    c1o3 * (mfacc + mfcac + mfcca) /rho 
-							-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
-							+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
-							+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
-							+ c1o27*((drho * drho - drho)/(rho*rho)));
-							//+ c1o27*(one -three/rho +two/(rho*rho)));
-			////////////////////////////////////////////////////////////////////////////////////
-
-			////////////////////////////////////////////////////////////////////////////////////
-			//the force be with you
-			mfbaa = -mfbaa;
-			mfaba = -mfaba;
-			mfaab = -mfaab;
-			////////////////////////////////////////////////////////////////////////////////////
-
-
-			////////////////////////////////////////////////////////////////////////////////////
-			//back
-			////////////////////////////////////////////////////////////////////////////////////
-			//mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
-			////////////////////////////////////////////////////////////////////////////////////
-			// Z - Dir
-			m0 =  mfaac * c1o2 +      mfaab * (vvz - c1o2) + (mfaaa + one* oMdrho) * (     vz2 - vvz) * c1o2; 
-			m1 = -mfaac        - two* mfaab *  vvz         +  mfaaa                * (one- vz2)              - one* oMdrho * vz2; 
-			m2 =  mfaac * c1o2 +      mfaab * (vvz + c1o2) + (mfaaa + one* oMdrho) * (     vz2 + vvz) * c1o2;
-			mfaaa = m0;
-			mfaab = m1;
-			mfaac = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfabc * c1o2 +      mfabb * (vvz - c1o2) + mfaba * (     vz2 - vvz) * c1o2; 
-			m1 = -mfabc        - two* mfabb *  vvz         + mfaba * (one- vz2); 
-			m2 =  mfabc * c1o2 +      mfabb * (vvz + c1o2) + mfaba * (     vz2 + vvz) * c1o2;
-			mfaba = m0;
-			mfabb = m1;
-			mfabc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfacc * c1o2 +      mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
-			m1 = -mfacc        - two* mfacb *  vvz         +  mfaca                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
-			m2 =  mfacc * c1o2 +      mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
-			mfaca = m0;
-			mfacb = m1;
-			mfacc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfbac * c1o2 +      mfbab * (vvz - c1o2) + mfbaa * (     vz2 - vvz) * c1o2; 
-			m1 = -mfbac        - two* mfbab *  vvz         + mfbaa * (one- vz2); 
-			m2 =  mfbac * c1o2 +      mfbab * (vvz + c1o2) + mfbaa * (     vz2 + vvz) * c1o2;
-			mfbaa = m0;
-			mfbab = m1;
-			mfbac = m2;
-			/////////b//////////////////////////////////////////////////////////////////////////
-			m0 =  mfbbc * c1o2 +      mfbbb * (vvz - c1o2) + mfbba * (     vz2 - vvz) * c1o2; 
-			m1 = -mfbbc        - two* mfbbb *  vvz         + mfbba * (one- vz2); 
-			m2 =  mfbbc * c1o2 +      mfbbb * (vvz + c1o2) + mfbba * (     vz2 + vvz) * c1o2;
-			mfbba = m0;
-			mfbbb = m1;
-			mfbbc = m2;
-			/////////b//////////////////////////////////////////////////////////////////////////
-			m0 =  mfbcc * c1o2 +      mfbcb * (vvz - c1o2) + mfbca * (     vz2 - vvz) * c1o2; 
-			m1 = -mfbcc        - two* mfbcb *  vvz         + mfbca * (one- vz2); 
-			m2 =  mfbcc * c1o2 +      mfbcb * (vvz + c1o2) + mfbca * (     vz2 + vvz) * c1o2;
-			mfbca = m0;
-			mfbcb = m1;
-			mfbcc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcac * c1o2 +      mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
-			m1 = -mfcac        - two* mfcab *  vvz         +  mfcaa                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
-			m2 =  mfcac * c1o2 +      mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
-			mfcaa = m0;
-			mfcab = m1;
-			mfcac = m2;
-			/////////c//////////////////////////////////////////////////////////////////////////
-			m0 =  mfcbc * c1o2 +      mfcbb * (vvz - c1o2) + mfcba * (     vz2 - vvz) * c1o2; 
-			m1 = -mfcbc        - two* mfcbb *  vvz         + mfcba * (one- vz2); 
-			m2 =  mfcbc * c1o2 +      mfcbb * (vvz + c1o2) + mfcba * (     vz2 + vvz) * c1o2;
-			mfcba = m0;
-			mfcbb = m1;
-			mfcbc = m2;
-			/////////c//////////////////////////////////////////////////////////////////////////
-			m0 =  mfccc * c1o2 +      mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 - vvz) * c1o2; 
-			m1 = -mfccc        - two* mfccb *  vvz         +  mfcca                  * (one- vz2)              - c1o9 * oMdrho * vz2; 
-			m2 =  mfccc * c1o2 +      mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 + vvz) * 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 * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfaca        - two* mfaba *  vvy         +  mfaaa                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
-			m2 =  mfaca * c1o2 +      mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
-			mfaaa = m0;
-			mfaba = m1;
-			mfaca = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfacb * c1o2 +      mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfacb        - two* mfabb *  vvy         +  mfaab                  * (one- vy2)              - c2o3 * oMdrho * vy2; 
-			m2 =  mfacb * c1o2 +      mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 + vvy) * c1o2;
-			mfaab = m0;
-			mfabb = m1;
-			mfacb = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfacc * c1o2 +      mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfacc        - two* mfabc *  vvy         +  mfaac                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
-			m2 =  mfacc * c1o2 +      mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
-			mfaac = m0;
-			mfabc = m1;
-			mfacc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfbca * c1o2 +      mfbba * (vvy - c1o2) + mfbaa * (     vy2 - vvy) * c1o2; 
-			m1 = -mfbca        - two* mfbba *  vvy         + mfbaa * (one- vy2); 
-			m2 =  mfbca * c1o2 +      mfbba * (vvy + c1o2) + mfbaa * (     vy2 + vvy) * c1o2;
-			mfbaa = m0;
-			mfbba = m1;
-			mfbca = m2;
-			/////////b//////////////////////////////////////////////////////////////////////////
-			m0 =  mfbcb * c1o2 +      mfbbb * (vvy - c1o2) + mfbab * (     vy2 - vvy) * c1o2; 
-			m1 = -mfbcb        - two* mfbbb *  vvy         + mfbab * (one- vy2); 
-			m2 =  mfbcb * c1o2 +      mfbbb * (vvy + c1o2) + mfbab * (     vy2 + vvy) * c1o2;
-			mfbab = m0;
-			mfbbb = m1;
-			mfbcb = m2;
-			/////////b//////////////////////////////////////////////////////////////////////////
-			m0 =  mfbcc * c1o2 +      mfbbc * (vvy - c1o2) + mfbac * (     vy2 - vvy) * c1o2; 
-			m1 = -mfbcc        - two* mfbbc *  vvy         + mfbac * (one- vy2); 
-			m2 =  mfbcc * c1o2 +      mfbbc * (vvy + c1o2) + mfbac * (     vy2 + vvy) * c1o2;
-			mfbac = m0;
-			mfbbc = m1;
-			mfbcc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcca * c1o2 +      mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfcca        - two* mfcba *  vvy         +  mfcaa                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
-			m2 =  mfcca * c1o2 +      mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 + vvy) * c1o2;
-			mfcaa = m0;
-			mfcba = m1;
-			mfcca = m2;
-			/////////c//////////////////////////////////////////////////////////////////////////
-			m0 =  mfccb * c1o2 +      mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfccb        - two* mfcbb *  vvy         +  mfcab                  * (one- vy2)              - c2o9 * oMdrho * vy2; 
-			m2 =  mfccb * c1o2 +      mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 + vvy) * c1o2;
-			mfcab = m0;
-			mfcbb = m1;
-			mfccb = m2;
-			/////////c//////////////////////////////////////////////////////////////////////////
-			m0 =  mfccc * c1o2 +      mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
-			m1 = -mfccc        - two* mfcbc *  vvy         +  mfcac                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
-			m2 =  mfccc * c1o2 +      mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 + vvy) * 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 * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcaa        - two* mfbaa *  vvx         +  mfaaa                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
-			m2 =  mfcaa * c1o2 +      mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfaaa = m0;
-			mfbaa = m1;
-			mfcaa = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcba * c1o2 +      mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcba        - two* mfbba *  vvx         +  mfaba                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
-			m2 =  mfcba * c1o2 +      mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfaba = m0;
-			mfbba = m1;
-			mfcba = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcca * c1o2 +      mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcca        - two* mfbca *  vvx         +  mfaca                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
-			m2 =  mfcca * c1o2 +      mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfaca = m0;
-			mfbca = m1;
-			mfcca = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcab * c1o2 +      mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcab        - two* mfbab *  vvx         +  mfaab                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
-			m2 =  mfcab * c1o2 +      mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfaab = m0;
-			mfbab = m1;
-			mfcab = m2;
-			///////////b////////////////////////////////////////////////////////////////////////
-			m0 =  mfcbb * c1o2 +      mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcbb        - two* mfbbb *  vvx         +  mfabb                  * (one- vx2)              - c4o9 * oMdrho * vx2; 
-			m2 =  mfcbb * c1o2 +      mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfabb = m0;
-			mfbbb = m1;
-			mfcbb = m2;
-			///////////b////////////////////////////////////////////////////////////////////////
-			m0 =  mfccb * c1o2 +      mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfccb        - two* mfbcb *  vvx         +  mfacb                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
-			m2 =  mfccb * c1o2 +      mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfacb = m0;
-			mfbcb = m1;
-			mfccb = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-			////////////////////////////////////////////////////////////////////////////////////
-			m0 =  mfcac * c1o2 +      mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcac        - two* mfbac *  vvx         +  mfaac                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
-			m2 =  mfcac * c1o2 +      mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfaac = m0;
-			mfbac = m1;
-			mfcac = m2;
-			///////////c////////////////////////////////////////////////////////////////////////
-			m0 =  mfcbc * c1o2 +      mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfcbc        - two* mfbbc *  vvx         +  mfabc                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
-			m2 =  mfcbc * c1o2 +      mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfabc = m0;
-			mfbbc = m1;
-			mfcbc = m2;
-			///////////c////////////////////////////////////////////////////////////////////////
-			m0 =  mfccc * c1o2 +      mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
-			m1 = -mfccc        - two* mfbcc *  vvx         +  mfacc                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
-			m2 =  mfccc * c1o2 +      mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
-			mfacc = m0;
-			mfbcc = m1;
-			mfccc = m2;
-			////////////////////////////////////////////////////////////////////////////////////
-
-			////////////////////////////////////////////////////////////////////////////////////
-			(D.f[ dirE   ])[k   ] = mfabb;                                                                     
-			(D.f[ dirW   ])[kw  ] = mfcbb;                                                                   
-			(D.f[ dirN   ])[k   ] = mfbab;
-			(D.f[ dirS   ])[ks  ] = mfbcb;
-			(D.f[ dirT   ])[k   ] = mfbba;
-			(D.f[ dirB   ])[kb  ] = mfbbc;
-			(D.f[ dirNE  ])[k   ] = mfaab;
-			(D.f[ dirSW  ])[ksw ] = mfccb;
-			(D.f[ dirSE  ])[ks  ] = mfacb;
-			(D.f[ dirNW  ])[kw  ] = mfcab;
-			(D.f[ dirTE  ])[k   ] = mfaba;
-			(D.f[ dirBW  ])[kbw ] = mfcbc;
-			(D.f[ dirBE  ])[kb  ] = mfabc;
-			(D.f[ dirTW  ])[kw  ] = mfcba;
-			(D.f[ dirTN  ])[k   ] = mfbaa;
-			(D.f[ dirBS  ])[kbs ] = mfbcc;
-			(D.f[ dirBN  ])[kb  ] = mfbac;
-			(D.f[ dirTS  ])[ks  ] = mfbca;
-			(D.f[ dirZERO])[k   ] = mfbbb;
-			(D.f[ dirTNE ])[k   ] = mfaaa;
-			(D.f[ dirTSE ])[ks  ] = mfaca;
-			(D.f[ dirBNE ])[kb  ] = mfaac;
-			(D.f[ dirBSE ])[kbs ] = mfacc;
-			(D.f[ dirTNW ])[kw  ] = mfcaa;
-			(D.f[ dirTSW ])[ksw ] = mfcca;
-			(D.f[ dirBNW ])[kbw ] = mfcac;
-			(D.f[ dirBSW ])[kbsw] = mfccc;
-			////////////////////////////////////////////////////////////////////////////////////
-		}                                                                                                                    
+		}
+		///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+		/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		////no correction
+		//mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz);//-magicBulk*OxxPyyPzz;
+		//mxxMyy    += -(-omega) * (-mxxMyy);
+		//mxxMzz    += -(-omega) * (-mxxMzz);
+		/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		////original
+		//mfabb     += omega * (-mfabb);
+		//mfbab     += omega * (-mfbab);
+		//mfbba     += omega * (-mfbba);
+		/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		//porous media M11
+		mfabb += omega * (rho * vvy * vvz * (one / porosity - one) - mfabb);
+		mfbab += omega * (rho * vvx * vvz * (one / porosity - one) - mfbab);
+		mfbba += omega * (rho * vvx * vvy * (one / porosity - one) - mfbba);
+
+		//////////////////////////////////////////////////////////////////////////
+		//limiter-Scheise Teil 3
+		//mfabb     += oyz * (-mfabb);
+		//mfbab     += oxz * (-mfbab);
+		//mfbba     += oxy * (-mfbba);
+		//////////////////////////////////////////////////////////////////////////
+
+		// linear combinations back
+		mfcaa = c1o3 * (       mxxMyy +      mxxMzz + mxxPyyPzz);
+		mfaca = c1o3 * (-two*  mxxMyy +      mxxMzz + mxxPyyPzz);
+		mfaac = c1o3 * (       mxxMyy - two* mxxMzz + mxxPyyPzz);
+
+		//3.
+		// linear combinations
+
+		doubflo mxxyPyzz = mfcba + mfabc;
+		doubflo mxxyMyzz = mfcba - mfabc;
+
+		doubflo mxxzPyyz = mfcab + mfacb;
+		doubflo mxxzMyyz = mfcab - mfacb;
+
+		doubflo mxyyPxzz = mfbca + mfbac;
+		doubflo mxyyMxzz = mfbca - mfbac;
+
+		//relax
+		//////////////////////////////////////////////////////////////////////////
+		//das ist der limiter
+		wadjust    = Oxyz+(one-Oxyz)*abs(mfbbb)/(abs(mfbbb)+qudricLimitD);
+		mfbbb     += wadjust * (-mfbbb);
+		wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxyPyzz)/(abs(mxxyPyzz)+qudricLimitP);
+		mxxyPyzz  += wadjust * (-mxxyPyzz);
+		wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxyMyzz)/(abs(mxxyMyzz)+qudricLimitM);
+		mxxyMyzz  += wadjust * (-mxxyMyzz);
+		wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxxzPyyz)/(abs(mxxzPyyz)+qudricLimitP);
+		mxxzPyyz  += wadjust * (-mxxzPyyz);
+		wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxxzMyyz)/(abs(mxxzMyyz)+qudricLimitM);
+		mxxzMyyz  += wadjust * (-mxxzMyyz);
+		wadjust    = OxyyPxzz+(one-OxyyPxzz)*abs(mxyyPxzz)/(abs(mxyyPxzz)+qudricLimitP);
+		mxyyPxzz  += wadjust * (-mxyyPxzz);
+		wadjust    = OxyyMxzz+(one-OxyyMxzz)*abs(mxyyMxzz)/(abs(mxyyMxzz)+qudricLimitM);
+		mxyyMxzz  += wadjust * (-mxyyMxzz);
+		//////////////////////////////////////////////////////////////////////////
+		//ohne limiter
+		//mfbbb     += OxyyMxzz * (-mfbbb);
+		//mxxyPyzz  += OxyyPxzz * (-mxxyPyzz);
+		//mxxyMyzz  += OxyyMxzz * (-mxxyMyzz);
+		//mxxzPyyz  += OxyyPxzz * (-mxxzPyyz);
+		//mxxzMyyz  += OxyyMxzz * (-mxxzMyyz);
+		//mxyyPxzz  += OxyyPxzz * (-mxyyPxzz);
+		//mxyyMxzz  += OxyyMxzz * (-mxyyMxzz);
+		//////////////////////////////////////////////////////////////////////////
+
+		//// 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.
+		//////////////////////////////////////////////////////////////////////////
+		//mit limiter
+	//	wadjust    = O4+(one-O4)*abs(CUMacc)/(abs(CUMacc)+qudricLimit);
+		//CUMacc    += wadjust * (-CUMacc);
+	//	wadjust    = O4+(one-O4)*abs(CUMcac)/(abs(CUMcac)+qudricLimit);
+		//CUMcac    += wadjust * (-CUMcac); 
+	//	wadjust    = O4+(one-O4)*abs(CUMcca)/(abs(CUMcca)+qudricLimit);
+		//CUMcca    += wadjust * (-CUMcca); 
+
+	//	wadjust    = O4+(one-O4)*abs(CUMbbc)/(abs(CUMbbc)+qudricLimit);
+		//CUMbbc    += wadjust * (-CUMbbc); 
+	//	wadjust    = O4+(one-O4)*abs(CUMbcb)/(abs(CUMbcb)+qudricLimit);
+		//CUMbcb    += wadjust * (-CUMbcb); 
+	//	wadjust    = O4+(one-O4)*abs(CUMcbb)/(abs(CUMcbb)+qudricLimit);
+		//CUMcbb    += wadjust * (-CUMcbb); 
+		//////////////////////////////////////////////////////////////////////////
+		//ohne limiter
+		CUMacc += O4 * (-CUMacc); 
+		CUMcac += O4 * (-CUMcac); 
+		CUMcca += O4 * (-CUMcca); 
+
+		CUMbbc += O4 * (-CUMbbc); 
+		CUMbcb += O4 * (-CUMbcb); 
+		CUMcbb += O4 * (-CUMcbb); 
+		//////////////////////////////////////////////////////////////////////////
+		
+				
+		//5.
+		CUMbcc += O5 * (-CUMbcc);
+		CUMcbc += O5 * (-CUMcbc);
+		CUMccb += O5 * (-CUMccb);
+
+		//6.
+		CUMccc += O6 * (-CUMccc);
+		
+
+
+		//back cumulants to central moments
+		//4.
+		mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + two * mfbba * mfbab) / rho;
+		mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + two * mfbba * mfabb) / rho;
+		mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + two * mfbab * mfabb) / rho;
+					   
+		mfcca = CUMcca + (((mfcaa * mfaca + two * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) / rho  - c1o9*(drho/rho));
+		mfcac = CUMcac + (((mfcaa * mfaac + two * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) / rho  - c1o9*(drho/rho));
+		mfacc = CUMacc + (((mfaac * mfaca + two * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) / rho  - c1o9*(drho/rho));
+
+		//5.
+		mfbcc = CUMbcc + ((mfaac * mfbca + mfaca * mfbac + four * mfabb * mfbbb + two * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) ) / rho ;
+		mfcbc = CUMcbc + ((mfaac * mfcba + mfcaa * mfabc + four * mfbab * mfbbb + two * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) ) / rho ;
+		mfccb = CUMccb + ((mfcaa * mfacb + mfaca * mfcab + four * mfbba * mfbbb + two * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) ) / rho ;
+		
+		//6.
+
+		mfccc = CUMccc - ((-four *  mfbbb * mfbbb  
+						-           (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+						-    four * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+						-     two * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) / rho
+						+(   four * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+						+     two * (mfcaa * mfaca * mfaac)
+						+ sixteen *  mfbba * mfbab * mfabb) / (rho * rho)
+						-    c1o3 * (mfacc + mfcac + mfcca) /rho 
+						-    c1o9 * (mfcaa + mfaca + mfaac) /rho 
+						+(    two * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) 
+						+           (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) / (rho * rho) * c2o3 
+						+ c1o27*((drho * drho - drho)/(rho*rho)));
+						//+ c1o27*(one -three/rho +two/(rho*rho)));
+		////////////////////////////////////////////////////////////////////////////////////
+
+		////////////////////////////////////////////////////////////////////////////////////
+		//the force be with you
+		mfbaa = -mfbaa;
+		mfaba = -mfaba;
+		mfaab = -mfaab;
+		////////////////////////////////////////////////////////////////////////////////////
+
+
+		////////////////////////////////////////////////////////////////////////////////////
+		//back
+		////////////////////////////////////////////////////////////////////////////////////
+		//mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9   Konditionieren
+		////////////////////////////////////////////////////////////////////////////////////
+		// Z - Dir
+		m0 =  mfaac * c1o2 +      mfaab * (vvz - c1o2) + (mfaaa + one* oMdrho) * (     vz2 - vvz) * c1o2; 
+		m1 = -mfaac        - two* mfaab *  vvz         +  mfaaa                * (one- vz2)              - one* oMdrho * vz2; 
+		m2 =  mfaac * c1o2 +      mfaab * (vvz + c1o2) + (mfaaa + one* oMdrho) * (     vz2 + vvz) * c1o2;
+		mfaaa = m0;
+		mfaab = m1;
+		mfaac = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfabc * c1o2 +      mfabb * (vvz - c1o2) + mfaba * (     vz2 - vvz) * c1o2; 
+		m1 = -mfabc        - two* mfabb *  vvz         + mfaba * (one- vz2); 
+		m2 =  mfabc * c1o2 +      mfabb * (vvz + c1o2) + mfaba * (     vz2 + vvz) * c1o2;
+		mfaba = m0;
+		mfabb = m1;
+		mfabc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfacc * c1o2 +      mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
+		m1 = -mfacc        - two* mfacb *  vvz         +  mfaca                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
+		m2 =  mfacc * c1o2 +      mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
+		mfaca = m0;
+		mfacb = m1;
+		mfacc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfbac * c1o2 +      mfbab * (vvz - c1o2) + mfbaa * (     vz2 - vvz) * c1o2; 
+		m1 = -mfbac        - two* mfbab *  vvz         + mfbaa * (one- vz2); 
+		m2 =  mfbac * c1o2 +      mfbab * (vvz + c1o2) + mfbaa * (     vz2 + vvz) * c1o2;
+		mfbaa = m0;
+		mfbab = m1;
+		mfbac = m2;
+		/////////b//////////////////////////////////////////////////////////////////////////
+		m0 =  mfbbc * c1o2 +      mfbbb * (vvz - c1o2) + mfbba * (     vz2 - vvz) * c1o2; 
+		m1 = -mfbbc        - two* mfbbb *  vvz         + mfbba * (one- vz2); 
+		m2 =  mfbbc * c1o2 +      mfbbb * (vvz + c1o2) + mfbba * (     vz2 + vvz) * c1o2;
+		mfbba = m0;
+		mfbbb = m1;
+		mfbbc = m2;
+		/////////b//////////////////////////////////////////////////////////////////////////
+		m0 =  mfbcc * c1o2 +      mfbcb * (vvz - c1o2) + mfbca * (     vz2 - vvz) * c1o2; 
+		m1 = -mfbcc        - two* mfbcb *  vvz         + mfbca * (one- vz2); 
+		m2 =  mfbcc * c1o2 +      mfbcb * (vvz + c1o2) + mfbca * (     vz2 + vvz) * c1o2;
+		mfbca = m0;
+		mfbcb = m1;
+		mfbcc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcac * c1o2 +      mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 - vvz) * c1o2; 
+		m1 = -mfcac        - two* mfcab *  vvz         +  mfcaa                  * (one- vz2)              - c1o3 * oMdrho * vz2; 
+		m2 =  mfcac * c1o2 +      mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (     vz2 + vvz) * c1o2;
+		mfcaa = m0;
+		mfcab = m1;
+		mfcac = m2;
+		/////////c//////////////////////////////////////////////////////////////////////////
+		m0 =  mfcbc * c1o2 +      mfcbb * (vvz - c1o2) + mfcba * (     vz2 - vvz) * c1o2; 
+		m1 = -mfcbc        - two* mfcbb *  vvz         + mfcba * (one- vz2); 
+		m2 =  mfcbc * c1o2 +      mfcbb * (vvz + c1o2) + mfcba * (     vz2 + vvz) * c1o2;
+		mfcba = m0;
+		mfcbb = m1;
+		mfcbc = m2;
+		/////////c//////////////////////////////////////////////////////////////////////////
+		m0 =  mfccc * c1o2 +      mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 - vvz) * c1o2; 
+		m1 = -mfccc        - two* mfccb *  vvz         +  mfcca                  * (one- vz2)              - c1o9 * oMdrho * vz2; 
+		m2 =  mfccc * c1o2 +      mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (     vz2 + vvz) * 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 * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfaca        - two* mfaba *  vvy         +  mfaaa                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
+		m2 =  mfaca * c1o2 +      mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
+		mfaaa = m0;
+		mfaba = m1;
+		mfaca = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfacb * c1o2 +      mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfacb        - two* mfabb *  vvy         +  mfaab                  * (one- vy2)              - c2o3 * oMdrho * vy2; 
+		m2 =  mfacb * c1o2 +      mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (     vy2 + vvy) * c1o2;
+		mfaab = m0;
+		mfabb = m1;
+		mfacb = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfacc * c1o2 +      mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfacc        - two* mfabc *  vvy         +  mfaac                  * (one- vy2)              - c1o6 * oMdrho * vy2; 
+		m2 =  mfacc * c1o2 +      mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (     vy2 + vvy) * c1o2;
+		mfaac = m0;
+		mfabc = m1;
+		mfacc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfbca * c1o2 +      mfbba * (vvy - c1o2) + mfbaa * (     vy2 - vvy) * c1o2; 
+		m1 = -mfbca        - two* mfbba *  vvy         + mfbaa * (one- vy2); 
+		m2 =  mfbca * c1o2 +      mfbba * (vvy + c1o2) + mfbaa * (     vy2 + vvy) * c1o2;
+		mfbaa = m0;
+		mfbba = m1;
+		mfbca = m2;
+		/////////b//////////////////////////////////////////////////////////////////////////
+		m0 =  mfbcb * c1o2 +      mfbbb * (vvy - c1o2) + mfbab * (     vy2 - vvy) * c1o2; 
+		m1 = -mfbcb        - two* mfbbb *  vvy         + mfbab * (one- vy2); 
+		m2 =  mfbcb * c1o2 +      mfbbb * (vvy + c1o2) + mfbab * (     vy2 + vvy) * c1o2;
+		mfbab = m0;
+		mfbbb = m1;
+		mfbcb = m2;
+		/////////b//////////////////////////////////////////////////////////////////////////
+		m0 =  mfbcc * c1o2 +      mfbbc * (vvy - c1o2) + mfbac * (     vy2 - vvy) * c1o2; 
+		m1 = -mfbcc        - two* mfbbc *  vvy         + mfbac * (one- vy2); 
+		m2 =  mfbcc * c1o2 +      mfbbc * (vvy + c1o2) + mfbac * (     vy2 + vvy) * c1o2;
+		mfbac = m0;
+		mfbbc = m1;
+		mfbcc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcca * c1o2 +      mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfcca        - two* mfcba *  vvy         +  mfcaa                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
+		m2 =  mfcca * c1o2 +      mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (     vy2 + vvy) * c1o2;
+		mfcaa = m0;
+		mfcba = m1;
+		mfcca = m2;
+		/////////c//////////////////////////////////////////////////////////////////////////
+		m0 =  mfccb * c1o2 +      mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfccb        - two* mfcbb *  vvy         +  mfcab                  * (one- vy2)              - c2o9 * oMdrho * vy2; 
+		m2 =  mfccb * c1o2 +      mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (     vy2 + vvy) * c1o2;
+		mfcab = m0;
+		mfcbb = m1;
+		mfccb = m2;
+		/////////c//////////////////////////////////////////////////////////////////////////
+		m0 =  mfccc * c1o2 +      mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 - vvy) * c1o2; 
+		m1 = -mfccc        - two* mfcbc *  vvy         +  mfcac                   * (one- vy2)              - c1o18 * oMdrho * vy2; 
+		m2 =  mfccc * c1o2 +      mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (     vy2 + vvy) * 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 * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcaa        - two* mfbaa *  vvx         +  mfaaa                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+		m2 =  mfcaa * c1o2 +      mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfaaa = m0;
+		mfbaa = m1;
+		mfcaa = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcba * c1o2 +      mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcba        - two* mfbba *  vvx         +  mfaba                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+		m2 =  mfcba * c1o2 +      mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfaba = m0;
+		mfbba = m1;
+		mfcba = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcca * c1o2 +      mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcca        - two* mfbca *  vvx         +  mfaca                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+		m2 =  mfcca * c1o2 +      mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfaca = m0;
+		mfbca = m1;
+		mfcca = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcab * c1o2 +      mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcab        - two* mfbab *  vvx         +  mfaab                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+		m2 =  mfcab * c1o2 +      mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfaab = m0;
+		mfbab = m1;
+		mfcab = m2;
+		///////////b////////////////////////////////////////////////////////////////////////
+		m0 =  mfcbb * c1o2 +      mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcbb        - two* mfbbb *  vvx         +  mfabb                  * (one- vx2)              - c4o9 * oMdrho * vx2; 
+		m2 =  mfcbb * c1o2 +      mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfabb = m0;
+		mfbbb = m1;
+		mfcbb = m2;
+		///////////b////////////////////////////////////////////////////////////////////////
+		m0 =  mfccb * c1o2 +      mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfccb        - two* mfbcb *  vvx         +  mfacb                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+		m2 =  mfccb * c1o2 +      mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfacb = m0;
+		mfbcb = m1;
+		mfccb = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+		////////////////////////////////////////////////////////////////////////////////////
+		m0 =  mfcac * c1o2 +      mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcac        - two* mfbac *  vvx         +  mfaac                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+		m2 =  mfcac * c1o2 +      mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfaac = m0;
+		mfbac = m1;
+		mfcac = m2;
+		///////////c////////////////////////////////////////////////////////////////////////
+		m0 =  mfcbc * c1o2 +      mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfcbc        - two* mfbbc *  vvx         +  mfabc                  * (one- vx2)              - c1o9 * oMdrho * vx2; 
+		m2 =  mfcbc * c1o2 +      mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfabc = m0;
+		mfbbc = m1;
+		mfcbc = m2;
+		///////////c////////////////////////////////////////////////////////////////////////
+		m0 =  mfccc * c1o2 +      mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 - vvx) * c1o2; 
+		m1 = -mfccc        - two* mfbcc *  vvx         +  mfacc                   * (one- vx2)              - c1o36 * oMdrho * vx2; 
+		m2 =  mfccc * c1o2 +      mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (     vx2 + vvx) * c1o2;
+		mfacc = m0;
+		mfbcc = m1;
+		mfccc = m2;
+		////////////////////////////////////////////////////////////////////////////////////
+
+		////////////////////////////////////////////////////////////////////////////////////
+		(D.f[ dirE   ])[k   ] = mfabb;                                                                     
+		(D.f[ dirW   ])[kw  ] = mfcbb;                                                                   
+		(D.f[ dirN   ])[k   ] = mfbab;
+		(D.f[ dirS   ])[ks  ] = mfbcb;
+		(D.f[ dirT   ])[k   ] = mfbba;
+		(D.f[ dirB   ])[kb  ] = mfbbc;
+		(D.f[ dirNE  ])[k   ] = mfaab;
+		(D.f[ dirSW  ])[ksw ] = mfccb;
+		(D.f[ dirSE  ])[ks  ] = mfacb;
+		(D.f[ dirNW  ])[kw  ] = mfcab;
+		(D.f[ dirTE  ])[k   ] = mfaba;
+		(D.f[ dirBW  ])[kbw ] = mfcbc;
+		(D.f[ dirBE  ])[kb  ] = mfabc;
+		(D.f[ dirTW  ])[kw  ] = mfcba;
+		(D.f[ dirTN  ])[k   ] = mfbaa;
+		(D.f[ dirBS  ])[kbs ] = mfbcc;
+		(D.f[ dirBN  ])[kb  ] = mfbac;
+		(D.f[ dirTS  ])[ks  ] = mfbca;
+		(D.f[ dirZERO])[k   ] = mfbbb;
+		(D.f[ dirTNE ])[k   ] = mfaaa;
+		(D.f[ dirTSE ])[ks  ] = mfaca;
+		(D.f[ dirBNE ])[kb  ] = mfaac;
+		(D.f[ dirBSE ])[kbs ] = mfacc;
+		(D.f[ dirTNW ])[kw  ] = mfcaa;
+		(D.f[ dirTSW ])[ksw ] = mfcca;
+		(D.f[ dirBNW ])[kbw ] = mfcac;
+		(D.f[ dirBSW ])[kbsw] = mfccc;
+		////////////////////////////////////////////////////////////////////////////////////
 	}
 }
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.cpp b/src/VirtualFluids_GPU/LBM/Simulation.cpp
index debcfc4367c5a1049b0f8f209a3a179af87b4188..4dd6ba12277fc35b87653c8cf8ec13b299817dc5 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -637,7 +637,6 @@ void Simulation::run()
 		{
 			KernelPMCumOneCompSP27( para->getParD(0)->numberofthreads,
 									para->getParD(0)->omega,			
-									para->getParD(0)->geoSP, 
 									para->getParD(0)->neighborX_SP, 
 									para->getParD(0)->neighborY_SP, 
 									para->getParD(0)->neighborZ_SP,
@@ -648,7 +647,8 @@ void Simulation::run()
 									pm0->getPorosity(),
 									pm0->getDarcyLBM(),
 									pm0->getForchheimerLBM(),
-									pm0->getGeoID(),
+									pm0->getSizePM(),
+									pm0->getHostNodeIDsPM(),
 									para->getParD(0)->evenOrOdd); 
 			getLastCudaError("KernelPMCumOneCompSP27 execution failed");
 		}
@@ -1333,12 +1333,12 @@ void Simulation::run()
 			//getLastCudaError("QDevComp27 (Geom) execution failed");
 
 
-		//QPressDevOld27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
-		//				para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
-		//				para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
-		//				para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-		//				para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-		//getLastCudaError("QPressDev27 execution failed");
+		QPressDevOld27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
+						para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
+						para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
+						para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+						para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+		getLastCudaError("QPressDev27 execution failed");
 
 		//QPressDevEQZ27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
 		//				para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
@@ -1480,12 +1480,12 @@ void Simulation::run()
 		  //getLastCudaError("QPressDevIncompNEQ27 execution failed");
 		  //////////////////////////////////////////////////////////////////////////////////
 		  //press NEQ comp
-		  QPressDevNEQ27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
-		  				  para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
-		  				  para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
-		  				  para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
-		  				  para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
-		  getLastCudaError("QPressDevNEQ27 execution failed");
+		  //QPressDevNEQ27( para->getParD(0)->numberofthreads, para->getParD(0)->QPress.RhoBC, 
+		  //				  para->getParD(0)->d0SP.f[0],       para->getParD(0)->QPress.k,  
+		  //				  para->getParD(0)->QPress.kN,       para->getParD(0)->QPress.kQ,    para->getParD(0)->omega,
+		  //				  para->getParD(0)->neighborX_SP,    para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+		  //				  para->getParD(0)->size_Mat_SP,     para->getParD(0)->evenOrOdd);
+		  //getLastCudaError("QPressDevNEQ27 execution failed");
 		  ////////////////////////////////////////////////////////////////////////////////
           //if (  myid == numprocs - 1)  
           //   PressSchlaffer27( para->getParD(0)->numberofthreads,  para->getParD(0)->Qoutflow.RhoBC,
@@ -2622,21 +2622,93 @@ void Simulation::run()
 }
 void Simulation::porousMedia()
 {
-	//Kondensator = porous media 0
-	double porosity = 0.7;
-	double darcySI = 137.36; //[1/s]
-	double forchheimerSI = 1037.8; //[1/m]
-	double dxLBM = 0.0004;
+	double porosity, darcySI, forchheimerSI;
+	double dxLBM = 0.004;
 	double dtLBM = 0.000007;
-	unsigned int level = 0;
-	output << "\nnew instance of PM \n";
-	pm0 = new PorousMedia(porosity, GEO_PM_0, darcySI, forchheimerSI, dxLBM, dtLBM, level);
-	output << "setStartCoordinates \n";
-	pm0->setStartCoordinates(20, 0, 0);
-	output << "setEndCoordinates \n";
-	pm0->setEndCoordinates(40, 20, 20);
-	output << "definePMarea \n";
+	unsigned int level, geo;
+	double startX, startY, startZ, endX, endY, endZ;
+	//////////////////////////////////////////////////////////////////////////
+
+	//////////////////////////////////////////////////////////////////////////
+	//Test = porous media 0
+	porosity = 0.7;
+	darcySI = 137.36; //[1/s]
+	forchheimerSI = 1037.8; //[1/m]
+	level = para->getFine();
+	geo = GEO_PM_0;
+	startX = 20.0;
+	startY =  0.0;
+	startZ =  0.0;
+	endX = 40.0;
+	endY = 22.0;
+	endZ = 22.0;
+	pm0 = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+	pm0->setStartCoordinates(startX, startY, startZ);
+	pm0->setEndCoordinates(endX, endY, endZ);
+	pm0->setResistanceLBM();
 	definePMarea(pm0);
+	//////////////////////////////////////////////////////////////////////////
+
+	////////////////////////////////////////////////////////////////////////////
+	////Kondensator = porous media 0
+	//porosity = 0.7;
+	//darcySI = 137.36; //[1/s]
+	//forchheimerSI = 1037.8; //[1/m]
+	//level = para->getFine();
+	//geo = GEO_PM_0;
+	//startX = -0.715882;
+	//startY = -0.260942;
+	//startZ = -0.031321;
+	//endX = -0.692484;
+	//endY =  0.277833;
+	//endZ =  0.360379;
+	//pm0 = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+	//pm0->setStartCoordinates(startX, startY, startZ);
+	//pm0->setEndCoordinates(endX, endY, endZ);
+	//pm0->setResistanceLBM();
+	//definePMarea(pm0);
+	////////////////////////////////////////////////////////////////////////////
+
+	////////////////////////////////////////////////////////////////////////////
+	////NT-Kuehler = porous media 1
+	//porosity = 0.6;
+	//darcySI = 149.98; //[1/s]
+	//forchheimerSI = 960.57; //[1/m]
+	//level = para->getFine();
+	//geo = GEO_PM_1;
+	//startX = -0.696146;
+	//startY = -0.32426;
+	//startZ = -0.0421345;
+	//endX = -0.651847;
+	//endY =  0.324822;
+	//endZ =  0.057098;
+	//pm1 = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+	//pm1->setStartCoordinates(startX, startY, startZ);
+	//pm1->setEndCoordinates(endX, endY, endZ);
+	//pm1->setResistanceLBM();
+	//definePMarea(pm1);
+	////////////////////////////////////////////////////////////////////////////
+
+	////////////////////////////////////////////////////////////////////////////
+	////Wasserkuehler = porous media 2
+	//porosity = 0.6;
+	//darcySI = 148.69; //[1/s]
+	//forchheimerSI = 629.45; //[1/m]
+	//level = para->getFine();
+	//geo = GEO_PM_2;
+	//startX = -0.692681;
+	//startY = -0.324954;
+	//startZ = 0.0789429;
+	//endX = -0.657262;
+	//endY =  0.32538;
+	//endZ =  0.400974;
+	//pm2 = new PorousMedia(porosity, geo, darcySI, forchheimerSI, dxLBM, dtLBM, level);
+	//pm2->setStartCoordinates(startX, startY, startZ);
+	//pm2->setEndCoordinates(endX, endY, endZ);
+	//pm2->setResistanceLBM();
+	//definePMarea(pm2);
+	////////////////////////////////////////////////////////////////////////////
+
 }
 void Simulation::definePMarea(PorousMedia* pm)
 {
@@ -2647,15 +2719,12 @@ void Simulation::definePMarea(PorousMedia* pm)
 
 	for (unsigned int i = 0; i < para->getParH(level)->size_Mat_SP; i++)
 	{
-		//output << "definePMarea....find nodes...for loop \n";
 		if (((para->getParH(level)->coordX_SP[i] >= pm->getStartX()) && (para->getParH(level)->coordX_SP[i] <= pm->getEndX())) &&
 			((para->getParH(level)->coordY_SP[i] >= pm->getStartY()) && (para->getParH(level)->coordY_SP[i] <= pm->getEndY())) &&
 			((para->getParH(level)->coordZ_SP[i] >= pm->getStartZ()) && (para->getParH(level)->coordZ_SP[i] <= pm->getEndZ())) )
 		{
-			//output << "definePMarea....find nodes...if area \n";
 			if (para->getParH(level)->geoSP[i] >= GEO_FLUID)
 			{
-				//output << "definePMarea....find nodes...if area and fluid \n";
 				para->getParH(level)->geoSP[i] = pm->getGeoID();
 				nodeIDsPorousMedia.push_back(i);
 				counter++;
@@ -2663,20 +2732,16 @@ void Simulation::definePMarea(PorousMedia* pm)
 		}
 	}
 
-	output << "definePMarea....cuda alloc \n";
+	para->cudaCopySP(level);
 	pm->setSizePM(counter);
 	para->cudaAllocPorousMedia(pm0, level);
 	unsigned int *tpmArrayIDs = pm->getHostNodeIDsPM();
-
-	output << "definePMarea....copy vector to array \n";
+	
 	for (unsigned int j = 0; j <= pm->getSizePM(); j++)
 	{
 		tpmArrayIDs[j] = nodeIDsPorousMedia[j];
 	}
-
-	output << "definePMarea....setHostNodeIDsPM \n";
+	
 	pm->setHostNodeIDsPM(tpmArrayIDs);
-
-	output << "definePMarea....cudaCopyPorousMedia \n";
 	para->cudaCopyPorousMedia(pm, level);
 }
diff --git a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
index bc3089bf6af180b71bcf1451ed2d7df67d68cf96..e0d3a91a924d85994707ceb00bf0151abdc95ea0 100644
--- a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
+++ b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp
@@ -498,7 +498,7 @@ namespace UnstrucuredGridWriter
 
 			for (unsigned int pos = startpos; pos < endpos; pos++)
 			{
-				if (/*para->getParH(level)->geoSP[pos] >= GEO_FLUID*/true)
+				if (/*(para->getParH(level)->geoSP[pos] >= GEO_FLUID) || ((para->getParH(level)->geoSP[pos] >= GEO_PM_0) && (para->getParH(level)->geoSP[pos] <= GEO_PM_2))*/true)
 				{
 					//////////////////////////////////////////////////////////////////////////
 					double x1 = para->getParH(level)->coordX_SP[pos];
@@ -528,13 +528,13 @@ namespace UnstrucuredGridWriter
 					number8 = para->getParH(level)->neighborZ_SP[number4];
 					//////////////////////////////////////////////////////////////////////////
 					//printf("\n test vor neighborsFluid... \n");
-					//if (para->getParH(level)->geoSP[number2] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number3] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number4] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number5] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number6] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number7] < GEO_FLUID ||
-					//	para->getParH(level)->geoSP[number8] < GEO_FLUID)  neighborsFluid = false;
+					if (((para->getParH(level)->geoSP[number2] != GEO_FLUID) && (para->getParH(level)->geoSP[number2] < GEO_PM_0) && (para->getParH(level)->geoSP[number2] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number3] != GEO_FLUID) && (para->getParH(level)->geoSP[number3] < GEO_PM_0) && (para->getParH(level)->geoSP[number3] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number4] != GEO_FLUID) && (para->getParH(level)->geoSP[number4] < GEO_PM_0) && (para->getParH(level)->geoSP[number4] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number5] != GEO_FLUID) && (para->getParH(level)->geoSP[number5] < GEO_PM_0) && (para->getParH(level)->geoSP[number5] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number6] != GEO_FLUID) && (para->getParH(level)->geoSP[number6] < GEO_PM_0) && (para->getParH(level)->geoSP[number6] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number7] != GEO_FLUID) && (para->getParH(level)->geoSP[number7] < GEO_PM_0) && (para->getParH(level)->geoSP[number7] > GEO_PM_2)) ||
+						((para->getParH(level)->geoSP[number8] != GEO_FLUID) && (para->getParH(level)->geoSP[number8] < GEO_PM_0) && (para->getParH(level)->geoSP[number8] > GEO_PM_2)))  neighborsFluid = false;
 					//////////////////////////////////////////////////////////////////////////
 					//if(neighborsFluid==false) counter++;
 					//////////////////////////////////////////////////////////////////////////