diff --git a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
index 969654fd013b30b0ce6d40ad6eeeea1e3b315a77..d862c857550b02ef7e6b59a99da5e7e512594c1e 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.cpp
@@ -109,6 +109,8 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolation(const D3
vf::lbm::calculateMomentsOnSourceNodes(icell.TSE, omegaF, drho_PMP, vx1_PMP, vx2_PMP, vx3_PMP,kxyFromfcNEQ_PMP, kyzFromfcNEQ_PMP, kxzFromfcNEQ_PMP, kxxMyyFromfcNEQ_PMP, kxxMzzFromfcNEQ_PMP);
vf::lbm::calculateMomentsOnSourceNodes(icell.BNE, omegaF, drho_PPM, vx1_PPM, vx2_PPM, vx3_PPM,kxyFromfcNEQ_PPM, kyzFromfcNEQ_PPM, kxzFromfcNEQ_PPM, kxxMyyFromfcNEQ_PPM, kxxMzzFromfcNEQ_PPM);
+
+
vf::lbm::interpolate_fc(icellC,
eps_new,
omegaC,
@@ -135,6 +137,406 @@ void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolation(const D3
kxyFromfcNEQ_PMM, kyzFromfcNEQ_PMM, kxzFromfcNEQ_PMM, kxxMyyFromfcNEQ_PMM, kxxMzzFromfcNEQ_PMM,
kxyFromfcNEQ_MMM, kyzFromfcNEQ_MMM, kxzFromfcNEQ_MMM, kxxMyyFromfcNEQ_MMM, kxxMzzFromfcNEQ_MMM);
}
+//////////////////////////////////////////////////////////////////////////
+void CompressibleOffsetMomentsInterpolationProcessor::calcMoments(const real* const f, real omega, real& press, real& vx1, real& vx2, real& vx3,
+ real& kxy, real& kyz, real& kxz, real& kxxMyy, real& kxxMzz)
+{
+ using namespace D3Q27System;
+ using namespace vf::lbm::dir;
+
+ real drho = 0.0;
+ D3Q27System::calcCompMacroscopicValues(f,drho,vx1,vx2,vx3);
+
+ press = drho; //interpolate rho!
+
+ kxy = -3.*omega*((((f[DIR_MMP]+f[DIR_PPM])-(f[DIR_MPP]+f[DIR_PMM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_MPM]+f[DIR_PMP])))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_MP0]+f[DIR_PM0]))/(c1o1 + drho)-(vx1*vx2));// might not be optimal MG 25.2.13
+ kyz = -3.*omega*((((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMP]+f[DIR_MPM]))+((f[DIR_PMM]+f[DIR_MPP])-(f[DIR_MMP]+f[DIR_PPM])))+((f[DIR_0MM]+f[DIR_0PP])-(f[DIR_0MP]+f[DIR_0PM]))/(c1o1 + drho)-(vx2*vx3));
+ kxz = -3.*omega*((((f[DIR_MPM]+f[DIR_PMP])-(f[DIR_MMP]+f[DIR_PPM]))+((f[DIR_MMM]+f[DIR_PPP])-(f[DIR_PMM]+f[DIR_MPP])))+((f[DIR_M0M]+f[DIR_P0P])-(f[DIR_M0P]+f[DIR_P0M]))/(c1o1 + drho)-(vx1*vx3));
+ kxxMyy = -3./2.*omega*((((f[DIR_M0M]+f[DIR_P0P])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_M0P]+f[DIR_P0M])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_0M0]+f[DIR_0P0]))/(c1o1 + drho)-(vx1*vx1-vx2*vx2));
+ kxxMzz = -3./2.*omega*((((f[DIR_MP0]+f[DIR_PM0])-(f[DIR_0MM]+f[DIR_0PP]))+((f[DIR_MM0]+f[DIR_PP0])-(f[DIR_0MP]+f[DIR_0PM])))+((f[DIR_M00]+f[DIR_P00])-(f[DIR_00M]+f[DIR_00P]))/(c1o1 + drho)-(vx1*vx1-vx3*vx3));
+}
+//////////////////////////////////////////////////////////////////////////
+void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedCoefficiets(const D3Q27ICell& icell, real omega, real eps_new)
+{
+ real vx1_SWT,vx2_SWT,vx3_SWT;
+ real vx1_NWT,vx2_NWT,vx3_NWT;
+ real vx1_NET,vx2_NET,vx3_NET;
+ real vx1_SET,vx2_SET,vx3_SET;
+ real vx1_SWB,vx2_SWB,vx3_SWB;
+ real vx1_NWB,vx2_NWB,vx3_NWB;
+ real vx1_NEB,vx2_NEB,vx3_NEB;
+ real vx1_SEB,vx2_SEB,vx3_SEB;
+
+ real kxyFromfcNEQ_SWT, kyzFromfcNEQ_SWT, kxzFromfcNEQ_SWT, kxxMyyFromfcNEQ_SWT, kxxMzzFromfcNEQ_SWT;
+ real kxyFromfcNEQ_NWT, kyzFromfcNEQ_NWT, kxzFromfcNEQ_NWT, kxxMyyFromfcNEQ_NWT, kxxMzzFromfcNEQ_NWT;
+ real kxyFromfcNEQ_NET, kyzFromfcNEQ_NET, kxzFromfcNEQ_NET, kxxMyyFromfcNEQ_NET, kxxMzzFromfcNEQ_NET;
+ real kxyFromfcNEQ_SET, kyzFromfcNEQ_SET, kxzFromfcNEQ_SET, kxxMyyFromfcNEQ_SET, kxxMzzFromfcNEQ_SET;
+ real kxyFromfcNEQ_SWB, kyzFromfcNEQ_SWB, kxzFromfcNEQ_SWB, kxxMyyFromfcNEQ_SWB, kxxMzzFromfcNEQ_SWB;
+ real kxyFromfcNEQ_NWB, kyzFromfcNEQ_NWB, kxzFromfcNEQ_NWB, kxxMyyFromfcNEQ_NWB, kxxMzzFromfcNEQ_NWB;
+ real kxyFromfcNEQ_NEB, kyzFromfcNEQ_NEB, kxzFromfcNEQ_NEB, kxxMyyFromfcNEQ_NEB, kxxMzzFromfcNEQ_NEB;
+ real kxyFromfcNEQ_SEB, kyzFromfcNEQ_SEB, kxzFromfcNEQ_SEB, kxxMyyFromfcNEQ_SEB, kxxMzzFromfcNEQ_SEB;
+
+ calcMoments(icell.TSW,omega,press_SWT,vx1_SWT,vx2_SWT,vx3_SWT, kxyFromfcNEQ_SWT, kyzFromfcNEQ_SWT, kxzFromfcNEQ_SWT, kxxMyyFromfcNEQ_SWT, kxxMzzFromfcNEQ_SWT);
+ calcMoments(icell.TNW,omega,press_NWT,vx1_NWT,vx2_NWT,vx3_NWT, kxyFromfcNEQ_NWT, kyzFromfcNEQ_NWT, kxzFromfcNEQ_NWT, kxxMyyFromfcNEQ_NWT, kxxMzzFromfcNEQ_NWT);
+ calcMoments(icell.TNE,omega,press_NET,vx1_NET,vx2_NET,vx3_NET, kxyFromfcNEQ_NET, kyzFromfcNEQ_NET, kxzFromfcNEQ_NET, kxxMyyFromfcNEQ_NET, kxxMzzFromfcNEQ_NET);
+ calcMoments(icell.TSE,omega,press_SET,vx1_SET,vx2_SET,vx3_SET, kxyFromfcNEQ_SET, kyzFromfcNEQ_SET, kxzFromfcNEQ_SET, kxxMyyFromfcNEQ_SET, kxxMzzFromfcNEQ_SET);
+ calcMoments(icell.BSW,omega,press_SWB,vx1_SWB,vx2_SWB,vx3_SWB, kxyFromfcNEQ_SWB, kyzFromfcNEQ_SWB, kxzFromfcNEQ_SWB, kxxMyyFromfcNEQ_SWB, kxxMzzFromfcNEQ_SWB);
+ calcMoments(icell.BNW,omega,press_NWB,vx1_NWB,vx2_NWB,vx3_NWB, kxyFromfcNEQ_NWB, kyzFromfcNEQ_NWB, kxzFromfcNEQ_NWB, kxxMyyFromfcNEQ_NWB, kxxMzzFromfcNEQ_NWB);
+ calcMoments(icell.BNE,omega,press_NEB,vx1_NEB,vx2_NEB,vx3_NEB, kxyFromfcNEQ_NEB, kyzFromfcNEQ_NEB, kxzFromfcNEQ_NEB, kxxMyyFromfcNEQ_NEB, kxxMzzFromfcNEQ_NEB);
+ calcMoments(icell.BSE,omega,press_SEB,vx1_SEB,vx2_SEB,vx3_SEB, kxyFromfcNEQ_SEB, kyzFromfcNEQ_SEB, kxzFromfcNEQ_SEB, kxxMyyFromfcNEQ_SEB, kxxMzzFromfcNEQ_SEB);
+
+ //LBMReal dxRho=c1o4*((press_NET-press_SWB)+(press_SET-press_NWB)+(press_NEB-press_SWT)+(press_SEB-press_NWT));
+ //LBMReal dyRho=c1o4*((press_NET-press_SWB)-(press_SET-press_NWB)+(press_NEB-press_SWT)-(press_SEB-press_NWT));
+ //LBMReal dzRho=c1o4*((press_NET-press_SWB)+(press_SET-press_NWB)-(press_NEB-press_SWT)-(press_SEB-press_NWT));
+
+ // kxyFromfcNEQ_SWT+=vx1_SWT*dyRho+vx2_SWT*dxRho;
+ // kxyFromfcNEQ_NWT+=vx1_NWT*dyRho+vx2_NWT*dxRho;
+ // kxyFromfcNEQ_NET+=vx1_NET*dyRho+vx2_NET*dxRho;
+ // kxyFromfcNEQ_SET+=vx1_SET*dyRho+vx2_SET*dxRho;
+ // kxyFromfcNEQ_SWB+=vx1_SWB*dyRho+vx2_SWB*dxRho;
+ // kxyFromfcNEQ_NWB+=vx1_NWB*dyRho+vx2_NWB*dxRho;
+ // kxyFromfcNEQ_NEB+=vx1_NEB*dyRho+vx2_NEB*dxRho;
+ // kxyFromfcNEQ_SEB+=vx1_SEB*dyRho+vx2_SEB*dxRho;
+
+ // kyzFromfcNEQ_SWT+=vx3_SWT*dyRho+vx2_SWT*dzRho;
+ // kyzFromfcNEQ_NWT+=vx3_NWT*dyRho+vx2_NWT*dzRho;
+ // kyzFromfcNEQ_NET+=vx3_NET*dyRho+vx2_NET*dzRho;
+ // kyzFromfcNEQ_SET+=vx3_SET*dyRho+vx2_SET*dzRho;
+ // kyzFromfcNEQ_SWB+=vx3_SWB*dyRho+vx2_SWB*dzRho;
+ // kyzFromfcNEQ_NWB+=vx3_NWB*dyRho+vx2_NWB*dzRho;
+ // kyzFromfcNEQ_NEB+=vx3_NEB*dyRho+vx2_NEB*dzRho;
+ // kyzFromfcNEQ_SEB+=vx3_SEB*dyRho+vx2_SEB*dzRho;
+
+ // kxzFromfcNEQ_SWT+=vx1_SWT*dzRho+vx3_SWT*dxRho;
+ // kxzFromfcNEQ_NWT+=vx1_NWT*dzRho+vx3_NWT*dxRho;
+ // kxzFromfcNEQ_NET+=vx1_NET*dzRho+vx3_NET*dxRho;
+ // kxzFromfcNEQ_SET+=vx1_SET*dzRho+vx3_SET*dxRho;
+ // kxzFromfcNEQ_SWB+=vx1_SWB*dzRho+vx3_SWB*dxRho;
+ // kxzFromfcNEQ_NWB+=vx1_NWB*dzRho+vx3_NWB*dxRho;
+ // kxzFromfcNEQ_NEB+=vx1_NEB*dzRho+vx3_NEB*dxRho;
+ // kxzFromfcNEQ_SEB+=vx1_SEB*dzRho+vx3_SEB*dxRho;
+
+ // kxxMyyFromfcNEQ_SWT+=vx1_SWT*dxRho-vx2_SWT*dyRho;
+ // kxxMyyFromfcNEQ_NWT+=vx1_NWT*dxRho-vx2_NWT*dyRho;
+ // kxxMyyFromfcNEQ_NET+=vx1_NET*dxRho-vx2_NET*dyRho;
+ // kxxMyyFromfcNEQ_SET+=vx1_SET*dxRho-vx2_SET*dyRho;
+ // kxxMyyFromfcNEQ_SWB+=vx1_SWB*dxRho-vx2_SWB*dyRho;
+ // kxxMyyFromfcNEQ_NWB+=vx1_NWB*dxRho-vx2_NWB*dyRho;
+ // kxxMyyFromfcNEQ_NEB+=vx1_NEB*dxRho-vx2_NEB*dyRho;
+ // kxxMyyFromfcNEQ_SEB+=vx1_SEB*dxRho-vx2_SEB*dyRho;
+
+ // kxxMzzFromfcNEQ_SWT+=vx1_SWT*dxRho-vx3_SWT*dzRho;
+ // kxxMzzFromfcNEQ_NWT+=vx1_NWT*dxRho-vx3_NWT*dzRho;
+ // kxxMzzFromfcNEQ_NET+=vx1_NET*dxRho-vx3_NET*dzRho;
+ // kxxMzzFromfcNEQ_SET+=vx1_SET*dxRho-vx3_SET*dzRho;
+ // kxxMzzFromfcNEQ_SWB+=vx1_SWB*dxRho-vx3_SWB*dzRho;
+ // kxxMzzFromfcNEQ_NWB+=vx1_NWB*dxRho-vx3_NWB*dzRho;
+ // kxxMzzFromfcNEQ_NEB+=vx1_NEB*dxRho-vx3_NEB*dzRho;
+ // kxxMzzFromfcNEQ_SEB+=vx1_SEB*dxRho-vx3_SEB*dzRho;
+
+
+ //kxxMzzFromfcNEQ_SWT=0.0;
+ //kxxMzzFromfcNEQ_NWT=0.0;
+ //kxxMzzFromfcNEQ_NET=0.0;
+ //kxxMzzFromfcNEQ_SET=0.0;
+ //kxxMzzFromfcNEQ_SWB=0.0;
+ //kxxMzzFromfcNEQ_NWB=0.0;
+ //kxxMzzFromfcNEQ_NEB=0.0;
+ //kxxMzzFromfcNEQ_SEB=0.0;
+
+
+
+
+
+ a0 = (-kxxMyyFromfcNEQ_NEB - kxxMyyFromfcNEQ_NET + kxxMyyFromfcNEQ_NWB + kxxMyyFromfcNEQ_NWT -
+ kxxMyyFromfcNEQ_SEB - kxxMyyFromfcNEQ_SET + kxxMyyFromfcNEQ_SWB + kxxMyyFromfcNEQ_SWT -
+ kxxMzzFromfcNEQ_NEB - kxxMzzFromfcNEQ_NET + kxxMzzFromfcNEQ_NWB + kxxMzzFromfcNEQ_NWT -
+ kxxMzzFromfcNEQ_SEB - kxxMzzFromfcNEQ_SET + kxxMzzFromfcNEQ_SWB + kxxMzzFromfcNEQ_SWT -
+ 2.*kxyFromfcNEQ_NEB - 2.*kxyFromfcNEQ_NET - 2.*kxyFromfcNEQ_NWB - 2.*kxyFromfcNEQ_NWT +
+ 2.*kxyFromfcNEQ_SEB + 2.*kxyFromfcNEQ_SET + 2.*kxyFromfcNEQ_SWB + 2.*kxyFromfcNEQ_SWT +
+ 2.*kxzFromfcNEQ_NEB - 2.*kxzFromfcNEQ_NET + 2.*kxzFromfcNEQ_NWB - 2.*kxzFromfcNEQ_NWT +
+ 2.*kxzFromfcNEQ_SEB - 2.*kxzFromfcNEQ_SET + 2.*kxzFromfcNEQ_SWB - 2.*kxzFromfcNEQ_SWT +
+ 8.*vx1_NEB + 8.*vx1_NET + 8.*vx1_NWB + 8.*vx1_NWT + 8.*vx1_SEB +
+ 8.*vx1_SET + 8.*vx1_SWB + 8.*vx1_SWT + 2.*vx2_NEB + 2.*vx2_NET -
+ 2.*vx2_NWB - 2.*vx2_NWT - 2.*vx2_SEB - 2.*vx2_SET + 2.*vx2_SWB +
+ 2.*vx2_SWT - 2.*vx3_NEB + 2.*vx3_NET + 2.*vx3_NWB - 2.*vx3_NWT -
+ 2.*vx3_SEB + 2.*vx3_SET + 2.*vx3_SWB - 2.*vx3_SWT)/64.;
+ b0 = (2.*kxxMyyFromfcNEQ_NEB + 2.*kxxMyyFromfcNEQ_NET + 2.*kxxMyyFromfcNEQ_NWB + 2.*kxxMyyFromfcNEQ_NWT -
+ 2.*kxxMyyFromfcNEQ_SEB - 2.*kxxMyyFromfcNEQ_SET - 2.*kxxMyyFromfcNEQ_SWB - 2.*kxxMyyFromfcNEQ_SWT -
+ kxxMzzFromfcNEQ_NEB - kxxMzzFromfcNEQ_NET - kxxMzzFromfcNEQ_NWB - kxxMzzFromfcNEQ_NWT +
+ kxxMzzFromfcNEQ_SEB + kxxMzzFromfcNEQ_SET + kxxMzzFromfcNEQ_SWB + kxxMzzFromfcNEQ_SWT -
+ 2.*kxyFromfcNEQ_NEB - 2.*kxyFromfcNEQ_NET + 2.*kxyFromfcNEQ_NWB + 2.*kxyFromfcNEQ_NWT -
+ 2.*kxyFromfcNEQ_SEB - 2.*kxyFromfcNEQ_SET + 2.*kxyFromfcNEQ_SWB + 2.*kxyFromfcNEQ_SWT +
+ 2.*kyzFromfcNEQ_NEB - 2.*kyzFromfcNEQ_NET + 2.*kyzFromfcNEQ_NWB - 2.*kyzFromfcNEQ_NWT +
+ 2.*kyzFromfcNEQ_SEB - 2.*kyzFromfcNEQ_SET + 2.*kyzFromfcNEQ_SWB - 2.*kyzFromfcNEQ_SWT +
+ 2.*vx1_NEB + 2.*vx1_NET - 2.*vx1_NWB - 2.*vx1_NWT -
+ 2.*vx1_SEB - 2.*vx1_SET + 2.*vx1_SWB + 2.*vx1_SWT +
+ 8.*vx2_NEB + 8.*vx2_NET + 8.*vx2_NWB + 8.*vx2_NWT +
+ 8.*vx2_SEB + 8.*vx2_SET + 8.*vx2_SWB + 8.*vx2_SWT -
+ 2.*vx3_NEB + 2.*vx3_NET - 2.*vx3_NWB + 2.*vx3_NWT +
+ 2.*vx3_SEB - 2.*vx3_SET + 2.*vx3_SWB - 2.*vx3_SWT)/64.;
+ c0 = (kxxMyyFromfcNEQ_NEB - kxxMyyFromfcNEQ_NET + kxxMyyFromfcNEQ_NWB - kxxMyyFromfcNEQ_NWT +
+ kxxMyyFromfcNEQ_SEB - kxxMyyFromfcNEQ_SET + kxxMyyFromfcNEQ_SWB - kxxMyyFromfcNEQ_SWT -
+ 2.*kxxMzzFromfcNEQ_NEB + 2.*kxxMzzFromfcNEQ_NET - 2.*kxxMzzFromfcNEQ_NWB + 2.*kxxMzzFromfcNEQ_NWT -
+ 2.*kxxMzzFromfcNEQ_SEB + 2.*kxxMzzFromfcNEQ_SET - 2.*kxxMzzFromfcNEQ_SWB + 2.*kxxMzzFromfcNEQ_SWT -
+ 2.*kxzFromfcNEQ_NEB - 2.*kxzFromfcNEQ_NET + 2.*kxzFromfcNEQ_NWB + 2.*kxzFromfcNEQ_NWT -
+ 2.*kxzFromfcNEQ_SEB - 2.*kxzFromfcNEQ_SET + 2.*kxzFromfcNEQ_SWB + 2.*kxzFromfcNEQ_SWT -
+ 2.*kyzFromfcNEQ_NEB - 2.*kyzFromfcNEQ_NET - 2.*kyzFromfcNEQ_NWB - 2.*kyzFromfcNEQ_NWT +
+ 2.*kyzFromfcNEQ_SEB + 2.*kyzFromfcNEQ_SET + 2.*kyzFromfcNEQ_SWB + 2.*kyzFromfcNEQ_SWT -
+ 2.*vx1_NEB + 2.*vx1_NET + 2.*vx1_NWB - 2.*vx1_NWT -
+ 2.*vx1_SEB + 2.*vx1_SET + 2.*vx1_SWB - 2.*vx1_SWT -
+ 2.*vx2_NEB + 2.*vx2_NET - 2.*vx2_NWB + 2.*vx2_NWT +
+ 2.*vx2_SEB - 2.*vx2_SET + 2.*vx2_SWB - 2.*vx2_SWT +
+ 8.*vx3_NEB + 8.*vx3_NET + 8.*vx3_NWB + 8.*vx3_NWT +
+ 8.*vx3_SEB + 8.*vx3_SET + 8.*vx3_SWB + 8.*vx3_SWT)/64.;
+ ax = (vx1_NEB + vx1_NET - vx1_NWB - vx1_NWT + vx1_SEB + vx1_SET - vx1_SWB - vx1_SWT)/4.;
+ bx = (vx2_NEB + vx2_NET - vx2_NWB - vx2_NWT + vx2_SEB + vx2_SET - vx2_SWB - vx2_SWT)/4.;
+ cx = (vx3_NEB + vx3_NET - vx3_NWB - vx3_NWT + vx3_SEB + vx3_SET - vx3_SWB - vx3_SWT)/4.;
+ axx= (kxxMyyFromfcNEQ_NEB + kxxMyyFromfcNEQ_NET - kxxMyyFromfcNEQ_NWB - kxxMyyFromfcNEQ_NWT +
+ kxxMyyFromfcNEQ_SEB + kxxMyyFromfcNEQ_SET - kxxMyyFromfcNEQ_SWB - kxxMyyFromfcNEQ_SWT +
+ kxxMzzFromfcNEQ_NEB + kxxMzzFromfcNEQ_NET - kxxMzzFromfcNEQ_NWB - kxxMzzFromfcNEQ_NWT +
+ kxxMzzFromfcNEQ_SEB + kxxMzzFromfcNEQ_SET - kxxMzzFromfcNEQ_SWB - kxxMzzFromfcNEQ_SWT +
+ 2.*vx2_NEB + 2.*vx2_NET - 2.*vx2_NWB - 2.*vx2_NWT -
+ 2.*vx2_SEB - 2.*vx2_SET + 2.*vx2_SWB + 2.*vx2_SWT -
+ 2.*vx3_NEB + 2.*vx3_NET + 2.*vx3_NWB - 2.*vx3_NWT -
+ 2.*vx3_SEB + 2.*vx3_SET + 2.*vx3_SWB - 2.*vx3_SWT)/16.;
+ bxx= (kxyFromfcNEQ_NEB + kxyFromfcNEQ_NET - kxyFromfcNEQ_NWB - kxyFromfcNEQ_NWT +
+ kxyFromfcNEQ_SEB + kxyFromfcNEQ_SET - kxyFromfcNEQ_SWB - kxyFromfcNEQ_SWT -
+ 2.*vx1_NEB - 2.*vx1_NET + 2.*vx1_NWB + 2.*vx1_NWT +
+ 2.*vx1_SEB + 2.*vx1_SET - 2.*vx1_SWB - 2.*vx1_SWT)/8.;
+ cxx= (kxzFromfcNEQ_NEB + kxzFromfcNEQ_NET - kxzFromfcNEQ_NWB - kxzFromfcNEQ_NWT +
+ kxzFromfcNEQ_SEB + kxzFromfcNEQ_SET - kxzFromfcNEQ_SWB - kxzFromfcNEQ_SWT +
+ 2.*vx1_NEB - 2.*vx1_NET - 2.*vx1_NWB + 2.*vx1_NWT +
+ 2.*vx1_SEB - 2.*vx1_SET - 2.*vx1_SWB + 2.*vx1_SWT)/8.;
+ ay = (vx1_NEB + vx1_NET + vx1_NWB + vx1_NWT - vx1_SEB - vx1_SET - vx1_SWB - vx1_SWT)/4.;
+ by = (vx2_NEB + vx2_NET + vx2_NWB + vx2_NWT - vx2_SEB - vx2_SET - vx2_SWB - vx2_SWT)/4.;
+ cy = (vx3_NEB + vx3_NET + vx3_NWB + vx3_NWT - vx3_SEB - vx3_SET - vx3_SWB - vx3_SWT)/4.;
+ ayy= (kxyFromfcNEQ_NEB + kxyFromfcNEQ_NET + kxyFromfcNEQ_NWB + kxyFromfcNEQ_NWT -
+ kxyFromfcNEQ_SEB - kxyFromfcNEQ_SET - kxyFromfcNEQ_SWB - kxyFromfcNEQ_SWT -
+ 2.*vx2_NEB - 2.*vx2_NET + 2.*vx2_NWB + 2.*vx2_NWT +
+ 2.*vx2_SEB + 2.*vx2_SET - 2.*vx2_SWB - 2.*vx2_SWT)/8.;
+ byy= (-2.*kxxMyyFromfcNEQ_NEB - 2.*kxxMyyFromfcNEQ_NET - 2.*kxxMyyFromfcNEQ_NWB - 2.*kxxMyyFromfcNEQ_NWT +
+ 2.*kxxMyyFromfcNEQ_SEB + 2.*kxxMyyFromfcNEQ_SET + 2.*kxxMyyFromfcNEQ_SWB + 2.*kxxMyyFromfcNEQ_SWT +
+ kxxMzzFromfcNEQ_NEB + kxxMzzFromfcNEQ_NET + kxxMzzFromfcNEQ_NWB + kxxMzzFromfcNEQ_NWT -
+ kxxMzzFromfcNEQ_SEB - kxxMzzFromfcNEQ_SET - kxxMzzFromfcNEQ_SWB - kxxMzzFromfcNEQ_SWT +
+ 2.*vx1_NEB + 2.*vx1_NET - 2.*vx1_NWB - 2.*vx1_NWT -
+ 2.*vx1_SEB - 2.*vx1_SET + 2.*vx1_SWB + 2.*vx1_SWT -
+ 2.*vx3_NEB + 2.*vx3_NET - 2.*vx3_NWB + 2.*vx3_NWT +
+ 2.*vx3_SEB - 2.*vx3_SET + 2.*vx3_SWB - 2.*vx3_SWT)/16.;
+ cyy= (kyzFromfcNEQ_NEB + kyzFromfcNEQ_NET + kyzFromfcNEQ_NWB + kyzFromfcNEQ_NWT -
+ kyzFromfcNEQ_SEB - kyzFromfcNEQ_SET - kyzFromfcNEQ_SWB - kyzFromfcNEQ_SWT +
+ 2.*vx2_NEB - 2.*vx2_NET + 2.*vx2_NWB - 2.*vx2_NWT -
+ 2.*vx2_SEB + 2.*vx2_SET - 2.*vx2_SWB + 2.*vx2_SWT)/8.;
+ az = (-vx1_NEB + vx1_NET - vx1_NWB + vx1_NWT - vx1_SEB + vx1_SET - vx1_SWB + vx1_SWT)/4.;
+ bz = (-vx2_NEB + vx2_NET - vx2_NWB + vx2_NWT - vx2_SEB + vx2_SET - vx2_SWB + vx2_SWT)/4.;
+ cz = (-vx3_NEB + vx3_NET - vx3_NWB + vx3_NWT - vx3_SEB + vx3_SET - vx3_SWB + vx3_SWT)/4.;
+ azz= (-kxzFromfcNEQ_NEB + kxzFromfcNEQ_NET - kxzFromfcNEQ_NWB + kxzFromfcNEQ_NWT -
+ kxzFromfcNEQ_SEB + kxzFromfcNEQ_SET - kxzFromfcNEQ_SWB + kxzFromfcNEQ_SWT +
+ 2.*vx3_NEB - 2.*vx3_NET - 2.*vx3_NWB + 2.*vx3_NWT +
+ 2.*vx3_SEB - 2.*vx3_SET - 2.*vx3_SWB + 2.*vx3_SWT)/8.;
+ bzz= (-kyzFromfcNEQ_NEB + kyzFromfcNEQ_NET - kyzFromfcNEQ_NWB + kyzFromfcNEQ_NWT -
+ kyzFromfcNEQ_SEB + kyzFromfcNEQ_SET - kyzFromfcNEQ_SWB + kyzFromfcNEQ_SWT +
+ 2.*vx3_NEB - 2.*vx3_NET + 2.*vx3_NWB - 2.*vx3_NWT -
+ 2.*vx3_SEB + 2.*vx3_SET - 2.*vx3_SWB + 2.*vx3_SWT)/8.;
+ czz= (-kxxMyyFromfcNEQ_NEB + kxxMyyFromfcNEQ_NET - kxxMyyFromfcNEQ_NWB + kxxMyyFromfcNEQ_NWT -
+ kxxMyyFromfcNEQ_SEB + kxxMyyFromfcNEQ_SET - kxxMyyFromfcNEQ_SWB + kxxMyyFromfcNEQ_SWT +
+ 2.*kxxMzzFromfcNEQ_NEB - 2.*kxxMzzFromfcNEQ_NET + 2.*kxxMzzFromfcNEQ_NWB - 2.*kxxMzzFromfcNEQ_NWT +
+ 2.*kxxMzzFromfcNEQ_SEB - 2.*kxxMzzFromfcNEQ_SET + 2.*kxxMzzFromfcNEQ_SWB - 2.*kxxMzzFromfcNEQ_SWT -
+ 2.*vx1_NEB + 2.*vx1_NET + 2.*vx1_NWB - 2.*vx1_NWT -
+ 2.*vx1_SEB + 2.*vx1_SET + 2.*vx1_SWB - 2.*vx1_SWT -
+ 2.*vx2_NEB + 2.*vx2_NET - 2.*vx2_NWB + 2.*vx2_NWT +
+ 2.*vx2_SEB - 2.*vx2_SET + 2.*vx2_SWB - 2.*vx2_SWT)/16.;
+ axy= (vx1_NEB + vx1_NET - vx1_NWB - vx1_NWT - vx1_SEB - vx1_SET + vx1_SWB + vx1_SWT)/2.;
+ bxy= (vx2_NEB + vx2_NET - vx2_NWB - vx2_NWT - vx2_SEB - vx2_SET + vx2_SWB + vx2_SWT)/2.;
+ cxy= (vx3_NEB + vx3_NET - vx3_NWB - vx3_NWT - vx3_SEB - vx3_SET + vx3_SWB + vx3_SWT)/2.;
+ axz= (-vx1_NEB + vx1_NET + vx1_NWB - vx1_NWT - vx1_SEB + vx1_SET + vx1_SWB - vx1_SWT)/2.;
+ bxz= (-vx2_NEB + vx2_NET + vx2_NWB - vx2_NWT - vx2_SEB + vx2_SET + vx2_SWB - vx2_SWT)/2.;
+ cxz= (-vx3_NEB + vx3_NET + vx3_NWB - vx3_NWT - vx3_SEB + vx3_SET + vx3_SWB - vx3_SWT)/2.;
+ ayz= (-vx1_NEB + vx1_NET - vx1_NWB + vx1_NWT + vx1_SEB - vx1_SET + vx1_SWB - vx1_SWT)/2.;
+ byz= (-vx2_NEB + vx2_NET - vx2_NWB + vx2_NWT + vx2_SEB - vx2_SET + vx2_SWB - vx2_SWT)/2.;
+ cyz= (-vx3_NEB + vx3_NET - vx3_NWB + vx3_NWT + vx3_SEB - vx3_SET + vx3_SWB - vx3_SWT)/2.;
+ axyz=-vx1_NEB + vx1_NET + vx1_NWB - vx1_NWT + vx1_SEB - vx1_SET - vx1_SWB + vx1_SWT;
+ bxyz=-vx2_NEB + vx2_NET + vx2_NWB - vx2_NWT + vx2_SEB - vx2_SET - vx2_SWB + vx2_SWT;
+ cxyz=-vx3_NEB + vx3_NET + vx3_NWB - vx3_NWT + vx3_SEB - vx3_SET - vx3_SWB + vx3_SWT;
+
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ kxyAverage =0;//(kxyFromfcNEQ_SWB+
+ //kxyFromfcNEQ_SWT+
+ //kxyFromfcNEQ_SET+
+ //kxyFromfcNEQ_SEB+
+ //kxyFromfcNEQ_NWB+
+ //kxyFromfcNEQ_NWT+
+ //kxyFromfcNEQ_NET+
+ //kxyFromfcNEQ_NEB)*c1o8-(ay+bx);
+ kyzAverage =0;//(kyzFromfcNEQ_SWB+
+ //kyzFromfcNEQ_SWT+
+ //kyzFromfcNEQ_SET+
+ //kyzFromfcNEQ_SEB+
+ //kyzFromfcNEQ_NWB+
+ //kyzFromfcNEQ_NWT+
+ //kyzFromfcNEQ_NET+
+ //kyzFromfcNEQ_NEB)*c1o8-(bz+cy);
+ kxzAverage =0;//(kxzFromfcNEQ_SWB+
+ //kxzFromfcNEQ_SWT+
+ //kxzFromfcNEQ_SET+
+ //kxzFromfcNEQ_SEB+
+ //kxzFromfcNEQ_NWB+
+ //kxzFromfcNEQ_NWT+
+ //kxzFromfcNEQ_NET+
+ //kxzFromfcNEQ_NEB)*c1o8-(az+cx);
+ kxxMyyAverage =0;//(kxxMyyFromfcNEQ_SWB+
+ //kxxMyyFromfcNEQ_SWT+
+ //kxxMyyFromfcNEQ_SET+
+ //kxxMyyFromfcNEQ_SEB+
+ //kxxMyyFromfcNEQ_NWB+
+ //kxxMyyFromfcNEQ_NWT+
+ //kxxMyyFromfcNEQ_NET+
+ //kxxMyyFromfcNEQ_NEB)*c1o8-(ax-by);
+ kxxMzzAverage =0;//(kxxMzzFromfcNEQ_SWB+
+ //kxxMzzFromfcNEQ_SWT+
+ //kxxMzzFromfcNEQ_SET+
+ //kxxMzzFromfcNEQ_SEB+
+ //kxxMzzFromfcNEQ_NWB+
+ //kxxMzzFromfcNEQ_NWT+
+ //kxxMzzFromfcNEQ_NET+
+ //kxxMzzFromfcNEQ_NEB)*c1o8-(ax-cz);
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //
+ // Bernd das Brot
+ //
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ a0 = a0 + xoff * ax + yoff * ay + zoff * az + xoff_sq * axx + yoff_sq * ayy + zoff_sq * azz + xoff*yoff*axy + xoff*zoff*axz + yoff*zoff*ayz + xoff*yoff*zoff*axyz ;
+ ax = ax + 2. * xoff * axx + yoff * axy + zoff * axz + yoff*zoff*axyz;
+ ay = ay + 2. * yoff * ayy + xoff * axy + zoff * ayz + xoff*zoff*axyz;
+ az = az + 2. * zoff * azz + xoff * axz + yoff * ayz + xoff*yoff*axyz;
+ b0 = b0 + xoff * bx + yoff * by + zoff * bz + xoff_sq * bxx + yoff_sq * byy + zoff_sq * bzz + xoff*yoff*bxy + xoff*zoff*bxz + yoff*zoff*byz + xoff*yoff*zoff*bxyz;
+ bx = bx + 2. * xoff * bxx + yoff * bxy + zoff * bxz + yoff*zoff*bxyz;
+ by = by + 2. * yoff * byy + xoff * bxy + zoff * byz + xoff*zoff*bxyz;
+ bz = bz + 2. * zoff * bzz + xoff * bxz + yoff * byz + xoff*yoff*bxyz;
+ c0 = c0 + xoff * cx + yoff * cy + zoff * cz + xoff_sq * cxx + yoff_sq * cyy + zoff_sq * czz + xoff*yoff*cxy + xoff*zoff*cxz + yoff*zoff*cyz + xoff*yoff*zoff*cxyz;
+ cx = cx + 2. * xoff * cxx + yoff * cxy + zoff * cxz + yoff*zoff*cxyz;
+ cy = cy + 2. * yoff * cyy + xoff * cxy + zoff * cyz + xoff*zoff*cxyz;
+ cz = cz + 2. * zoff * czz + xoff * cxz + yoff * cyz + xoff*yoff*cxyz;
+ axy= axy + zoff*axyz;
+ axz= axz + yoff*axyz;
+ ayz= ayz + xoff*axyz;
+ bxy= bxy + zoff*bxyz;
+ bxz= bxz + yoff*bxyz;
+ byz= byz + xoff*bxyz;
+ cxy= cxy + zoff*cxyz;
+ cxz= cxz + yoff*cxyz;
+ cyz= cyz + xoff*cxyz;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ const real o = omega;
+
+ f_E = eps_new*((2*(-2*ax + by + cz-kxxMzzAverage-kxxMyyAverage))/(27.*o));
+ f_N = eps_new*((2*(ax - 2*by + cz+2*kxxMyyAverage-kxxMzzAverage))/(27.*o));
+ f_T = eps_new*((2*(ax + by - 2*cz-kxxMyyAverage+2*kxxMzzAverage))/(27.*o));
+ f_NE = eps_new*(-(ax + 3*ay + 3*bx + by - 2*cz+2*kxxMyyAverage-kxxMyyAverage+3*kxyAverage)/(54.*o));
+ f_SE = eps_new*(-(ax - 3*ay - 3*bx + by - 2*cz+2*kxxMyyAverage-kxxMyyAverage-3*kxyAverage)/(54.*o));
+ f_TE = eps_new*(-(ax + 3*az - 2*by + 3*cx + cz+2*kxxMyyAverage-kxxMzzAverage+3*kxzAverage)/(54.*o));
+ f_BE = eps_new*(-(ax - 3*az - 2*by - 3*cx + cz+2*kxxMyyAverage-kxxMzzAverage-3*kxzAverage)/(54.*o));
+ f_TN = eps_new*(-(-2*ax + by + 3*bz + 3*cy + cz-kxxMyyAverage-kxxMzzAverage+3*kyzAverage)/(54.*o));
+ f_BN = eps_new*(-(-2*ax + by - 3*bz - 3*cy + cz-kxxMyyAverage-kxxMzzAverage-3*kyzAverage)/(54.*o));
+ f_ZERO = 0.;
+ f_TNE = eps_new*(-(ay + az + bx + bz + cx + cy+kxyAverage+kxzAverage+kyzAverage)/(72.*o));
+ f_TSW = eps_new*((-ay + az - bx + bz + cx + cy-kxyAverage+kxzAverage+kyzAverage)/(72.*o));
+ f_TSE = eps_new*((ay - az + bx + bz - cx + cy+kxyAverage-kxzAverage+kyzAverage)/(72.*o));
+ f_TNW = eps_new*((ay + az + bx - bz + cx - cy+kxyAverage+kxzAverage-kyzAverage)/(72.*o));
+
+ x_E = 0.25*eps_new*((2*(-4*axx + bxy + cxz))/(27.*o));
+ x_N = 0.25*eps_new*((2*(2*axx - 2*bxy + cxz))/(27.*o));
+ x_T = 0.25*eps_new*((2*(2*axx + bxy - 2*cxz))/(27.*o));
+ x_NE = 0.25*eps_new*(-((2*axx + 3*axy + 6*bxx + bxy - 2*cxz))/(54.*o));
+ x_SE = 0.25*eps_new*(-((2*axx - 3*axy - 6*bxx + bxy - 2*cxz))/(54.*o));
+ x_TE = 0.25*eps_new*(-((2*axx + 3*axz - 2*bxy + 6*cxx + cxz))/(54.*o));
+ x_BE = 0.25*eps_new*(-((2*axx - 3*axz - 2*bxy - 6*cxx + cxz))/(54.*o));
+ x_TN = 0.25*eps_new*(-((-4*axx + bxy + 3*bxz + 3*cxy + cxz))/(54.*o));
+ x_BN = 0.25*eps_new*(-((-4*axx + bxy - 3*bxz - 3*cxy + cxz))/(54.*o));
+ x_ZERO = 0.;
+ x_TNE = 0.25*eps_new*(-((axy + axz + 2*bxx + bxz + 2*cxx + cxy))/(72.*o));
+ x_TSW = 0.25*eps_new*(((-axy + axz - 2*bxx + bxz + 2*cxx + cxy))/(72.*o));
+ x_TSE = 0.25*eps_new*(((axy - axz + 2*bxx + bxz - 2*cxx + cxy))/(72.*o));
+ x_TNW = 0.25*eps_new*(((axy + axz + 2*bxx - bxz + 2*cxx - cxy))/(72.*o));
+
+ y_E = 0.25*eps_new*(2*(-2*axy + 2*byy + cyz))/(27.*o);
+ y_N = 0.25*eps_new*(2*(axy - 4*byy + cyz))/(27.*o);
+ y_T = 0.25*eps_new*(2*(axy + 2*byy - 2*cyz))/(27.*o);
+ y_NE = 0.25*eps_new*(-((axy + 6*ayy + 3*bxy + 2*byy - 2*cyz))/(54.*o));
+ y_SE = 0.25*eps_new*(-((axy - 6*ayy - 3*bxy + 2*byy - 2*cyz))/(54.*o));
+ y_TE = 0.25*eps_new*(-((axy + 3*ayz - 4*byy + 3*cxy + cyz))/(54.*o));
+ y_BE = 0.25*eps_new*(-((axy - 3*ayz - 4*byy - 3*cxy + cyz))/(54.*o));
+ y_TN = 0.25*eps_new*(-((-2*axy + 2*byy + 3*byz + 6*cyy + cyz))/(54.*o));
+ y_BN = 0.25*eps_new*(-((-2*axy + 2*byy - 3*byz - 6*cyy + cyz))/(54.*o));
+ y_ZERO = 0.;
+ y_TNE = 0.25*eps_new*(-((2*ayy + ayz + bxy + byz + cxy + 2*cyy))/(72.*o));
+ y_TSW = 0.25*eps_new*(((-2*ayy + ayz - bxy + byz + cxy + 2*cyy))/(72.*o));
+ y_TSE = 0.25*eps_new*(((2*ayy - ayz + bxy + byz - cxy + 2*cyy))/(72.*o));
+ y_TNW = 0.25*eps_new*(((2*ayy + ayz + bxy - byz + cxy - 2*cyy))/(72.*o));
+
+ z_E = 0.25*eps_new*((2*(-2*axz + byz + 2*czz))/(27.*o));
+ z_N = 0.25*eps_new*((2*(axz - 2*byz + 2*czz))/(27.*o));
+ z_T = 0.25*eps_new*((2*(axz + byz - 4*czz))/(27.*o));
+ z_NE = 0.25*eps_new*(-((axz + 3*ayz + 3*bxz + byz - 4*czz))/(54.*o));
+ z_SE = 0.25*eps_new*(-((axz - 3*ayz - 3*bxz + byz - 4*czz))/(54.*o));
+ z_TE = 0.25*eps_new*(-((axz + 6*azz - 2*byz + 3*cxz + 2*czz))/(54.*o));
+ z_BE = 0.25*eps_new*(-((axz - 6*azz - 2*byz - 3*cxz + 2*czz))/(54.*o));
+ z_TN = 0.25*eps_new*(-((-2*axz + byz + 6*bzz + 3*cyz + 2*czz))/(54.*o));
+ z_BN = 0.25*eps_new*(-((-2*axz + byz - 6*bzz - 3*cyz + 2*czz))/(54.*o));
+ z_ZERO = 0.;
+ z_TNE = 0.25*eps_new*(-((ayz + 2*azz + bxz + 2*bzz + cxz + cyz))/(72.*o));
+ z_TSW = 0.25*eps_new*(((-ayz + 2*azz - bxz + 2*bzz + cxz + cyz))/(72.*o));
+ z_TSE = 0.25*eps_new*(((ayz - 2*azz + bxz + 2*bzz - cxz + cyz))/(72.*o));
+ z_TNW = 0.25*eps_new*(((ayz + 2*azz + bxz - 2*bzz + cxz - cyz))/(72.*o));
+
+ xy_E = 0.0625*eps_new *(( 2.*cxyz)/(27.*o));
+ xy_N = 0.0625*eps_new *(( 2.*cxyz)/(27.*o));
+ xy_T = -(0.0625*eps_new *(( 4.*cxyz)/(27.*o)));
+ xy_NE = 0.0625*eps_new *( cxyz /(27.*o));
+ xy_SE = 0.0625*eps_new *( cxyz /(27.*o));
+ xy_TE = -(0.0625*eps_new *(( 3.*axyz + cxyz)/(54.*o)));
+ xy_BE = -(0.0625*eps_new *((-3.*axyz + cxyz)/(54.*o)));
+ xy_TN = -(0.0625*eps_new *(( 3.*bxyz + cxyz)/(54.*o)));
+ xy_BN = -(0.0625*eps_new *(( - 3.*bxyz + cxyz)/(54.*o)));
+ //xy_ZERO= 0.0625*eps_new;
+ xy_TNE = -(0.0625*eps_new *(( axyz + bxyz )/(72.*o)));
+ xy_TSW = 0.0625*eps_new *(( axyz + bxyz )/(72.*o));
+ xy_TSE = 0.0625*eps_new *((- axyz + bxyz )/(72.*o));
+ xy_TNW = 0.0625*eps_new *(( axyz - bxyz )/(72.*o));
+
+ xz_E = 0.0625*eps_new *(( 2.*bxyz )/(27.*o));
+ xz_N = -(0.0625*eps_new *(( 4.*bxyz )/(27.*o)));
+ xz_T = 0.0625*eps_new *(( 2.*bxyz )/(27.*o));
+ xz_NE = -(0.0625*eps_new *(( 3.*axyz + bxyz )/(54.*o)));
+ xz_SE = -(0.0625*eps_new *((-3.*axyz + bxyz )/(54.*o)));
+ xz_TE = 0.0625*eps_new *(( bxyz )/(27.*o));
+ xz_BE = 0.0625*eps_new *(( bxyz )/(27.*o));
+ xz_TN = -(0.0625*eps_new *(( bxyz + 3.*cxyz)/(54.*o)));
+ xz_BN = -(0.0625*eps_new *(( bxyz - 3.*cxyz)/(54.*o)));
+ //xz_ZERO= 0.0625*eps_new;
+ xz_TNE = -(0.0625*eps_new *(( axyz + cxyz)/(72.*o)));
+ xz_TSW = 0.0625*eps_new *((- axyz + cxyz)/(72.*o));
+ xz_TSE = 0.0625*eps_new *(( axyz + cxyz)/(72.*o));
+ xz_TNW = 0.0625*eps_new *(( axyz - cxyz)/(72.*o));
+
+ yz_E = -(0.0625*eps_new *(( 4.*axyz )/(27.*o)));
+ yz_N = 0.0625*eps_new *(( 2.*axyz )/(27.*o));
+ yz_T = 0.0625*eps_new *(( 2.*axyz )/(27.*o));
+ yz_NE = -(0.0625*eps_new *(( axyz + 3.*bxyz )/(54.*o)));
+ yz_SE = -(0.0625*eps_new *(( axyz - 3.*bxyz )/(54.*o)));
+ yz_TE = -(0.0625*eps_new *(( axyz + 3.*cxyz)/(54.*o)));
+ yz_BE = -(0.0625*eps_new *(( axyz - 3.*cxyz)/(54.*o)));
+ yz_TN = 0.0625*eps_new *(( axyz )/(27.*o));
+ yz_BN = 0.0625*eps_new *(( axyz )/(27.*o));
+ //yz_ZERO= 0.0625*eps_new;
+ yz_TNE = -(0.0625*eps_new *(( bxyz + cxyz)/(72.*o)));
+ yz_TSW = 0.0625*eps_new *(( - bxyz + cxyz)/(72.*o));
+ yz_TSE = 0.0625*eps_new *(( bxyz - cxyz)/(72.*o));
+ yz_TNW = 0.0625*eps_new *(( bxyz + cxyz)/(72.*o));
+}
+
//////////////////////////////////////////////////////////////////////////
void CompressibleOffsetMomentsInterpolationProcessor::calcInterpolatedNodeCF(real* f, real omega, real x, real y, real z, real press, real xs, real ys, real zs)
{
diff --git a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
index 174256f0250a01b01a6d1072f94962a2d3122991..4d684b0ddd6bc705e22bfc9aac6be940ebbfafdb 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
+++ b/src/cpu/VirtualFluidsCore/LBM/CompressibleOffsetMomentsInterpolationProcessor.h
@@ -61,7 +61,9 @@ private:
real calcPressBNE();
void calcInterpolatedVelocity(real x, real y, real z,real& vx1, real& vx2, real& vx3) override;
void calcInterpolatedShearStress(real x, real y, real z,real& tauxx, real& tauyy, real& tauzz,real& tauxy, real& tauxz, real& tauyz) override;
-
+ void calcMoments(const real* const f, real omega, real& rho, real& vx1, real& vx2, real& vx3,
+ real& kxy, real& kyz, real& kxz, real& kxxMyy, real& kxxMzz);
+ void calcInterpolatedCoefficiets(const D3Q27ICell& icell, real omega, real eps_new) override;
void calcInterpolation(const D3Q27ICell& icell, real* icellC);
};
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
index 23dcb1554f49dfe53f52dc7624bcd173ca22cdae..21d126496c4671ae87bac8988487ae578a4f28d3 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
@@ -167,10 +167,15 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
+
readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_MMP, vx1_MMP, vx2_MMP, vx3_MMP,
kxyFromfcNEQ_MMP, kyzFromfcNEQ_MMP, kxzFromfcNEQ_MMP, kxxMyyFromfcNEQ_MMP, kxxMzzFromfcNEQ_MMP);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_MMP, vx1_MMP, vx2_MMP, vx3_MMP,
+ // kxyFromfcNEQ_MMP, kyzFromfcNEQ_MMP, kxzFromfcNEQ_MMP, kxxMyyFromfcNEQ_MMP, kxxMzzFromfcNEQ_MMP);
+
//////////////////////////////////////////////////////////////////////////
// source node TSE = PMP
//////////////////////////////////////////////////////////////////////////
@@ -186,10 +191,14 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PMP, vx1_PMP, vx2_PMP, vx3_PMP,
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_PMP, vx1_PMP, vx2_PMP, vx3_PMP,
kxyFromfcNEQ_PMP, kyzFromfcNEQ_PMP, kxzFromfcNEQ_PMP, kxxMyyFromfcNEQ_PMP, kxxMzzFromfcNEQ_PMP);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PMP, vx1_PMP, vx2_PMP, vx3_PMP,
+ // kxyFromfcNEQ_PMP, kyzFromfcNEQ_PMP, kxzFromfcNEQ_PMP, kxxMyyFromfcNEQ_PMP, kxxMzzFromfcNEQ_PMP);
+
//////////////////////////////////////////////////////////////////////////
// source node BSE = PMM
//////////////////////////////////////////////////////////////////////////
@@ -205,10 +214,14 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PMM, vx1_PMM, vx2_PMM, vx3_PMM,
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_PMM, vx1_PMM, vx2_PMM, vx3_PMM,
kxyFromfcNEQ_PMM, kyzFromfcNEQ_PMM, kxzFromfcNEQ_PMM, kxxMyyFromfcNEQ_PMM, kxxMzzFromfcNEQ_PMM);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PMM, vx1_PMM, vx2_PMM, vx3_PMM,
+ // kxyFromfcNEQ_PMM, kyzFromfcNEQ_PMM, kxzFromfcNEQ_PMM, kxxMyyFromfcNEQ_PMM, kxxMzzFromfcNEQ_PMM);
+
//////////////////////////////////////////////////////////////////////////
// source node BNW = MPM
//////////////////////////////////////////////////////////////////////////
@@ -234,10 +247,14 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_MPM, vx1_MPM, vx2_MPM, vx3_MPM,
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_MPM, vx1_MPM, vx2_MPM, vx3_MPM,
kxyFromfcNEQ_MPM, kyzFromfcNEQ_MPM, kxzFromfcNEQ_MPM, kxxMyyFromfcNEQ_MPM, kxxMzzFromfcNEQ_MPM);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_MPM, vx1_MPM, vx2_MPM, vx3_MPM,
+ // kxyFromfcNEQ_MPM, kyzFromfcNEQ_MPM, kxzFromfcNEQ_MPM, kxxMyyFromfcNEQ_MPM, kxxMzzFromfcNEQ_MPM);
+
//////////////////////////////////////////////////////////////////////////
// source node TNW = MPP
//////////////////////////////////////////////////////////////////////////
@@ -252,10 +269,14 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
k_MMM = neighborZfine[k_MMM];
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
-
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_MPP, vx1_MPP, vx2_MPP, vx3_MPP,
+
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_MPP, vx1_MPP, vx2_MPP, vx3_MPP,
kxyFromfcNEQ_MPP, kyzFromfcNEQ_MPP, kxzFromfcNEQ_MPP, kxxMyyFromfcNEQ_MPP, kxxMzzFromfcNEQ_MPP);
+
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_MPP, vx1_MPP, vx2_MPP, vx3_MPP,
+ // kxyFromfcNEQ_MPP, kyzFromfcNEQ_MPP, kxzFromfcNEQ_MPP, kxxMyyFromfcNEQ_MPP, kxxMzzFromfcNEQ_MPP);
//////////////////////////////////////////////////////////////////////////
// source node TNE = PPP
@@ -272,10 +293,15 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PPP, vx1_PPP, vx2_PPP, vx3_PPP,
+
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_PPP, vx1_PPP, vx2_PPP, vx3_PPP,
kxyFromfcNEQ_PPP, kyzFromfcNEQ_PPP, kxzFromfcNEQ_PPP, kxxMyyFromfcNEQ_PPP, kxxMzzFromfcNEQ_PPP);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PPP, vx1_PPP, vx2_PPP, vx3_PPP,
+ // kxyFromfcNEQ_PPP, kyzFromfcNEQ_PPP, kxzFromfcNEQ_PPP, kxxMyyFromfcNEQ_PPP, kxxMzzFromfcNEQ_PPP);
+
//////////////////////////////////////////////////////////////////////////
// source node BNE = PPM
//////////////////////////////////////////////////////////////////////////
@@ -291,10 +317,15 @@ template<bool hasTurbulentViscosity> __global__ void scaleFC_compressible(
if(hasTurbulentViscosity) omegaF = omegaFine/ (c1o1 + c3o1*omegaFine*turbulentViscosityFine[k_000]);
- calculateMomentsOnSourceNodes( distFine, omegaF,
- k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PPM, vx1_PPM, vx2_PPM, vx3_PPM,
+
+ readDistributionFromList(distribution, distFine, k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM);
+ vf::lbm::calculateMomentsOnSourceNodes(distribution.f, omegaF, drho_PPM, vx1_PPM, vx2_PPM, vx3_PPM,
kxyFromfcNEQ_PPM, kyzFromfcNEQ_PPM, kxzFromfcNEQ_PPM, kxxMyyFromfcNEQ_PPM, kxxMzzFromfcNEQ_PPM);
+ // calculateMomentsOnSourceNodes( distFine, omegaF,
+ // k_000, k_M00, k_0M0, k_00M, k_MM0, k_M0M, k_0MM, k_MMM, drho_PPM, vx1_PPM, vx2_PPM, vx3_PPM,
+ // kxyFromfcNEQ_PPM, kyzFromfcNEQ_PPM, kxzFromfcNEQ_PPM, kxxMyyFromfcNEQ_PPM, kxxMzzFromfcNEQ_PPM);
+
////////////////////////////////////////////////////////////////////////////////
//! - Set the relative position of the offset cell {-1, 0, 1}
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
index 72a2b2bc19bfd420897db6292539b2ff2b71ec50..ca40fb3038ea589d8c201cf3cd2100b1c547af67 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
@@ -51,33 +51,33 @@ __device__ __inline__ void readDistributionFromList(vf::lbm::Distribution27 &dis
unsigned int &k_MM0, unsigned int &k_M0M, unsigned int &k_0MM,
unsigned int &k_MMM)
{
- distribution.f[vf::lbm::dir::PZZ] = (dist.f[DIR_000])[k_000];
- distribution.f[vf::lbm::dir::MZZ] = (dist.f[DIR_P00])[k_000];
- distribution.f[vf::lbm::dir::ZPZ] = (dist.f[DIR_M00])[k_M00];
- distribution.f[vf::lbm::dir::ZMZ] = (dist.f[DIR_0P0])[k_000];
- distribution.f[vf::lbm::dir::ZZP] = (dist.f[DIR_0M0])[k_0M0];
- distribution.f[vf::lbm::dir::ZZM] = (dist.f[DIR_00P])[k_000];
- distribution.f[vf::lbm::dir::PPZ] = (dist.f[DIR_00M])[k_00M];
- distribution.f[vf::lbm::dir::MMZ] = (dist.f[DIR_PP0])[k_000];
- distribution.f[vf::lbm::dir::PMZ] = (dist.f[DIR_MM0])[k_MM0];
- distribution.f[vf::lbm::dir::MPZ] = (dist.f[DIR_PM0])[k_0M0];
- distribution.f[vf::lbm::dir::PZP] = (dist.f[DIR_MP0])[k_M00];
- distribution.f[vf::lbm::dir::MZM] = (dist.f[DIR_P0P])[k_000];
- distribution.f[vf::lbm::dir::PZM] = (dist.f[DIR_M0M])[k_M0M];
- distribution.f[vf::lbm::dir::MZP] = (dist.f[DIR_P0M])[k_00M];
- distribution.f[vf::lbm::dir::ZPP] = (dist.f[DIR_M0P])[k_M00];
- distribution.f[vf::lbm::dir::ZMM] = (dist.f[DIR_0PP])[k_000];
- distribution.f[vf::lbm::dir::ZPM] = (dist.f[DIR_0MM])[k_0MM];
- distribution.f[vf::lbm::dir::ZMP] = (dist.f[DIR_0PM])[k_00M];
- distribution.f[vf::lbm::dir::PPP] = (dist.f[DIR_0MP])[k_0M0];
- distribution.f[vf::lbm::dir::MPP] = (dist.f[DIR_PPP])[k_000];
- distribution.f[vf::lbm::dir::PMP] = (dist.f[DIR_MPP])[k_M00];
- distribution.f[vf::lbm::dir::MMP] = (dist.f[DIR_PMP])[k_0M0];
- distribution.f[vf::lbm::dir::PPM] = (dist.f[DIR_MMP])[k_MM0];
- distribution.f[vf::lbm::dir::MPM] = (dist.f[DIR_PPM])[k_00M];
- distribution.f[vf::lbm::dir::PMM] = (dist.f[DIR_MPM])[k_M0M];
- distribution.f[vf::lbm::dir::MMM] = (dist.f[DIR_PMM])[k_0MM];
- distribution.f[vf::lbm::dir::ZZZ] = (dist.f[DIR_MMM])[k_MMM];
+ distribution.f[DIR_000] = (dist.f[DIR_000])[k_000];
+ distribution.f[DIR_P00] = (dist.f[DIR_P00])[k_000];
+ distribution.f[DIR_M00] = (dist.f[DIR_M00])[k_M00];
+ distribution.f[DIR_0P0] = (dist.f[DIR_0P0])[k_000];
+ distribution.f[DIR_0M0] = (dist.f[DIR_0M0])[k_0M0];
+ distribution.f[DIR_00P] = (dist.f[DIR_00P])[k_000];
+ distribution.f[DIR_00M] = (dist.f[DIR_00M])[k_00M];
+ distribution.f[DIR_PP0] = (dist.f[DIR_PP0])[k_000];
+ distribution.f[DIR_MM0] = (dist.f[DIR_MM0])[k_MM0];
+ distribution.f[DIR_PM0] = (dist.f[DIR_PM0])[k_0M0];
+ distribution.f[DIR_MP0] = (dist.f[DIR_MP0])[k_M00];
+ distribution.f[DIR_P0P] = (dist.f[DIR_P0P])[k_000];
+ distribution.f[DIR_M0M] = (dist.f[DIR_M0M])[k_M0M];
+ distribution.f[DIR_P0M] = (dist.f[DIR_P0M])[k_00M];
+ distribution.f[DIR_M0P] = (dist.f[DIR_M0P])[k_M00];
+ distribution.f[DIR_0PP] = (dist.f[DIR_0PP])[k_000];
+ distribution.f[DIR_0MM] = (dist.f[DIR_0MM])[k_0MM];
+ distribution.f[DIR_0PM] = (dist.f[DIR_0PM])[k_00M];
+ distribution.f[DIR_0MP] = (dist.f[DIR_0MP])[k_0M0];
+ distribution.f[DIR_PPP] = (dist.f[DIR_PPP])[k_000];
+ distribution.f[DIR_MPP] = (dist.f[DIR_MPP])[k_M00];
+ distribution.f[DIR_PMP] = (dist.f[DIR_PMP])[k_0M0];
+ distribution.f[DIR_MMP] = (dist.f[DIR_MMP])[k_MM0];
+ distribution.f[DIR_PPM] = (dist.f[DIR_PPM])[k_00M];
+ distribution.f[DIR_MPM] = (dist.f[DIR_MPM])[k_M0M];
+ distribution.f[DIR_PMM] = (dist.f[DIR_PMM])[k_0MM];
+ distribution.f[DIR_MMM] = (dist.f[DIR_MMM])[k_MMM];
}