From 2dda2729f3df7512584d7855471aead6020c3877 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Wed, 3 Jan 2018 15:43:31 +0100
Subject: [PATCH] SpongeLayerBlockVisitor expended only for X1 direction

---
 .../Applications/DLR-F16-Solid/f16-solid.cfg  | 10 ++---
 source/Applications/DLR-F16-Solid/f16.cpp     | 44 ++++++++++++-------
 .../Visitors/SpongeLayerBlockVisitor.cpp      | 20 ++++++++-
 3 files changed, 51 insertions(+), 23 deletions(-)

diff --git a/source/Applications/DLR-F16-Solid/f16-solid.cfg b/source/Applications/DLR-F16-Solid/f16-solid.cfg
index 6db9ce199..f159e1056 100644
--- a/source/Applications/DLR-F16-Solid/f16-solid.cfg
+++ b/source/Applications/DLR-F16-Solid/f16-solid.cfg
@@ -1,4 +1,4 @@
-pathOut = d:/temp/DLR-F16-Solid2
+pathOut = d:/temp/DLR-F16-Solid3
 pathGeo = d:/Projects/SFB880/DLR-F16/A1_Forschungsdaten_Profilgeometrie_STL_CATIA_Rossian
 #fngFileWhole = f16-ascii.stl
 fngFileWhole = grundgeometrie-direkter-export.stl
@@ -48,16 +48,16 @@ refineDistance = 0.6
 
 writeBlocks = true
 
-newStart = flase
-restartStep = 1000
+newStart = true
+restartStep = 19000
 
 cpStep = 1000
-cpStart = 2000
+cpStart = 21000
 
 outTimeStep = 100
 outTimeStart = 100
 
-endTime = 2000
+endTime = 40000
 
 #Cp
 pcpStart = 1000000
diff --git a/source/Applications/DLR-F16-Solid/f16.cpp b/source/Applications/DLR-F16-Solid/f16.cpp
index 7d632b34a..69e1eed8b 100644
--- a/source/Applications/DLR-F16-Solid/f16.cpp
+++ b/source/Applications/DLR-F16-Solid/f16.cpp
@@ -140,6 +140,11 @@ void run(string configname)
       BCAdapterPtr velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
       velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
 
+      fct.SetExpr("U");
+      fct.DefineConst("U", 0.01);
+      BCAdapterPtr velBCAdapterOut(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
+      velBCAdapterOut->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithm()));
+
       BCAdapterPtr outflowBCAdapter(new DensityBCAdapter(rhoLB));
       outflowBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NonReflectingOutflowBCAlgorithm()));
 
@@ -208,15 +213,15 @@ void run(string configname)
 
         
 
-         GbObject3DPtr fngMeshWhole(new GbCylinder3D(15.0, 0.0, 0.0, 15.0, 100.0, 0.0, 25.0));
-         GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut + "/geo/fngMeshWholeCylinder", WbWriterVtkXmlBinary::getInstance());
+         //GbObject3DPtr fngMeshWhole(new GbCylinder3D(15.0, 0.0, 0.0, 15.0, 100.0, 0.0, 25.0));
+         //GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut + "/geo/fngMeshWholeCylinder", WbWriterVtkXmlBinary::getInstance());
 
-         //GbTriFaceMesh3DPtr fngMeshWhole;
-         //if (myid==0) UBLOG(logINFO, "Read fngFileWhole:start");
-         //fngMeshWhole = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(pathGeo+"/"+fngFileWhole, "fngMeshWhole", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
-         //if (myid==0) UBLOG(logINFO, "Read fngFileWhole:end");
-         //fngMeshWhole->rotate(0.0, 0.5, 0.0);
-         //if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut+"/geo/fngMeshWhole2", WbWriterVtkXmlBinary::getInstance());
+         GbTriFaceMesh3DPtr fngMeshWhole;
+         if (myid==0) UBLOG(logINFO, "Read fngFileWhole:start");
+         fngMeshWhole = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(pathGeo+"/"+fngFileWhole, "fngMeshWhole", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
+         if (myid==0) UBLOG(logINFO, "Read fngFileWhole:end");
+         fngMeshWhole->rotate(0.0, 0.5, 0.0);
+         if (myid==0) GbSystem3D::writeGeoObject(fngMeshWhole.get(), pathOut+"/geo/fngMeshWhole2", WbWriterVtkXmlBinary::getInstance());
 
          ////////////////////////////////////////////////////////////////////////////
          //// Zackenband
@@ -260,8 +265,8 @@ void run(string configname)
          }
          //////////////////////////////////////////////////////////////////////////
          Interactor3DPtr fngIntrWhole;
