From 5865f93f6f735ea2a1cc5869abf30429422f54bc Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Fri, 28 Aug 2020 17:15:00 +0200
Subject: [PATCH] add grid refinement for thixotropic fluid
---
apps/cpu/Applications.cmake | 2 +-
apps/cpu/HerschelBulkleySphere/CMakeLists.txt | 17 -
apps/cpu/HerschelBulkleySphere/hbsphere.cfg | 15 +-
apps/cpu/HerschelBulkleySphere/hbsphere.cpp | 27 +-
apps/cpu/sphere/CMakeLists.txt | 19 +-
apps/cpu/sphere/config.txt | 8 +-
apps/cpu/sphere/sphere.cpp | 99 +--
src/cpu/VirtualFluids.h | 1 +
.../WriteBoundaryConditionsCoProcessor.cpp | 48 +-
.../WriteMacroscopicQuantitiesCoProcessor.cpp | 92 +--
src/cpu/VirtualFluidsCore/LBM/Thixotropy.h | 45 +-
.../LBM/ThixotropyInterpolationProcessor.cpp | 666 ++++++++++++++++++
.../LBM/ThixotropyInterpolationProcessor.h | 106 +++
.../LBM/ThixotropyModelLBMKernel.h | 2 +-
14 files changed, 929 insertions(+), 218 deletions(-)
create mode 100644 src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
create mode 100644 src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 0b92f834e..1e444b79b 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 162bf5604..1f0641fa5 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 1f636ec21..ebcc5d1f8 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 deb44fc21..b437a2736 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 77d7e0a41..022df23ea 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 cbd03fd4b..aedb74bf0 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 872e7b81f..c33a3d696 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 31671c3ed..31937b7ef 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 bae45335e..b7b6246d9 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 1d71e4bbf..e5976c682 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 bb38862b0..235de154c 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 000000000..6a09a1f64
--- /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 000000000..c0ea158b9
--- /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 d3210a084..b0f492a17 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
{
--
GitLab