diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 0b92f834e538f8d027b40adf7bbfffb469164951..1e444b79b734285653c6f805bbbc33018b394365 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -2,7 +2,7 @@
#add_subdirectory(Applications/gridRf)
#add_subdirectory(Applications/greenvortex)
# add_subdirectory(Applications/micropart)
-#add_subdirectory(Applications/sphere)
+add_subdirectory(${APPS_ROOT_CPU}/sphere)
#add_subdirectory(Applications/vfscript)
#add_subdirectory(Applications/reefer)
#add_subdirectory(Applications/bananas)
diff --git a/apps/cpu/HerschelBulkleySphere/CMakeLists.txt b/apps/cpu/HerschelBulkleySphere/CMakeLists.txt
index 162bf56047faa176ced95b7b7ad780c40f666be7..1f0641fa5992f811dd3f6f7ce5e0a827c46aaa24 100644
--- a/apps/cpu/HerschelBulkleySphere/CMakeLists.txt
+++ b/apps/cpu/HerschelBulkleySphere/CMakeLists.txt
@@ -7,21 +7,4 @@ PROJECT(HerschelBulkleySphere)
INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
-################################################################
-## LOCAL FILES ###
-################################################################
-# FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h
- # ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp
- # ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp )
-
-# SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES})
-# SOURCE_GROUP(src FILES ${SPECIFIC_FILES})
-
-# SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids)
-
-################################################################
-## CREATE PROJECT ###
-################################################################
-# CREATE_CAB_PROJECT(hbsphere BINARY)
-
vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES hbsphere.cpp )
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cfg b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
index 1f636ec212b5cd8500777006e01eb2f7795e2a27..ebcc5d1f8cdffa35905949c64297750994b587ea 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
@@ -1,24 +1,27 @@
-outputPath = /work/koskuche/Herschel-BulkleySphere
+outputPath = d:/temp/Herschel-BulkleySphere
numOfThreads = 1
availMem = 8e9
logToFile = false
-blocknx = 10 10 10
-boundingBox = 600 600 600 #30*20=600**3=216000000
+blocknx = 5 5 5
+boundingBox = 40 40 40 #30*20=600**3=216000000
deltax = 1
-radius = 15
+refineLevel = 1
+
+radius = 5
velocity = 1e-3
n = 0.3
Re = 1
Bn = 0.01
+
newStart = true
restartStep = 100000
cpStart = 10000
cpStep = 10000
-outTime = 10000
-endTime = 1000000
\ No newline at end of file
+outTime = 1
+endTime = 1000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
index deb44fc212b7ad202d4d1ba0a03a9f2c3e2de244..b437a2736072ecf9c7242c9f10e688c37a6fcd31 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
@@ -21,7 +21,7 @@ void bflow(string configname)
double endTime = config.getValue<double>("endTime");
double outTime = config.getValue<double>("outTime");
double availMem = config.getValue<double>("availMem");
- //int refineLevel = config.getValue<int>("refineLevel");
+ int refineLevel = config.getValue<int>("refineLevel");
bool logToFile = config.getValue<bool>("logToFile");
double restartStep = config.getValue<double>("restartStep");
double deltax = config.getValue<double>("deltax");
@@ -164,6 +164,19 @@ void bflow(string configname)
GenBlocksGridVisitor genBlocks(gridCube);
grid->accept(genBlocks);
+ if (refineLevel > 0)
+ {
+ GbCuboid3DPtr refCube(new GbCuboid3D(-10, -10, -10, 0, 0, 0));
+ if (myid == 0) GbSystem3D::writeGeoObject(refCube.get(), outputPath + "/geo/refCube", WbWriterVtkXmlASCII::getInstance());
+
+ if (myid == 0) UBLOG(logINFO, "Refinement - start");
+ RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel, comm);
+ //refineHelper.addGbObject(sphere, refineLevel);
+ refineHelper.addGbObject(refCube, refineLevel);
+ refineHelper.refine();
+ if (myid == 0) UBLOG(logINFO, "Refinement - end");
+ }
+
//walls
GbCuboid3DPtr wallZmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_minX3));
if (myid == 0) GbSystem3D::writeGeoObject(wallZmin.get(), outputPath + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
@@ -206,7 +219,7 @@ void bflow(string configname)
intHelper.addInteractor(wallYmaxInt);
intHelper.addInteractor(wallXminInt);
intHelper.addInteractor(wallXmaxInt);
- intHelper.addInteractor(sphereInt);
+ //intHelper.addInteractor(sphereInt);
intHelper.selectBlocks();
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
//////////////////////////////////////
@@ -242,6 +255,12 @@ void bflow(string configname)
SetKernelBlockVisitor kernelVisitor(kernel, k, availMem, needMem);
grid->accept(kernelVisitor);
+ if (refineLevel > 0)
+ {
+ SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ grid->accept(undefNodesVisitor);
+ }
+
//BC
intHelper.setBC();
@@ -261,14 +280,14 @@ void bflow(string configname)
omp_set_num_threads(numOfThreads);
//set connectors
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ InterpolationProcessorPtr iProcessor(new ThixotropyInterpolationProcessor());
SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, k, iProcessor);
grid->accept(setConnsVisitor);
grid->accept(bcVisitor);
SPtr<UbScheduler> geoSch(new UbScheduler(1));
- WriteBoundaryConditionsCoProcessor ppgeo = WriteBoundaryConditionsCoProcessor(grid, geoSch, outputPath, WbWriterVtkXmlBinary::getInstance(), comm);
+ WriteBoundaryConditionsCoProcessor ppgeo = WriteBoundaryConditionsCoProcessor(grid, geoSch, outputPath, WbWriterVtkXmlASCII::getInstance(), comm);
ppgeo.process(0);
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
diff --git a/apps/cpu/sphere/CMakeLists.txt b/apps/cpu/sphere/CMakeLists.txt
index 77d7e0a41ea36fef0813ace18ee699da481d2ddd..022df23ea789b9a5c62be699c9b3119ef33d8a4b 100644
--- a/apps/cpu/sphere/CMakeLists.txt
+++ b/apps/cpu/sphere/CMakeLists.txt
@@ -5,21 +5,6 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
########################################################
PROJECT(sphere)
-INCLUDE(${APPS_ROOT}/IncludsList.cmake)
+INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
-#################################################################
-### LOCAL FILES ###
-#################################################################
-FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h
- ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp
- ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp )
-
-SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES})
-SOURCE_GROUP(src FILES ${SPECIFIC_FILES})
-
-SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids)
-
-#################################################################
-### CREATE PROJECT ###
-#################################################################
-CREATE_CAB_PROJECT(sphere BINARY)
+vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES sphere.cpp )
\ No newline at end of file
diff --git a/apps/cpu/sphere/config.txt b/apps/cpu/sphere/config.txt
index cbd03fd4bb66449608a1ea07a79295b3bf340acb..aedb74bf054561257683e86c8991c6a3d8800d7a 100644
--- a/apps/cpu/sphere/config.txt
+++ b/apps/cpu/sphere/config.txt
@@ -1,7 +1,7 @@
-#Ordner für Simulationsergebnisse
+#Ordner f�r Simulationsergebnisse
path=d:/temp/sphere
-#Verfügbare Arbeitsspeicher in Byte
+#Verf�gbare Arbeitsspeicher in Byte
memory=3e9
#Pfad zum Metafile
@@ -11,12 +11,12 @@ metafile=d:/Data/insituDemo/metafile.csv
outstep=1
#maximale Anzahl Simulationszeitschritte
-endstep=100000
+endstep=10
#Anzahl von Threads
threads=1
#max refierment level (1 - 5)
-level=4
+level=1
test = true
\ No newline at end of file
diff --git a/apps/cpu/sphere/sphere.cpp b/apps/cpu/sphere/sphere.cpp
index 872e7b81f402f5288eaf8b0f56678d7dfb4c34d9..c33a3d696474958961bee538ae8b9417d6f38d0f 100644
--- a/apps/cpu/sphere/sphere.cpp
+++ b/apps/cpu/sphere/sphere.cpp
@@ -20,21 +20,7 @@ void run(string configname)
int mybundle = comm->getBundleID();
int root = comm->getRoot();
- //ConfigFileReader cf(cstr);
- //if ( !cf.read() )
- //{
- // std::string exceptionText = "Unable to read configuration file\n";
- // throw exceptionText;
- //}
-
- //pathname = cf.getValue("path");
- //availMem = UbSystem::stringTo<double>(cf.getValue("memory"));
- //string metafile = cf.getValue("metafile");
- //double outstep = UbSystem::stringTo<double>(cf.getValue("outstep"));
- //double endstep = UbSystem::stringTo<double>(cf.getValue("endstep"));
- //int numOfThreads = UbSystem::stringTo<int>(cf.getValue("threads"));
-
- ConfigurationFile config;
+ ConfigurationFile config;
config.load(configname);
string pathname = config.getValue<string>("path");
@@ -47,8 +33,8 @@ void run(string configname)
bool test = config.getValue<bool>("test");
- LBMReal radius = 4;
- LBMReal uLB = 0.1;
+ LBMReal radius = 10;
+ LBMReal uLB = 0.01;
LBMReal Re = 1;
LBMReal rhoLB = 0.0;
//LBMReal nuLB = (uLB*2.0*radius)/Re;
@@ -79,13 +65,13 @@ void run(string configname)
double dx = 1;
- const int blocknx1 = 8;
- const int blocknx2 = 8;
- const int blocknx3 = 8;
+ const int blocknx1 = 5;
+ const int blocknx2 = 5;
+ const int blocknx3 = 5;
- const int gridNx1 = 4;//18;
- const int gridNx2 = 4;// 11;
- const int gridNx3 = 4;// 11;
+ const int gridNx1 = 8;//18;
+ const int gridNx2 = 8;// 11;
+ const int gridNx3 = 8;// 11;
//const int blocknx1 = 40;
//const int blocknx2 = 40;
@@ -106,11 +92,6 @@ void run(string configname)
grid->setDeltaX(dx);
grid->setBlockNX(blocknx1, blocknx2, blocknx3);
- //////////////////////////////////////////////////////////////////////////
- //restart
- SPtr<UbScheduler> restartSch(new UbScheduler(100000, 100000, 100000));
- RestartCoProcessor rp(grid, restartSch, comm, pathname, RestartCoProcessor::BINARY);
- //////////////////////////////////////////////////////////////////////////
if (grid->getTimeStep() == 0)
{
@@ -153,7 +134,7 @@ void run(string configname)
//sphere
- SPtr<GbObject3D> sphere(new GbSphere3D(L1 / 2.0, L2*0.5, L3*0.5, radius));
+ SPtr<GbObject3D> sphere(new GbSphere3D(L1*0.5, L2*0.5, L3*0.5, radius));
//SPtr<GbObject3D> sphere(new GbSphere3D(L1/2.0-4.0, L2*0.5+4.0, L3*0.5+4.0, radius));
//SPtr<GbObject3D> sphere(new GbCuboid3D(L1/4.0-radius, L2/2.0-radius, L3/2.0-radius, L1/4.0+radius, L2/2.0+radius, L3/2.0+radius));
GbSystem3D::writeGeoObject(sphere.get(), pathname + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
@@ -194,7 +175,7 @@ void run(string configname)
GbCuboid3DPtr geoOutflow(new GbCuboid3D(d_maxX1, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength));
if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
- WriteBlocksSPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
@@ -232,18 +213,13 @@ void run(string configname)
intHelper.addInteractor(outflowInt);
intHelper.selectBlocks();
- //domain decomposition for threads
- PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
- grid->accept(pqPartVisitor);
+ ////domain decomposition for threads
+ //PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
+ //grid->accept(pqPartVisitor);
-
- //set connectors
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- grid->accept(setConnsVisitor);
-
- //Block3DSPtr<ConnectorFactory> factory(new Block3DConnectorFactory());
- //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory);
+ ////set connectors
+ //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ //SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
//grid->accept(setConnsVisitor);
ppblocks->process(0);
@@ -265,7 +241,8 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
}
- SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel(blocknx1, blocknx2, blocknx3, IncompressibleCumulantLBMKernel::NORMAL));
+ //SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel());
+ SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel());
SPtr<BCProcessor> bcProcessor(new BCProcessor());
@@ -292,17 +269,17 @@ void run(string configname)
fctRoh.DefineConst("l", d_maxX1 - d_minX1);
//initialization of distributions
- InitDistributionsBlockVisitor initVisitor(nuLB, rhoLB);
- initVisitor.setVx1(fct);
+ InitDistributionsBlockVisitor initVisitor;
+ //initVisitor.setVx1(fct);
//initVisitor.setRho(fctRoh);
grid->accept(initVisitor);
//Postrozess
SPtr<UbScheduler> geoSch(new UbScheduler(1));
- WriteBoundaryConditionsSPtr<CoProcessor> ppgeo(
- new WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ SPtr<CoProcessor> ppgeo(
+ new WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
ppgeo->process(0);
- ppgeo.reset();;
+ ppgeo.reset();
if (myid == 0) UBLOG(logINFO, "Preprozess - end");
}
@@ -311,10 +288,11 @@ void run(string configname)
UBLOG(logINFO, "SetConnectors - start, id=" << myid);
//set connectors
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- //D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- SPtr<ConnectorFactory> cFactory(new Block3DConnectorFactory());
- ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, cFactory);
+ //SPtr<InterpolationProcessor> iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ SPtr<CompressibleOffsetMomentsInterpolationProcessor> iProcessor(new CompressibleOffsetMomentsInterpolationProcessor());
+ SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
+ //SPtr<ConnectorFactory> cFactory(new Block3DConnectorFactory());
+ //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, cFactory);
grid->accept(setConnsVisitor);
UBLOG(logINFO, "SetConnectors - end, id=" << myid);
@@ -322,21 +300,20 @@ void run(string configname)
SPtr<UbScheduler> stepSch(new UbScheduler(outstep));
//stepSch->addSchedule(10000, 0, 1000000);
- WriteMacroscopicQuantitiesCoProcessor pp(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv,comm);
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv,comm));
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
- NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
-
- const SPtr<ConcreteCalculatorFactory> calculatorFactory = std::make_shared<ConcreteCalculatorFactory>(stepSch);
- CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endstep, calculatorFactory, CalculatorType::HYBRID));
+ SPtr<CoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
- if (myid == 0)
- UBLOG(logINFO, "Simulation-start");
+ SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
+ SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endstep));
+ calculator->addCoProcessor(npr);
+ calculator->addCoProcessor(writeMQCoProcessor);
- calculation->calculate();
- if (myid == 0)
- UBLOG(logINFO, "Simulation-end");
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculator->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
}
catch (std::exception& e)
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index 31671c3ed7fc7e5c619753d21313b2af5178288d..31937b7ef94da983263a73e86f65b3c22aeaf33b 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -221,6 +221,7 @@
#include <LBM/ThixotropyModelLBMKernel.h>
#include <LBM/BinghamModelLBMKernel.h>
#include <LBM/HerschelBulkleyModelLBMKernel.h>
+#include <LBM/ThixotropyInterpolationProcessor.h>
#include <geometry3d/CoordinateTransformation3D.h>
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
index bae45335e4c98bb60b6f67c65de3ec95ec19a3a0..b7b6246d9e7ceb99c5adb288cc4c01ead504777b 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
@@ -148,8 +148,13 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
SPtr<ILBMKernel> kernel = block->getKernel();
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+<<<<<<< HEAD
// knotennummerierung faengt immer bei 0 an!
unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+=======
+ //knotennummerierung faengt immer bei 0 an!
+ int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+>>>>>>> add grid refinement for thixotropic fluid
int minX1 = 0;
int minX2 = 0;
@@ -209,24 +214,31 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
//}
}
}
- }
- }
-
- maxX1 -= 1;
- maxX2 -= 1;
- maxX3 -= 1;
-
- // cell vector erstellen
- for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
- for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
- for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
- if ((SWB = nodeNumbers(ix1, ix2, ix3)) >= 0 && (SEB = nodeNumbers(ix1 + 1, ix2, ix3)) >= 0 &&
- (NEB = nodeNumbers(ix1 + 1, ix2 + 1, ix3)) >= 0 && (NWB = nodeNumbers(ix1, ix2 + 1, ix3)) >= 0 &&
- (SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
- (NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
- (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
- cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
- }
+ }
+ }
+ }
+
+ maxX1 -= 1;
+ maxX2 -= 1;
+ maxX3 -= 1;
+
+ //cell vector erstellen
+ for (int ix3 = minX3; ix3<=maxX3; ix3++)
+ {
+ for (int ix2 = minX2; ix2<=maxX2; ix2++)
+ {
+ for (int ix1 = minX1; ix1<=maxX1; ix1++)
+ {
+ if ((SWB = nodeNumbers(ix1, ix2, ix3))>=0
+ &&(SEB = nodeNumbers(ix1+1, ix2, ix3))>=0
+ &&(NEB = nodeNumbers(ix1+1, ix2+1, ix3))>=0
+ &&(NWB = nodeNumbers(ix1, ix2+1, ix3))>=0
+ &&(SWT = nodeNumbers(ix1, ix2, ix3+1))>=0
+ &&(SET = nodeNumbers(ix1+1, ix2, ix3+1))>=0
+ &&(NET = nodeNumbers(ix1+1, ix2+1, ix3+1))>=0
+ &&(NWT = nodeNumbers(ix1, ix2+1, ix3+1))>=0)
+ {
+ cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB, (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET, (unsigned int)NET, (unsigned int)NWT));
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
index 1d71e4bbf809153dfec33c33bb3a1fd4a6f468ae..e5976c68276572681ef7d3480fc8f0e2f1311ebf 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
@@ -155,7 +155,7 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
datanames.push_back("Vy");
datanames.push_back("Vz");
//datanames.push_back("Press");
- //datanames.push_back("Level");
+ datanames.push_back("Level");
//datanames.push_back("BlockID");
//datanames.push_back("gamma");
//datanames.push_back("collFactor");
@@ -170,7 +170,7 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
LBMReal vx1,vx2,vx3,rho;
//knotennummerierung faengt immer bei 0 an!
- unsigned int SWB,SEB,NEB,NWB,SWT,SET,NET,NWT;
+ int SWB,SEB,NEB,NWB,SWT,SET,NET,NWT;
if(block->getKernel()->getCompressible())
{
@@ -198,11 +198,13 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
//int maxX3 = (int)(distributions->getNX3());
//nummern vergeben und node vector erstellen + daten sammeln
- CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3,-1);
+ CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
maxX1 -= 2;
maxX2 -= 2;
maxX3 -= 2;
+
+
//D3Q27BoundaryConditionPtr bcPtr;
int nr = (int)nodes.size();
@@ -226,25 +228,20 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
double press = D3Q27System::getPressure(f); //D3Q27System::calcPress(f,rho,vx1,vx2,vx3);
if (UbMath::isNaN(rho) || UbMath::isInfinity(rho))
- UB_THROW( UbException(UB_EXARGS,"rho is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
- ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
- //rho=999.0;
+ //UB_THROW( UbException(UB_EXARGS,"rho is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ rho=999.0;
if (UbMath::isNaN(press) || UbMath::isInfinity(press))
- UB_THROW( UbException(UB_EXARGS,"press is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
- ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
- //press=999.0;
+ //UB_THROW( UbException(UB_EXARGS,"press is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ press=999.0;
if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
- UB_THROW( UbException(UB_EXARGS,"vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
- ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
- //vx1=999.0;
+ //UB_THROW( UbException(UB_EXARGS,"vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx1=999.0;
if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
- UB_THROW( UbException(UB_EXARGS,"vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
- ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
- //vx2=999.0;
+ //UB_THROW( UbException(UB_EXARGS,"vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx2=999.0;
if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
- UB_THROW( UbException(UB_EXARGS,"vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
- ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
- //vx3 = 999.0;
+ //UB_THROW( UbException(UB_EXARGS,"vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx3 = 999.0;
data[index++].push_back(rho);
data[index++].push_back(vx1);
@@ -263,43 +260,32 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
//data[index++].push_back(vx2 * conv->getFactorVelocityLbToW());
//data[index++].push_back(vx3 * conv->getFactorVelocityLbToW());
//data[index++].push_back((press * conv->getFactorPressureLbToW()) / ((rho+1.0) * conv->getFactorDensityLbToW()));
- //data[index++].push_back(level);
+ data[index++].push_back(level);
//data[index++].push_back(blockID);
}
- }
- }
- maxX1 -= 1;
- maxX2 -= 1;
- maxX3 -= 1;
-
- unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
- int SWBi, SEBi, NEBi, NWBi, SWTi, SETi, NETi, NWTi;
- // cell vector erstellen
- for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
- for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
- for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
- if (
- ( SWBi = nodeNumbers(ix1, ix2, ix3)) >= 0
- && (SEBi = nodeNumbers(ix1 + 1, ix2, ix3)) >= 0
- && (NEBi = nodeNumbers(ix1 + 1, ix2 + 1, ix3)) >= 0
- && (NWBi = nodeNumbers(ix1, ix2 + 1, ix3)) >= 0
- && (SWTi = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0
- && (SETi = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0
- && (NETi = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0
- && (NWTi = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0
- )
- {
- SWB =SWBi;
- SEB =SEBi;
- NEB =NEBi;
- NWB =NWBi;
- SWT =SWTi;
- SET =SETi;
- NET =NETi;
- NWT =NWTi;
-
- cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
- }
+ }
+ }
+ }
+ maxX1 -= 1;
+ maxX2 -= 1;
+ maxX3 -= 1;
+ //cell vector erstellen
+ for(int ix3=minX3; ix3<=maxX3; ix3++)
+ {
+ for(int ix2=minX2; ix2<=maxX2; ix2++)
+ {
+ for(int ix1=minX1; ix1<=maxX1; ix1++)
+ {
+ if( (SWB=nodeNumbers( ix1 , ix2, ix3 )) >= 0
+ && (SEB=nodeNumbers( ix1+1, ix2, ix3 )) >= 0
+ && (NEB=nodeNumbers( ix1+1, ix2+1, ix3 )) >= 0
+ && (NWB=nodeNumbers( ix1 , ix2+1, ix3 )) >= 0
+ && (SWT=nodeNumbers( ix1 , ix2, ix3+1 )) >= 0
+ && (SET=nodeNumbers( ix1+1, ix2, ix3+1 )) >= 0
+ && (NET=nodeNumbers( ix1+1, ix2+1, ix3+1 )) >= 0
+ && (NWT=nodeNumbers( ix1 , ix2+1, ix3+1 )) >= 0 )
+ {
+ cells.push_back( makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB, (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET, (unsigned int)NET, (unsigned int)NWT) );
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index bb38862b0fa93e19dda7431a12db9f1afe174d11..235de154c83444f47ab03a4b8ea1ac5926a4d3bb 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -5,42 +5,6 @@
#include <LBMSystem.h>
#include <UbMath.h>
-//namespace Thixotropy
-//{
-// //////////////////////////////////////////////////////////////////////////
-// inline LBMReal getBinghamCollFactor(LBMReal omegaInf, LBMReal yieldSterss, LBMReal shearRate, LBMReal drho)
-// {
-// LBMReal cs2 = one_over_sqrt3 * one_over_sqrt3;
-// LBMReal rho = one + drho;
-// LBMReal omega = omegaInf * (one - (omegaInf * yieldSterss) / (shearRate * cs2 * rho + UbMath::Epsilon<LBMReal>::val()));
-// return omega;
-// }
-// //////////////////////////////////////////////////////////////////////////
-// inline LBMReal getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMReal yieldSterss, LBMReal shearRate, LBMReal drho, LBMReal k, LBMReal n)
-// {
-// LBMReal cs2 = one_over_sqrt3 * one_over_sqrt3;
-// LBMReal rho = one + drho;
-// LBMReal gammaDot = shearRate;
-// LBMReal tau0 = yieldSterss;
-// LBMReal omega = omegaInf;
-// LBMReal epsilon = 1;
-//
-// while (epsilon > 1e-10)
-// {
-// LBMReal omegaOld = omega;
-// LBMReal gammaDotPowN = std::pow(gammaDot, n);
-// LBMReal omegaByOmegaInfPowN = std::pow(omega / omegaInf, n);
-// LBMReal numerator = (2.0 * gammaDotPowN * k * omegaByOmegaInfPowN * omegaInf + cs2 * gammaDot * (omega - 2.0) * rho + 2.0 * omegaInf * tau0);
-// LBMReal denominator = (2.0 * k * n * gammaDotPowN * omegaByOmegaInfPowN * omegaInf + cs2 * gammaDot * rho * omega) + UbMath::Epsilon<LBMReal>::val();
-// omega = omega - omega * numerator / denominator;
-// omega = (omega < zeroReal) ? c1o2 * omegaOld : omega;
-// epsilon = std::abs(omega - omegaOld);
-// }
-//
-// return omega;
-// }
-//}
-
class Thixotropy
{
public:
@@ -58,6 +22,7 @@ public:
static LBMReal getBinghamCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
static LBMReal getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
+ static LBMReal getHerschelBulkleyCollFactorBackward(LBMReal shearRate, LBMReal drho);
private:
Thixotropy();
@@ -99,5 +64,13 @@ inline LBMReal Thixotropy::getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMRea
return omega;
}
+//////////////////////////////////////////////////////////////////////////
+inline LBMReal Thixotropy::getHerschelBulkleyCollFactorBackward(LBMReal shearRate, LBMReal drho)
+{
+ LBMReal rho = UbMath::one + drho;
+ LBMReal gamma = shearRate + UbMath::Epsilon<LBMReal>::val();
+ LBMReal cs2 = UbMath::one_over_sqrt3 * UbMath::one_over_sqrt3;
+ return 1.0 / ((tau0 + k * std::pow(gamma, n)) / (cs2 * rho * gamma) + UbMath::c1o2);
+}
#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a09a1f64d2acd2d71b18556381331740ee63ae5
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
@@ -0,0 +1,666 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ThixotropyInterpolationProcessor.cpp
+//! \ingroup BoundarConditions
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#include "ThixotropyInterpolationProcessor.h"
+#include "D3Q27System.h"
+#include "Thixotropy.h"
+
+
+ThixotropyInterpolationProcessor::ThixotropyInterpolationProcessor()
+ : omegaC(0.0), omegaF(0.0)
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+ThixotropyInterpolationProcessor::ThixotropyInterpolationProcessor(LBMReal omegaC, LBMReal omegaF)
+ : omegaC(omegaC), omegaF(omegaF)
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+ThixotropyInterpolationProcessor::~ThixotropyInterpolationProcessor()
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+InterpolationProcessorPtr ThixotropyInterpolationProcessor::clone()
+{
+ InterpolationProcessorPtr iproc = InterpolationProcessorPtr (new ThixotropyInterpolationProcessor(this->omegaC, this->omegaF));
+ return iproc;
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::setOmegas( LBMReal omegaC, LBMReal omegaF )
+{
+ this->omegaC = omegaC;
+ this->omegaF = omegaF;
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::setOffsets(LBMReal xoff, LBMReal yoff, LBMReal zoff)
+{
+ this->xoff = xoff;
+ this->yoff = yoff;
+ this->zoff = zoff;
+ this->xoff_sq = xoff * xoff;
+ this->yoff_sq = yoff * yoff;
+ this->zoff_sq = zoff * zoff;
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::interpolateCoarseToFine(D3Q27ICell& icellC, D3Q27ICell& icellF, LBMReal xoff, LBMReal yoff, LBMReal zoff)
+{
+ setOffsets(xoff, yoff, zoff);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, -1, -1, -1);
+ calcInterpolatedNode(icellF.BSW, /*omegaF,*/ -0.25, -0.25, -0.25, calcPressBSW(), -1, -1, -1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, 0.25, -0.25, 1, 1, -1);
+ calcInterpolatedNode(icellF.BNE, /*omegaF,*/ 0.25, 0.25, -0.25, calcPressBNE(), 1, 1, -1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, 0.25, 0.25, -1, 1, 1);
+ calcInterpolatedNode(icellF.TNW, /*omegaF,*/ -0.25, 0.25, 0.25, calcPressTNW(), -1, 1, 1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, 0.25, 1, -1, 1);
+ calcInterpolatedNode(icellF.TSE, /*omegaF,*/ 0.25, -0.25, 0.25, calcPressTSE(), 1, -1, 1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, 0.25, -0.25, -1, 1, -1);
+ calcInterpolatedNode(icellF.BNW, /*omegaF,*/ -0.25, 0.25, -0.25, calcPressBNW(), -1, 1, -1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, 1, -1, -1);
+ calcInterpolatedNode(icellF.BSE, /*omegaF,*/ 0.25, -0.25, -0.25, calcPressBSE(), 1, -1, -1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, -0.25, 0.25, -1, -1, 1);
+ calcInterpolatedNode(icellF.TSW, /*omegaF,*/ -0.25, -0.25, 0.25, calcPressTSW(), -1, -1, 1);
+ calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, 0.25, 0.25, 1, 1, 1);
+ calcInterpolatedNode(icellF.TNE, /*omegaF,*/ 0.25, 0.25, 0.25, calcPressTNE(), 1, 1, 1);
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::interpolateFineToCoarse(D3Q27ICell& icellF, LBMReal* icellC, LBMReal xoff, LBMReal yoff, LBMReal zoff)
+{
+ setOffsets(xoff, yoff, zoff);
+ calcInterpolatedCoefficiets(icellF, omegaF, 2.0, 0, 0, 0, 0, 0, 0);
+ calcInterpolatedNodeFC(icellC, omegaC);
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::calcMoments(const LBMReal* const f, LBMReal omegaInf, LBMReal& press, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3, LBMReal& kxy, LBMReal& kyz, LBMReal& kxz, LBMReal& kxxMyy, LBMReal& kxxMzz)
+{
+ using namespace D3Q27System;
+
+ rho = 0.0;
+ D3Q27System::calcIncompMacroscopicValues(f,rho,vx1,vx2,vx3);
+
+ shearRate = D3Q27System::getShearRate(f, omegaInf);
+
+ LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(omegaInf, shearRate, rho);
+
+ press = rho; //interpolate rho!
+
+ kxy = -3.*omega*((((f[TSW]+f[BNE])-(f[TNW]+f[BSE]))+((f[BSW]+f[TNE])-(f[BNW]+f[TSE])))+((f[SW]+f[NE])-(f[NW]+f[SE]))-(vx1*vx2));// might not be optimal MG 25.2.13
+ kyz = -3.*omega*((((f[BSW]+f[TNE])-(f[TSE]+f[BNW]))+((f[BSE]+f[TNW])-(f[TSW]+f[BNE])))+((f[BS]+f[TN])-(f[TS]+f[BN]))-(vx2*vx3));
+ kxz = -3.*omega*((((f[BNW]+f[TSE])-(f[TSW]+f[BNE]))+((f[BSW]+f[TNE])-(f[BSE]+f[TNW])))+((f[BW]+f[TE])-(f[TW]+f[BE]))-(vx1*vx3));
+ kxxMyy = -3./2.*omega*((((f[D3Q27System::BW]+f[TE])-(f[BS]+f[TN]))+((f[TW]+f[BE])-(f[TS]+f[BN])))+((f[W]+f[E])-(f[S]+f[N]))-(vx1*vx1-vx2*vx2));
+ kxxMzz = -3./2.*omega*((((f[NW]+f[SE])-(f[BS]+f[TN]))+((f[SW]+f[NE])-(f[TS]+f[BN])))+((f[W]+f[E])-(f[B]+f[T]))-(vx1*vx1-vx3*vx3));
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::calcInterpolatedCoefficiets(const D3Q27ICell& icell, LBMReal omega, LBMReal eps_new, LBMReal x, LBMReal y, LBMReal z, LBMReal xs, LBMReal ys, LBMReal zs)
+{
+ LBMReal vx1_SWT,vx2_SWT,vx3_SWT;
+ LBMReal vx1_NWT,vx2_NWT,vx3_NWT;
+ LBMReal vx1_NET,vx2_NET,vx3_NET;
+ LBMReal vx1_SET,vx2_SET,vx3_SET;
+ LBMReal vx1_SWB,vx2_SWB,vx3_SWB;
+ LBMReal vx1_NWB,vx2_NWB,vx3_NWB;
+ LBMReal vx1_NEB,vx2_NEB,vx3_NEB;
+ LBMReal vx1_SEB,vx2_SEB,vx3_SEB;
+
+ LBMReal kxyFromfcNEQ_SWT, kyzFromfcNEQ_SWT, kxzFromfcNEQ_SWT, kxxMyyFromfcNEQ_SWT, kxxMzzFromfcNEQ_SWT;
+ LBMReal kxyFromfcNEQ_NWT, kyzFromfcNEQ_NWT, kxzFromfcNEQ_NWT, kxxMyyFromfcNEQ_NWT, kxxMzzFromfcNEQ_NWT;
+ LBMReal kxyFromfcNEQ_NET, kyzFromfcNEQ_NET, kxzFromfcNEQ_NET, kxxMyyFromfcNEQ_NET, kxxMzzFromfcNEQ_NET;
+ LBMReal kxyFromfcNEQ_SET, kyzFromfcNEQ_SET, kxzFromfcNEQ_SET, kxxMyyFromfcNEQ_SET, kxxMzzFromfcNEQ_SET;
+ LBMReal kxyFromfcNEQ_SWB, kyzFromfcNEQ_SWB, kxzFromfcNEQ_SWB, kxxMyyFromfcNEQ_SWB, kxxMzzFromfcNEQ_SWB;
+ LBMReal kxyFromfcNEQ_NWB, kyzFromfcNEQ_NWB, kxzFromfcNEQ_NWB, kxxMyyFromfcNEQ_NWB, kxxMzzFromfcNEQ_NWB;
+ LBMReal kxyFromfcNEQ_NEB, kyzFromfcNEQ_NEB, kxzFromfcNEQ_NEB, kxxMyyFromfcNEQ_NEB, kxxMzzFromfcNEQ_NEB;
+ LBMReal 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);
+
+ 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;
+ kyzAverage =0;
+ kxzAverage =0;
+ kxxMyyAverage =0;
+ kxxMzzAverage =0;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //
+ // 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;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ LBMReal dxux = ax + 0.5*axx*xs+ 0.25*(axy*ys+axz*zs)+0.0625*axyz*ys*zs;
+ LBMReal dyuy = by + 0.5 * byy * ys + 0.25 * (bxy * xs + byz * zs) + 0.0625 * bxyz * xs * zs;
+ LBMReal dzuz = cz + 0.5 * czz * zs + 0.25 * (cxz * xs + cyz * ys) + 0.0625 * cxyz * xs * ys;
+
+ LBMReal Dxy = bx + 0.5 * bxx * xs + 0.25 * (bxy * ys + bxz * zs) + 0.0625 * bxyz * ys * zs + ay + 0.5 * ayy * ys + 0.25 * (axy * xs + ayz * zs) + 0.0625 * axyz * xs * zs;
+ LBMReal Dxz = cx + 0.5 * cxx * xs + 0.25 * (cxy * ys + cxz * zs) + 0.0625 * cxyz * ys * zs + az + 0.5 * azz * zs + 0.25 * (axz * xs + ayz * ys) + 0.0625 * axyz * xs * ys;
+ LBMReal Dyz = cy + 0.5 * cyy * ys + 0.25 * (cxy * xs + cyz * zs) + 0.0625 * cxyz * xs * zs + bz + 0.5 * bzz * zs + 0.25 * (bxz * xs + byz * ys) + 0.0625 * bxyz * xs * ys;
+
+ shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz);
+
+
+ LBMReal o = Thixotropy::getHerschelBulkleyCollFactorBackward(shearRate, rho); //omega;
+
+ if (o < 0.5)
+ o = 0.5;
+
+ 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_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_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_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 ThixotropyInterpolationProcessor::calcInterpolatedNode(LBMReal* f, /*LBMReal omega,*/ LBMReal x, LBMReal y, LBMReal z, LBMReal press, LBMReal xs, LBMReal ys, LBMReal zs)
+{
+ using namespace D3Q27System;
+
+ LBMReal rho = press ;
+ LBMReal vx1 = a0 + 0.25*( xs*ax + ys*ay + zs*az) + 0.0625*(axx + xs*ys*axy + xs*zs*axz + ayy + ys*zs*ayz + azz) + 0.015625*(xs*ys*zs*axyz);
+ LBMReal vx2 = b0 + 0.25*( xs*bx + ys*by + zs*bz) + 0.0625*(bxx + xs*ys*bxy + xs*zs*bxz + byy + ys*zs*byz + bzz) + 0.015625*(xs*ys*zs*bxyz);
+ LBMReal vx3 = c0 + 0.25*( xs*cx + ys*cy + zs*cz) + 0.0625*(cxx + xs*ys*cxy + xs*zs*cxz + cyy + ys*zs*cyz + czz) + 0.015625*(xs*ys*zs*cxyz);
+
+ LBMReal feq[ENDF+1];
+ D3Q27System::calcIncompFeq(feq,rho,vx1,vx2,vx3);
+
+ f[E] = f_E + xs*x_E + ys*y_E + zs*z_E + xs*ys*xy_E + xs*zs*xz_E + ys*zs*yz_E + feq[E];
+ f[W] = f_E + xs*x_E + ys*y_E + zs*z_E + xs*ys*xy_E + xs*zs*xz_E + ys*zs*yz_E + feq[W];
+ f[N] = f_N + xs*x_N + ys*y_N + zs*z_N + xs*ys*xy_N + xs*zs*xz_N + ys*zs*yz_N + feq[N];
+ f[S] = f_N + xs*x_N + ys*y_N + zs*z_N + xs*ys*xy_N + xs*zs*xz_N + ys*zs*yz_N + feq[S];
+ f[T] = f_T + xs*x_T + ys*y_T + zs*z_T + xs*ys*xy_T + xs*zs*xz_T + ys*zs*yz_T + feq[T];
+ f[B] = f_T + xs*x_T + ys*y_T + zs*z_T + xs*ys*xy_T + xs*zs*xz_T + ys*zs*yz_T + feq[B];
+ f[NE] = f_NE + xs*x_NE + ys*y_NE + zs*z_NE + xs*ys*xy_NE + xs*zs*xz_NE + ys*zs*yz_NE + feq[NE];
+ f[SW] = f_NE + xs*x_NE + ys*y_NE + zs*z_NE + xs*ys*xy_NE + xs*zs*xz_NE + ys*zs*yz_NE + feq[SW];
+ f[SE] = f_SE + xs*x_SE + ys*y_SE + zs*z_SE + xs*ys*xy_SE + xs*zs*xz_SE + ys*zs*yz_SE + feq[SE];
+ f[NW] = f_SE + xs*x_SE + ys*y_SE + zs*z_SE + xs*ys*xy_SE + xs*zs*xz_SE + ys*zs*yz_SE + feq[NW];
+ f[TE] = f_TE + xs*x_TE + ys*y_TE + zs*z_TE + xs*ys*xy_TE + xs*zs*xz_TE + ys*zs*yz_TE + feq[TE];
+ f[BW] = f_TE + xs*x_TE + ys*y_TE + zs*z_TE + xs*ys*xy_TE + xs*zs*xz_TE + ys*zs*yz_TE + feq[BW];
+ f[BE] = f_BE + xs*x_BE + ys*y_BE + zs*z_BE + xs*ys*xy_BE + xs*zs*xz_BE + ys*zs*yz_BE + feq[BE];
+ f[TW] = f_BE + xs*x_BE + ys*y_BE + zs*z_BE + xs*ys*xy_BE + xs*zs*xz_BE + ys*zs*yz_BE + feq[TW];
+ f[TN] = f_TN + xs*x_TN + ys*y_TN + zs*z_TN + xs*ys*xy_TN + xs*zs*xz_TN + ys*zs*yz_TN + feq[TN];
+ f[BS] = f_TN + xs*x_TN + ys*y_TN + zs*z_TN + xs*ys*xy_TN + xs*zs*xz_TN + ys*zs*yz_TN + feq[BS];
+ f[BN] = f_BN + xs*x_BN + ys*y_BN + zs*z_BN + xs*ys*xy_BN + xs*zs*xz_BN + ys*zs*yz_BN + feq[BN];
+ f[TS] = f_BN + xs*x_BN + ys*y_BN + zs*z_BN + xs*ys*xy_BN + xs*zs*xz_BN + ys*zs*yz_BN + feq[TS];
+ f[TNE] = f_TNE + xs*x_TNE + ys*y_TNE + zs*z_TNE + xs*ys*xy_TNE + xs*zs*xz_TNE + ys*zs*yz_TNE + feq[TNE];
+ f[TSW] = f_TSW + xs*x_TSW + ys*y_TSW + zs*z_TSW + xs*ys*xy_TSW + xs*zs*xz_TSW + ys*zs*yz_TSW + feq[TSW];
+ f[TSE] = f_TSE + xs*x_TSE + ys*y_TSE + zs*z_TSE + xs*ys*xy_TSE + xs*zs*xz_TSE + ys*zs*yz_TSE + feq[TSE];
+ f[TNW] = f_TNW + xs*x_TNW + ys*y_TNW + zs*z_TNW + xs*ys*xy_TNW + xs*zs*xz_TNW + ys*zs*yz_TNW + feq[TNW];
+ f[BNE] = f_TSW + xs*x_TSW + ys*y_TSW + zs*z_TSW + xs*ys*xy_TSW + xs*zs*xz_TSW + ys*zs*yz_TSW + feq[BNE];
+ f[BSW] = f_TNE + xs*x_TNE + ys*y_TNE + zs*z_TNE + xs*ys*xy_TNE + xs*zs*xz_TNE + ys*zs*yz_TNE + feq[BSW];
+ f[BSE] = f_TNW + xs*x_TNW + ys*y_TNW + zs*z_TNW + xs*ys*xy_TNW + xs*zs*xz_TNW + ys*zs*yz_TNW + feq[BSE];
+ f[BNW] = f_TSE + xs*x_TSE + ys*y_TSE + zs*z_TSE + xs*ys*xy_TSE + xs*zs*xz_TSE + ys*zs*yz_TSE + feq[BNW];
+ f[ZERO] = f_ZERO + xs*x_ZERO + ys*y_ZERO + zs*z_ZERO + feq[ZERO];
+}
+//////////////////////////////////////////////////////////////////////////
+//Position SWB -0.25, -0.25, -0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressBSW()
+{
+ return press_SWT * (0.140625 + 0.1875 * xoff + 0.1875 * yoff - 0.5625 * zoff) +
+ press_NWT * (0.046875 + 0.0625 * xoff - 0.1875 * yoff - 0.1875 * zoff) +
+ press_SET * (0.046875 - 0.1875 * xoff + 0.0625 * yoff - 0.1875 * zoff) +
+ press_NET * (0.015625 - 0.0625 * xoff - 0.0625 * yoff - 0.0625 * zoff) +
+ press_NEB * (0.046875 - 0.1875 * xoff - 0.1875 * yoff + 0.0625 * zoff) +
+ press_NWB * (0.140625 + 0.1875 * xoff - 0.5625 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.140625 - 0.5625 * xoff + 0.1875 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.421875 + 0.5625 * xoff + 0.5625 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position SWT -0.25, -0.25, 0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressTSW()
+{
+ return press_SWT * (0.421875 + 0.5625 * xoff + 0.5625 * yoff - 0.5625 * zoff) +
+ press_NWT * (0.140625 + 0.1875 * xoff - 0.5625 * yoff - 0.1875 * zoff) +
+ press_SET * (0.140625 - 0.5625 * xoff + 0.1875 * yoff - 0.1875 * zoff) +
+ press_NET * (0.046875 - 0.1875 * xoff - 0.1875 * yoff - 0.0625 * zoff) +
+ press_NEB * (0.015625 - 0.0625 * xoff - 0.0625 * yoff + 0.0625 * zoff) +
+ press_NWB * (0.046875 + 0.0625 * xoff - 0.1875 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.046875 - 0.1875 * xoff + 0.0625 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.140625 + 0.1875 * xoff + 0.1875 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position SET 0.25, -0.25, 0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressTSE()
+{
+ return press_SET * (0.421875 - 0.5625 * xoff + 0.5625 * yoff - 0.5625 * zoff) +
+ press_NET * (0.140625 - 0.1875 * xoff - 0.5625 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.140625 + 0.5625 * xoff + 0.1875 * yoff - 0.1875 * zoff) +
+ press_NWT * (0.046875 + 0.1875 * xoff - 0.1875 * yoff - 0.0625 * zoff) +
+ press_NWB * (0.015625 + 0.0625 * xoff - 0.0625 * yoff + 0.0625 * zoff) +
+ press_NEB * (0.046875 - 0.0625 * xoff - 0.1875 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.046875 + 0.1875 * xoff + 0.0625 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.140625 - 0.1875 * xoff + 0.1875 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position SEB 0.25, -0.25, -0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressBSE()
+{
+ return press_SET * (0.140625 - 0.1875 * xoff + 0.1875 * yoff - 0.5625 * zoff) +
+ press_NET * (0.046875 - 0.0625 * xoff - 0.1875 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.046875 + 0.1875 * xoff + 0.0625 * yoff - 0.1875 * zoff) +
+ press_NWT * (0.015625 + 0.0625 * xoff - 0.0625 * yoff - 0.0625 * zoff) +
+ press_NWB * (0.046875 + 0.1875 * xoff - 0.1875 * yoff + 0.0625 * zoff) +
+ press_NEB * (0.140625 - 0.1875 * xoff - 0.5625 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.140625 + 0.5625 * xoff + 0.1875 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.421875 - 0.5625 * xoff + 0.5625 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position NWB -0.25, 0.25, -0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressBNW()
+{
+ return press_NWT * (0.140625 + 0.1875 * xoff - 0.1875 * yoff - 0.5625 * zoff) +
+ press_NET * (0.046875 - 0.1875 * xoff - 0.0625 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.046875 + 0.0625 * xoff + 0.1875 * yoff - 0.1875 * zoff) +
+ press_SET * (0.015625 - 0.0625 * xoff + 0.0625 * yoff - 0.0625 * zoff) +
+ press_SEB * (0.046875 - 0.1875 * xoff + 0.1875 * yoff + 0.0625 * zoff) +
+ press_NEB * (0.140625 - 0.5625 * xoff - 0.1875 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.140625 + 0.1875 * xoff + 0.5625 * yoff + 0.1875 * zoff) +
+ press_NWB * (0.421875 + 0.5625 * xoff - 0.5625 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position NWT -0.25, 0.25, 0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressTNW()
+{
+ return press_NWT * (0.421875 + 0.5625 * xoff - 0.5625 * yoff - 0.5625 * zoff) +
+ press_NET * (0.140625 - 0.5625 * xoff - 0.1875 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.140625 + 0.1875 * xoff + 0.5625 * yoff - 0.1875 * zoff) +
+ press_SET * (0.046875 - 0.1875 * xoff + 0.1875 * yoff - 0.0625 * zoff) +
+ press_SEB * (0.015625 - 0.0625 * xoff + 0.0625 * yoff + 0.0625 * zoff) +
+ press_NEB * (0.046875 - 0.1875 * xoff - 0.0625 * yoff + 0.1875 * zoff) +
+ press_SWB * (0.046875 + 0.0625 * xoff + 0.1875 * yoff + 0.1875 * zoff) +
+ press_NWB * (0.140625 + 0.1875 * xoff - 0.1875 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position NET 0.25, 0.25, 0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressTNE()
+{
+ return press_NET * (0.421875 - 0.5625 * xoff - 0.5625 * yoff - 0.5625 * zoff) +
+ press_NWT * (0.140625 + 0.5625 * xoff - 0.1875 * yoff - 0.1875 * zoff) +
+ press_SET * (0.140625 - 0.1875 * xoff + 0.5625 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.046875 + 0.1875 * xoff + 0.1875 * yoff - 0.0625 * zoff) +
+ press_SWB * (0.015625 + 0.0625 * xoff + 0.0625 * yoff + 0.0625 * zoff) +
+ press_NWB * (0.046875 + 0.1875 * xoff - 0.0625 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.046875 - 0.0625 * xoff + 0.1875 * yoff + 0.1875 * zoff) +
+ press_NEB * (0.140625 - 0.1875 * xoff - 0.1875 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position NEB 0.25, 0.25, -0.25
+LBMReal ThixotropyInterpolationProcessor::calcPressBNE()
+{
+ return press_NET * (0.140625 - 0.1875 * xoff - 0.1875 * yoff - 0.5625 * zoff) +
+ press_NWT * (0.046875 + 0.1875 * xoff - 0.0625 * yoff - 0.1875 * zoff) +
+ press_SET * (0.046875 - 0.0625 * xoff + 0.1875 * yoff - 0.1875 * zoff) +
+ press_SWT * (0.015625 + 0.0625 * xoff + 0.0625 * yoff - 0.0625 * zoff) +
+ press_SWB * (0.046875 + 0.1875 * xoff + 0.1875 * yoff + 0.0625 * zoff) +
+ press_NWB * (0.140625 + 0.5625 * xoff - 0.1875 * yoff + 0.1875 * zoff) +
+ press_SEB * (0.140625 - 0.1875 * xoff + 0.5625 * yoff + 0.1875 * zoff) +
+ press_NEB * (0.421875 - 0.5625 * xoff - 0.5625 * yoff + 0.5625 * zoff);
+}
+//////////////////////////////////////////////////////////////////////////
+//Position C 0.0, 0.0, 0.0
+void ThixotropyInterpolationProcessor::calcInterpolatedNodeFC(LBMReal* f, LBMReal omega)
+{
+ using namespace D3Q27System;
+
+ LBMReal press = press_NET * (0.125 - 0.25 * xoff - 0.25 * yoff - 0.25 * zoff) +
+ press_NWT * (0.125 + 0.25 * xoff - 0.25 * yoff - 0.25 * zoff) +
+ press_SET * (0.125 - 0.25 * xoff + 0.25 * yoff - 0.25 * zoff) +
+ press_SWT * (0.125 + 0.25 * xoff + 0.25 * yoff - 0.25 * zoff) +
+ press_NEB * (0.125 - 0.25 * xoff - 0.25 * yoff + 0.25 * zoff) +
+ press_NWB * (0.125 + 0.25 * xoff - 0.25 * yoff + 0.25 * zoff) +
+ press_SEB * (0.125 - 0.25 * xoff + 0.25 * yoff + 0.25 * zoff) +
+ press_SWB * (0.125 + 0.25 * xoff + 0.25 * yoff + 0.25 * zoff);
+ LBMReal vx1 = a0;
+ LBMReal vx2 = b0;
+ LBMReal vx3 = c0;
+
+ LBMReal rho = press ;
+
+ LBMReal feq[ENDF+1];
+ D3Q27System::calcIncompFeq(feq,rho,vx1,vx2,vx3);
+
+ LBMReal eps_new = 2.;
+
+
+ LBMReal dxux = ax;
+ LBMReal dyuy = by;
+ LBMReal dzuz = cz;
+
+ LBMReal Dxy = bx + ay;
+ LBMReal Dxz = cx + az;
+ LBMReal Dyz = cy + bz;
+
+ shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz);
+
+
+ LBMReal o = Thixotropy::getHerschelBulkleyCollFactorBackward(shearRate, rho); //omega;
+
+ if (o < 0.5)
+ o = 0.5;
+
+ 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));
+
+ f[E] = f_E + feq[E];
+ f[W] = f_E + feq[W];
+ f[N] = f_N + feq[N];
+ f[S] = f_N + feq[S];
+ f[T] = f_T + feq[T];
+ f[B] = f_T + feq[B];
+ f[NE] = f_NE + feq[NE];
+ f[SW] = f_NE + feq[SW];
+ f[SE] = f_SE + feq[SE];
+ f[NW] = f_SE + feq[NW];
+ f[TE] = f_TE + feq[TE];
+ f[BW] = f_TE + feq[BW];
+ f[BE] = f_BE + feq[BE];
+ f[TW] = f_BE + feq[TW];
+ f[TN] = f_TN + feq[TN];
+ f[BS] = f_TN + feq[BS];
+ f[BN] = f_BN + feq[BN];
+ f[TS] = f_BN + feq[TS];
+ f[TNE] = f_TNE + feq[TNE];
+ f[TNW] = f_TNW + feq[TNW];
+ f[TSE] = f_TSE + feq[TSE];
+ f[TSW] = f_TSW + feq[TSW];
+ f[BNE] = f_TSW + feq[BNE];
+ f[BNW] = f_TSE + feq[BNW];
+ f[BSE] = f_TNW + feq[BSE];
+ f[BSW] = f_TNE + feq[BSW];
+ f[ZERO] = f_ZERO + feq[ZERO];
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::calcInterpolatedVelocity(LBMReal x, LBMReal y, LBMReal z, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3)
+{
+ vx1 = a0 + ax*x + ay*y + az*z + axx*x*x + ayy*y*y + azz*z*z + axy*x*y + axz*x*z + ayz*y*z+axyz*x*y*z;
+ vx2 = b0 + bx*x + by*y + bz*z + bxx*x*x + byy*y*y + bzz*z*z + bxy*x*y + bxz*x*z + byz*y*z+bxyz*x*y*z;
+ vx3 = c0 + cx*x + cy*y + cz*z + cxx*x*x + cyy*y*y + czz*z*z + cxy*x*y + cxz*x*z + cyz*y*z+cxyz*x*y*z;
+}
+//////////////////////////////////////////////////////////////////////////
+void ThixotropyInterpolationProcessor::calcInterpolatedShearStress(LBMReal x, LBMReal y, LBMReal z,LBMReal& tauxx, LBMReal& tauyy, LBMReal& tauzz,LBMReal& tauxy, LBMReal& tauxz, LBMReal& tauyz)
+{
+ tauxx=ax+2*axx*x+axy*y+axz*z+axyz*y*z;
+ tauyy=by+2*byy*y+bxy*x+byz*z+bxyz*x*z;
+ tauzz=cz+2*czz*z+cxz*x+cyz*y+cxyz*x*y;
+ tauxy=0.5*((ay+2.0*ayy*y+axy*x+ayz*z+axyz*x*z)+(bx+2.0*bxx*x+bxy*y+bxz*z+bxyz*y*z));
+ tauxz=0.5*((az+2.0*azz*z+axz*x+ayz*y+axyz*x*y)+(cx+2.0*cxx*x+cxy*y+cxz*z+cxyz*y*z));
+ tauyz=0.5*((bz+2.0*bzz*z+bxz*x+byz*y+bxyz*x*y)+(cy+2.0*cyy*y+cxy*x+cyz*z+cxyz*x*z));
+}
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0ea158b905edbaa3d216da748f384df4de05c3c
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
@@ -0,0 +1,106 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ThixotropyInterpolationProcessor.h
+//! \ingroup BoundarConditions
+//! \author Konstantin Kutscher
+//=======================================================================================
+#ifndef ThixotropyInterpolationProcessor_H_
+#define ThixotropyInterpolationProcessor_H_
+
+#include "InterpolationProcessor.h"
+#include "D3Q27System.h"
+
+//! \brief A class implements an interpolation function of grid refinement for thixotropic fluid.
+
+class ThixotropyInterpolationProcessor : public InterpolationProcessor
+{
+public:
+ ThixotropyInterpolationProcessor();
+ ThixotropyInterpolationProcessor(LBMReal omegaC, LBMReal omegaF);
+ virtual ~ThixotropyInterpolationProcessor();
+ InterpolationProcessorPtr clone();
+ void setOmegas(LBMReal omegaC, LBMReal omegaF);
+ void interpolateCoarseToFine(D3Q27ICell& icellC, D3Q27ICell& icellF);
+ void interpolateCoarseToFine(D3Q27ICell& icellC, D3Q27ICell& icellF, LBMReal xoff, LBMReal yoff, LBMReal zoff);
+ void interpolateFineToCoarse(D3Q27ICell& icellF, LBMReal* icellC);
+ void interpolateFineToCoarse(D3Q27ICell& icellF, LBMReal* icellC, LBMReal xoff, LBMReal yoff, LBMReal zoff);
+ //LBMReal forcingC, forcingF;
+protected:
+private:
+ LBMReal omegaC, omegaF;
+ LBMReal a0, ax, ay, az, axx, ayy, azz, axy, axz, ayz, b0, bx, by, bz, bxx, byy, bzz, bxy, bxz, byz, c0, cx, cy, cz, cxx, cyy, czz, cxy, cxz, cyz, axyz, bxyz, cxyz;
+ LBMReal xoff, yoff, zoff;
+ LBMReal xoff_sq, yoff_sq, zoff_sq;
+ LBMReal press_SWT, press_NWT, press_NET, press_SET, press_SWB, press_NWB, press_NEB, press_SEB;
+
+ LBMReal f_E, f_N, f_T, f_NE, f_SE, f_BE, f_TE, f_TN, f_BN, f_TNE, f_TNW, f_TSE, f_TSW, f_ZERO;
+ LBMReal x_E, x_N, x_T, x_NE, x_SE, x_BE, x_TE, x_TN, x_BN, x_TNE, x_TNW, x_TSE, x_TSW, x_ZERO;
+ LBMReal y_E, y_N, y_T, y_NE, y_SE, y_BE, y_TE, y_TN, y_BN, y_TNE, y_TNW, y_TSE, y_TSW, y_ZERO;
+ LBMReal z_E, z_N, z_T, z_NE, z_SE, z_BE, z_TE, z_TN, z_BN, z_TNE, z_TNW, z_TSE, z_TSW, z_ZERO;
+ LBMReal xy_E, xy_N, xy_T, xy_NE, xy_SE, xy_BE, xy_TE, xy_TN, xy_BN, xy_TNE, xy_TNW, xy_TSE, xy_TSW/*, xy_ZERO*/;
+ LBMReal xz_E, xz_N, xz_T, xz_NE, xz_SE, xz_BE, xz_TE, xz_TN, xz_BN, xz_TNE, xz_TNW, xz_TSE, xz_TSW/*, xz_ZERO*/;
+ LBMReal yz_E, yz_N, yz_T, yz_NE, yz_SE, yz_BE, yz_TE, yz_TN, yz_BN, yz_TNE, yz_TNW, yz_TSE, yz_TSW/*, yz_ZERO*/;
+
+ LBMReal kxyAverage, kyzAverage, kxzAverage, kxxMyyAverage, kxxMzzAverage;
+
+ LBMReal a,b,c;
+
+ LBMReal rho;
+ LBMReal shearRate;
+
+ void setOffsets(LBMReal xoff, LBMReal yoff, LBMReal zoff);
+ void calcMoments(const LBMReal* const f, LBMReal omegaInf, LBMReal& rho, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3,
+ LBMReal& kxy, LBMReal& kyz, LBMReal& kxz, LBMReal& kxxMyy, LBMReal& kxxMzz);
+ void calcInterpolatedCoefficiets(const D3Q27ICell& icell, LBMReal omega, LBMReal eps_new, LBMReal x, LBMReal y, LBMReal z, LBMReal xs, LBMReal ys, LBMReal zs);
+ void calcInterpolatedNode(LBMReal* f, /*LBMReal omega,*/ LBMReal x, LBMReal y, LBMReal z, LBMReal press, LBMReal xs, LBMReal ys, LBMReal zs);
+ LBMReal calcPressBSW();
+ LBMReal calcPressTSW();
+ LBMReal calcPressTSE();
+ LBMReal calcPressBSE();
+ LBMReal calcPressBNW();
+ LBMReal calcPressTNW();
+ LBMReal calcPressTNE();
+ LBMReal calcPressBNE();
+ void calcInterpolatedNodeFC(LBMReal* f, LBMReal omega);
+ void calcInterpolatedVelocity(LBMReal x, LBMReal y, LBMReal z,LBMReal& vx1, LBMReal& vx2, LBMReal& vx3);
+ void calcInterpolatedShearStress(LBMReal x, LBMReal y, LBMReal z,LBMReal& tauxx, LBMReal& tauyy, LBMReal& tauzz,LBMReal& tauxy, LBMReal& tauxz, LBMReal& tauyz);
+};
+
+//////////////////////////////////////////////////////////////////////////
+inline void ThixotropyInterpolationProcessor::interpolateCoarseToFine(D3Q27ICell& icellC, D3Q27ICell& icellF)
+{
+ this->interpolateCoarseToFine(icellC, icellF, 0.0, 0.0, 0.0);
+}
+//////////////////////////////////////////////////////////////////////////
+inline void ThixotropyInterpolationProcessor::interpolateFineToCoarse(D3Q27ICell& icellF, LBMReal* icellC)
+{
+ this->interpolateFineToCoarse(icellF, icellC, 0.0, 0.0, 0.0);
+}
+
+#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.h
index d3210a08467339b937f5e81cfdcbb20b4c0c6c73..b0f492a170a886d23ab1f56d897ac9adf698d8d9 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.h
@@ -10,7 +10,7 @@
class ThixotropyModelLBMKernel;
-//! \brief abstract class for Template Method use Cumulant LBM.
+//! \brief Base class for model of thixotropy based on K16. Use Template Method design pattern for Implementation of different models.
//! \author K. Kutscher, M. Geier
class ThixotropyModelLBMKernel : public LBMKernel
{