-         fngIntrWhole = D3Q27InteractorPtr(new D3Q27Interactor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
-         //fngIntrWhole = D3Q27TriFaceMeshInteractorPtr(new D3Q27TriFaceMeshInteractor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
+         //fngIntrWhole = D3Q27InteractorPtr(new D3Q27Interactor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
+         fngIntrWhole = D3Q27TriFaceMeshInteractorPtr(new D3Q27TriFaceMeshInteractor(fngMeshWhole, grid, noSlipBCAdapter, Interactor3D::SOLID));//, Interactor3D::POINTS));
 
          //D3Q27TriFaceMeshInteractorPtr triBand1Interactor(new D3Q27TriFaceMeshInteractor(meshBand1, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
          //D3Q27TriFaceMeshInteractorPtr triBand2Interactor(new D3Q27TriFaceMeshInteractor(meshBand2, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
@@ -550,16 +555,17 @@ void run(string configname)
 
          ////sponge layer
          ////////////////////////////////////////////////////////////////////////////
+
+         GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-5.0*blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+         if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance());
+         SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, 1.0);
+         grid->accept(slVisitorX1max);
+
          //GbCuboid3DPtr spongeLayerX1min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_minX1+75, g_maxX2+blockLength, g_maxX3+blockLength));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1min.get(), pathOut+"/geo/spongeLayerX1min", WbWriterVtkXmlASCII::getInstance());
          //SpongeLayerBlockVisitor slVisitorX1min(spongeLayerX1min);
          //grid->accept(slVisitorX1min);
 
-         //GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-75, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
-         //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance());
-         //SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max);
-         //grid->accept(slVisitorX1max);
-
          //GbCuboid3DPtr spongeLayerX3min(new GbCuboid3D(g_minX1-blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_minX3+75));
          //if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX3min.get(), pathOut+"/geo/spongeLayerX3min", WbWriterVtkXmlASCII::getInstance());
          //SpongeLayerBlockVisitor slVisitorX3min(spongeLayerX3min);
@@ -587,6 +593,10 @@ void run(string configname)
          PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
          grid->accept(pqPartVisitor);
 
+         GbCuboid3DPtr spongeLayerX1max(new GbCuboid3D(g_maxX1-5.0*blockLength, g_minX2-blockLength, g_minX3-blockLength, g_maxX1+blockLength, g_maxX2+blockLength, g_maxX3+blockLength));
+         if (myid==0) GbSystem3D::writeGeoObject(spongeLayerX1max.get(), pathOut+"/geo/spongeLayerX1max", WbWriterVtkXmlASCII::getInstance());
+         SpongeLayerBlockVisitor slVisitorX1max(spongeLayerX1max, 1.0);
+         grid->accept(slVisitorX1max);
       }
 
       UbSchedulerPtr nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2]));
@@ -613,9 +623,9 @@ void run(string configname)
       //TimeseriesCoProcessor tsp1(grid, stepMV, mic1, pathOut+"/mic/mic1", comm);
 
       const std::shared_ptr<ConcreteCalculatorFactory> calculatorFactory = std::make_shared<ConcreteCalculatorFactory>(stepSch);
-      CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::MPI));
+      //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::MPI));
       //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, stepSch));
-      //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, tavSch, CalculationManager::PrePostBc));
+      CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, calculatorFactory, CalculatorType::PREPOSTBC));
 
       if (myid==0) UBLOG(logINFO, "Simulation-start");
       calculation->calculate();
diff --git a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
index 1dc76eda5..21510e9e7 100644
--- a/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
+++ b/source/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
@@ -40,7 +40,25 @@ void SpongeLayerBlockVisitor::visit(Grid3DPtr grid, Block3DPtr block)
       if (boundingBox->isCellInsideGbObject3D(minX1, minX2, minX3, maxX1, maxX2, maxX3))
       {
          ILBMKernelPtr kernel = block->getKernel();
-         kernel->setCollisionFactor(collFactor);
+         double oldCollFactor = kernel->getCollisionFactor();
+
+         int ibX1 = block->getX1();
+
+         UbTupleInt3 ixMin = grid->getBlockIndexes(boundingBox->getX1Minimum(),boundingBox->getX2Minimum(),boundingBox->getX3Minimum());
+         UbTupleInt3 ixMax = grid->getBlockIndexes(boundingBox->getX1Maximum(),boundingBox->getX2Maximum(),boundingBox->getX3Maximum());
+
+         int ibMax = val<1>(ixMax)-val<1>(ixMin)+1;
+         double index = ibX1-val<1>(ixMin)+1;
+         
+         //double a = oldCollFactor-1.0 / (double)ibMax;
+         //double b = val<1>(ixMin)-ibX1+1;
+         //double a = oldCollFactor-1.0 / (double)ibMax;
+         //double a = oldCollFactor-1.0 / (double)ibMax;
+
+
+         double newCollFactor = oldCollFactor - (oldCollFactor-1.0)/ibMax*index;
+
+         kernel->setCollisionFactor(newCollFactor);
       }
    }
 }
-- 
GitLab