diff --git a/CMake/cmake_config_files/BOMBADIL.config.cmake b/CMake/cmake_config_files/BOMBADIL.config.cmake
index 9c4bd4ecffab1e63161343ecc493eb9d9bc951a4..0534bf1c87cb54bf6d8839b881fe6b772e2b7365 100644
--- a/CMake/cmake_config_files/BOMBADIL.config.cmake
+++ b/CMake/cmake_config_files/BOMBADIL.config.cmake
@@ -48,25 +48,6 @@ set(LIGGGHTS_RELEASE_LIBRARY "d:/Tools/LIGGGHTS/build/Release/liggghts.lib")
# SET(METIS_RELEASE_LIBRARY "/mnt/d/Tools/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a")
#ENDIF()
-#################################################################################
-# PE
-#################################################################################
-IF(${USE_DEM_COUPLING})
- SET(PE_BINARY_DIR "d:/Tools/waLBerla/walberlaGit/build" CACHE PATH "pe binary dir")
- SET(PE_ROOT "d:/Tools/waLBerla/walberlaGit" CACHE PATH "pe root")
-
- SET(PE_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/pe/Debug/pe.lib)
- SET(PE_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/pe/Release/pe.lib)
- SET(BLOCKFOREST_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/blockforest/Debug/blockforest.lib)
- SET(BLOCKFOREST_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/blockforest/Release/blockforest.lib)
- SET(DOMAIN_DECOMPOSITION_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/domain_decomposition/Debug/domain_decomposition.lib)
- SET(DOMAIN_DECOMPOSITION_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/domain_decomposition/Release/domain_decomposition.lib)
- SET(GEOMETRY_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/geometry/Debug/geometry.lib)
- SET(GEOMETRY_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/geometry/Release/geometry.lib)
- SET(CORE_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/core/Debug/core.lib)
- SET(CORE_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/core/Release/core.lib)
-
- ENDIF()
##################################################################################
# FETOL
diff --git a/CMake/cmake_config_files/PHOENIX.config.cmake b/CMake/cmake_config_files/PHOENIX.config.cmake
index d31d8684a53a769e48408ad5febe7d2c6b22c623..2112bd6aa50e9335bc6b23bda0f0e9fda3ef7533 100644
--- a/CMake/cmake_config_files/PHOENIX.config.cmake
+++ b/CMake/cmake_config_files/PHOENIX.config.cmake
@@ -4,24 +4,6 @@
# OS: CentOS 7.3
#################################################################################
-#################################################################################
-# PE (legacy)
-#################################################################################
-IF(${USE_DEM_COUPLING})
- SET(PE_BINARY_DIR "/home/irmb/walberla-git/build" CACHE PATH "pe binary dir")
- SET(PE_ROOT "/home/irmb/walberla-git" CACHE PATH "pe root")
-
- SET(PE_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/pe/libpe.a)
- SET(PE_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/pe/libpe.a)
- SET(BLOCKFOREST_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/blockforest/libblockforest.a)
- SET(BLOCKFOREST_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/blockforest/libblockforest.a)
- SET(DOMAIN_DECOMPOSITION_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/domain_decomposition/libdomain_decomposition.a)
- SET(DOMAIN_DECOMPOSITION_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/domain_decomposition/libdomain_decomposition.a)
- SET(GEOMETRY_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/geometry/libgeometry.a)
- SET(GEOMETRY_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/geometry/libgeometry.a)
- SET(CORE_DEBUG_LIBRARY ${PE_BINARY_DIR}/src/core/libcore.a)
- SET(CORE_RELEASE_LIBRARY ${PE_BINARY_DIR}/src/core/libcore.a)
-ENDIF()
## nvidia
set(CMAKE_CUDA_ARCHITECTURES 60) # NVIDIA Tesla P100
diff --git a/apps/cpu/Thermoplast/CMakeLists.txt b/apps/cpu/Thermoplast/CMakeLists.txt
deleted file mode 100644
index 5624b03136a7c901d1fe69fa224464104d3a08a6..0000000000000000000000000000000000000000
--- a/apps/cpu/Thermoplast/CMakeLists.txt
+++ /dev/null
@@ -1,29 +0,0 @@
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
-
-########################################################
-## C++ PROJECT ###
-########################################################
-PROJECT(thermoplast)
-IF(${USE_DEM_COUPLING})
- INCLUDE(${APPS_ROOT}/IncludsList.cmake)
- INCLUDE(${SOURCE_ROOT}/DemCoupling/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 ${CAB_ADDITIONAL_LINK_LIBRARIES} VirtualFluids)
-
- #message("CAB_ADDITIONAL_LINK_LIBRARIES: " ${CAB_ADDITIONAL_LINK_LIBRARIES})
-
- #################################################################
- ### CREATE PROJECT ###
- #################################################################
- CREATE_CAB_PROJECT(thermoplast BINARY)
-ENDIF()
\ No newline at end of file
diff --git a/apps/cpu/Thermoplast/config.txt b/apps/cpu/Thermoplast/config.txt
deleted file mode 100644
index 18e7e7ea896acdc20527de1d6af9154fbac16e6b..0000000000000000000000000000000000000000
--- a/apps/cpu/Thermoplast/config.txt
+++ /dev/null
@@ -1,57 +0,0 @@
-#simulation parameters
-
-#x1min x2min x3min x1max x2max x3max
-#boundingBox = 0 0 0 300 1520 2320
-
-#boundingBox = 60 1370 130 190 1530 320 #test bb
-
-boundingBox = 60 20 130 190 170 320 #test bb
-
-#boundingBox = 60 0 10 190 1530 750 #test bb 2
-
-#boundingBox = 60 0 10 190 1530 2320 #production bb
-
-blocknx = 10 10 10
-#blocknx = 300 420 320
-availMem = 25e9
-#uLB = 0.1
-uLB = 0.03
-Re = 300
-
-#PE parameters
-#test pe offset
-peMinOffset = 46 2 2
-peMaxOffset = -8 -25 -2
-
-#production pe offset
-#peMinOffset = 46 18 14
-#peMaxOffset = -8 -25 -23
-
-sphereTime = 10
-
-#geometry files
-pathGeo = d:/Projects/ThermoPlast/SimPerfMS
-michel = /Werkzeug_Michel_MS.stl
-plexiglas = /plexiglas.stl
-
-#obstacle
-obstacle = true
-obstacleGeo1 = /QuaderMS.stl #/DreieckMS_2.stl # DreieckSchoen.iges.stl #/QuaderMS.stl
-obstacleGeo2 = /KugelMS_2.stl
-obstacleGeo3 = /DreieckMS.stl
-
-pathOut = g:/temp/thermoplastObst
-
-logToFile = false
-
-#restart
-restart = false
-restartStep = 1000
-
-#timing
-nupsTime = 100 100 1000000
-cpStart = 1000
-cpStep = 1000
-outTime = 1000
-endTime = 100000
-
diff --git a/apps/cpu/Thermoplast/thermoplast.cpp b/apps/cpu/Thermoplast/thermoplast.cpp
deleted file mode 100644
index b004543254e4788ffa6a7314536b1811c217700b..0000000000000000000000000000000000000000
--- a/apps/cpu/Thermoplast/thermoplast.cpp
+++ /dev/null
@@ -1,763 +0,0 @@
-#include <iostream>
-#include <string>
-
-#include "PointerDefinitions.h"
-
-#include <iostream>
-#include <string>
-#include <memory>
-#include <array>
-
-#include "VirtualFluids.h"
-#include <muParser.h>
-#include "ForceCalculator.h"
-
-
-#include <MovableObjectInteractor.h>
-#include <DemCoProcessor.h>
-#include <PePartitioningGridVisitor.h>
-
-#include <PePhysicsEngineMaterialAdapter.h>
-#include <PePhysicsEngineGeometryAdapter.h>
-#include <PePhysicsEngineSolverAdapter.h>
-#include "PeLoadBalancerAdapter.h"
-
-#include <VelocityBcReconstructor.h>
-#include <EquilibriumReconstructor.h>
-#include <ExtrapolationReconstructor.h>
-
-#include <DummyPhysicsEngineSolverAdapter.h>
-#include <DummyPhysicsEngineMaterialAdapter.h>
-#include <DummyPhysicsEngineGeometryAdapter.h>
-#include <WriteDemObjectsCoProcessor.h>
-#include <WritePeBlocksCoProcessor.h>
-
-#include "CreateDemObjectsCoProcessor.h"
-#include "RestartDemObjectsCoProcessor.h"
-
-using namespace std;
-
-//simulation bounding box
-double g_minX1 = 0;
-double g_minX2 = 0;
-double g_minX3 = 0;
-
-double g_maxX1 = 0;
-double g_maxX2 = 0;
-double g_maxX3 = 0;
-
-vector<double> peMinOffset;
-vector<double> peMaxOffset;
-
-string pathOut;// = "d:/temp/thermoplastCluster";
-string pathGeo;// = "d:/Projects/ThermoPlast/Geometrie";
-
-void addNozzle(SPtr<Grid3D> grid, SPtr<vf::mpi::Communicator> comm, SPtr<BCAdapter> noSlipBCAdapter/*, InteractorsHelper& intHelper*/)
-{
- int myid = comm->getProcessID();
- if (myid==0) UBLOG(logINFO, "Add nozzles:start");
-
- SPtr<UbScheduler> sch(new UbScheduler(1));
- WriteGbObjectsCoProcessor gbObjectsCoProcessor(grid, sch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
-
- std::vector< SPtr<Interactor3D> > interactors;
-
- for (int i = 0; i <= 55; i++)
- {
- SPtr<GbTriFaceMesh3D> bbGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/n_bb_new/bb_new"+UbSystem::toString(i)+".stl", "bb", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<Interactor3D> bbInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(bbGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
- //GbSystem3D::writeGeoObject(bbGeo.get(), pathOut+"/ns/bbGeo"+UbSystem::toString(i), WbWriterVtkXmlBinary::getInstance());
- //intHelper.addInteractor(bbInt);
- if (myid==0) gbObjectsCoProcessor.addGbObject(bbGeo);
- interactors.push_back(bbInt);
- }
-
- for (int i = 0; i <= 334; i++)
- {
- SPtr<GbTriFaceMesh3D> bbGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/n_bb/bb"+UbSystem::toString(i)+".stl", "bb", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<Interactor3D> bbInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(bbGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
- //GbSystem3D::writeGeoObject(bbGeo.get(), pathOut+"/ns/bbGeo"+UbSystem::toString(i), WbWriterVtkXmlBinary::getInstance());
- //intHelper.addInteractor(bbInt);
- if (myid==0) gbObjectsCoProcessor.addGbObject(bbGeo);
- interactors.push_back(bbInt);
- }
-
- for (int i = 0; i <= 51; i++)
- {
- SPtr<GbTriFaceMesh3D> bsGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/n_bs/bs"+UbSystem::toString(i)+".stl", "bs", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<Interactor3D> bsInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(bsGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
- //intHelper.addInteractor(bsInt);
- if (myid==0) gbObjectsCoProcessor.addGbObject(bsGeo);
- interactors.push_back(bsInt);
- }
-
- std::array<int, 6> n ={ 0,1,3,4,6,7 };
-
- for (int i = 0; i < n.size(); i++)
- {
- SPtr<GbTriFaceMesh3D> biGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/n_bi/bi"+UbSystem::toString(n[i])+".stl", "bi", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- SPtr<Interactor3D> biInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(biGeo, grid, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES));
- //intHelper.addInteractor(biInt);
- if (myid==0) gbObjectsCoProcessor.addGbObject(biGeo);
- interactors.push_back(biInt);
- }
-
- if (myid==0) gbObjectsCoProcessor.process(0);
-
-
- for (SPtr<Interactor3D> interactor : interactors)
- {
- std::vector< std::shared_ptr<Block3D> > blockVector;
- UbTupleInt3 blockNX=grid->getBlockNX();
- SPtr<GbObject3D> geoObject(interactor->getGbObject3D());
- double ext = 0.0;
- std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
- grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
- for (std::shared_ptr<Block3D> block : blockVector)
- {
- if (block->getKernel())
- {
- interactor->setBCBlock(block);
- }
- }
- interactor->initInteractor();
- }
-
- if (myid==0) UBLOG(logINFO, "Add nozzles:end");
-}
-
-std::shared_ptr<DemCoProcessor> makePeCoProcessor(SPtr<Grid3D> grid, SPtr<vf::mpi::Communicator> comm, const SPtr<UbScheduler> peScheduler, const std::shared_ptr<LBMUnitConverter> lbmUnitConverter, int maxpeIterations)
-{
- double peRelaxtion = 0.7;
- //int maxpeIterations = 10000;
- //Beschleunigung g
- double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
- //Vector3D globalLinearAcc(0.0, -g, 0.0);
- //Vector3D globalLinearAcc(0.0, 0.0, -g);
- Vector3D globalLinearAcc(0.0, 0.0, 0.0);
-
- std::shared_ptr<PePhysicsEngineMaterialAdapter> planeMaterial = std::make_shared<PePhysicsEngineMaterialAdapter>("granular", 1.0, 0, 0.1 / 2, 0.1 / 2, 0.5, 1, 1, 0, 0);
-
- const int gridNX1 = val<1>(grid->getBlockNX()) * grid->getNX1();
- const int gridNX2 = val<2>(grid->getBlockNX()) * grid->getNX2();
- const int gridNX3 = val<3>(grid->getBlockNX()) * grid->getNX3();
-
- //UbTupleInt3 simulationDomain(gridNx, gridNy, gridNz);
- //std::array<double, 6> simulationDomain = {1, 1, 1, 30, 30, 30};
- std::array<double, 6> simulationDomain ={ g_minX1, g_minX2, g_minX3, g_minX1+gridNX1, g_minX2+gridNX2, g_minX3+gridNX3 };
- UbTupleInt3 numberOfBlocks(grid->getNX1(), grid->getNX2(), grid->getNX3());
- //UbTupleInt3 numberOfBlocks((simulationDomain[3]-simulationDomain[0])/val<1>(grid->getBlockNX()), (simulationDomain[4]-simulationDomain[1])/val<2>(grid->getBlockNX()), (simulationDomain[5]-simulationDomain[2])/val<3>(grid->getBlockNX()));
- UbTupleBool3 isPeriodic(grid->isPeriodicX1(), grid->isPeriodicX2(), grid->isPeriodicX3());
- Vector3D minOffset(peMinOffset[0], peMinOffset[1], peMinOffset[2]);
- Vector3D maxOffset(peMaxOffset[0], peMaxOffset[1], peMaxOffset[2]);
-
- SPtr<GbObject3D> boxPE(new GbCuboid3D(simulationDomain[0]+minOffset[0], simulationDomain[1]+minOffset[1], simulationDomain[2]+minOffset[2], simulationDomain[3]+maxOffset[0], simulationDomain[4]+maxOffset[1], simulationDomain[5]+maxOffset[2]));
- GbSystem3D::writeGeoObject(boxPE.get(), pathOut + "/geo/boxPE", WbWriterVtkXmlBinary::getInstance());
-
- std::shared_ptr<PeParameter> peParamter = std::make_shared<PeParameter>(peRelaxtion, maxpeIterations, globalLinearAcc,
- planeMaterial, simulationDomain, numberOfBlocks, isPeriodic, minOffset, maxOffset);
- std::shared_ptr<PeLoadBalancerAdapter> loadBalancer(new PeLoadBalancerAdapter(grid, comm->getNumberOfProcesses(), comm->getProcessID()));
- std::shared_ptr<PhysicsEngineSolverAdapter> peSolver = std::make_shared<PePhysicsEngineSolverAdapter>(peParamter, loadBalancer);
- //create obstacle
- //test
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->createObstacle(Vector3D( 90, 260, 472), Vector3D( 115, 320, 460));
- //production
- //std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->createObstacle(Vector3D( 90, 430, 472), Vector3D( 115, 320, 460));
- //std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->createObstacle(Vector3D( 100, 430, 1840), Vector3D( 130, 320, 470));
- //std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->createObstacle(Vector3D( 100, 821, 1159), Vector3D( 125, 625, 625));
- //walberla::pe::createSphere(*globalBodyStorage, *forest, *storageId, 0, walberla::pe::Vec3( -720, 820, 1150), 900, global, communicating, infiniteMass);
- //walberla::pe::createSphere(*globalBodyStorage, *forest, *storageId, 0, walberla::pe::Vec3( -720, 220, 472), 900, material, global, communicating, infiniteMass);
-
- SPtr<CoProcessor> peblocks(new WritePeBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm, std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(peSolver)->getBlockForest()));
- peblocks->process(0);
- peblocks.reset();
-
- const std::shared_ptr<ForceCalculator> forceCalculator = std::make_shared<ForceCalculator>(comm);
-
- return std::make_shared<DemCoProcessor>(grid, peScheduler, comm, forceCalculator, peSolver);
-}
-
-void createSpheres(double radius, Vector3D origin, int maxX2, int maxX3, double uLB, SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor)
-{
- double d = 2.0*radius;
- double dividerX2 = (double)maxX2/2.0;
- double dividerX3 = (double)maxX3/2.0;
- for (int x3 = 0; x3 < maxX3; x3++)
- for (int x2 = 0; x2 < maxX2; x2++)
- //for (int x1 = 0; x1 < 1; x1++)
- {
- //SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+2.0*d*(double)x1, origin[1]+(double)x2*1.0*d, origin[2]+(double)x3*1.0*d, radius));
- SPtr<GbObject3D> sphere(new GbSphere3D(origin[0]+2.0*d, origin[1]+(double)x2*1.0*d, origin[2]+(double)x3*1.0*d, radius));
- createSphereCoProcessor->addGeoObject(sphere, Vector3D(uLB, -uLB+uLB/dividerX2*(double)x2, -uLB+uLB/dividerX3*(double)x3));
- }
-}
-
-void thermoplast(string configname)
-{
- SPtr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
- int myid = comm->getProcessID();
-
- vf::basics::ConfigurationFile config;
- config.load(configname);
-
- vector<int> blocknx = config.getVector<int>("blocknx");
- vector<double> boundingBox = config.getVector<double>("boundingBox");
-
- int endTime = config.getValue<int>("endTime");
- double outTime = config.getValue<double>("outTime");
- double availMem = config.getValue<double>("availMem");
- double uLB = config.getValue<double>("uLB");
- double Re = config.getValue<double>("Re");
-
- string michel = config.getValue<string>("michel");
- string plexiglas = config.getValue<string>("plexiglas");
- double sphereTime = config.getValue<double>("sphereTime");
-
- double cpStart = config.getValue<double>("cpStart");
- double cpStep = config.getValue<double>("cpStep");
- bool restart = config.getValue<bool>("restart");
- int restartStep = config.getValue<int>("restartStep");
-
- peMinOffset = config.getVector<double>("peMinOffset");
- peMaxOffset = config.getVector<double>("peMaxOffset");
-
- pathOut = config.getValue<string>("pathOut");
- pathGeo = config.getValue<string>("pathGeo");
-
- vector<int> nupsTime = config.getVector<int>("nupsTime");
-
- bool logToFile = config.getValue<bool>("logToFile");
- if (logToFile)
- {
-#if defined(__unix__)
- if (myid==0)
- {
- const char* str = pathOut.c_str();
- mkdir(str, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH);
- }
-#endif
-
- if (myid==0)
- {
- stringstream logFilename;
- logFilename<<pathOut+"/logfile"+UbSystem::getTimeStamp()+".txt";
- UbLog::output_policy::setStream(logFilename.str());
- }
- }
-
- bool obstacle = config.getValue<bool>("obstacle");
- string obstacleGeo1 = config.getValue<string>("obstacleGeo1");
- string obstacleGeo2 = config.getValue<string>("obstacleGeo2");
- string obstacleGeo3 = config.getValue<string>("obstacleGeo3");
-
- if (myid==0) UBLOG(logINFO, "BEGIN LOGGING - " << UbSystem::getTimeStamp());
-
- //parameters
- //string pathOut = "d:/temp/thermoplast3";
- //string pathGeo = "d:/Projects/ThermoPlast/Geometrie";
- int numOfThreads = 1;
- //int blocknx[3] ={ 10,10,10 };
- //double endTime = 1000000;
- //double outTime = 300;
- //double availMem = 8e9;
- double deltax = 1;
- double rhoLB = 0.0;
- //double uLB = 0.1;
- double radiusLB = 7.5;
- double radiusWorld = 1.5e-3;
- //double nuLB = 0.000333333;
- //double Re = (uLB*2.0*radiusLB)/nuLB;
- //double Re = 900;
- double nuLB = (uLB*2.0*radiusLB)/Re;
-
- //geometry definition
-
- //simulation bounding box
- g_minX1 = boundingBox[0];
- g_minX2 = boundingBox[1];
- g_minX3 = boundingBox[2];
-
- g_maxX1 = boundingBox[3];
- g_maxX2 = boundingBox[4];
- g_maxX3 = boundingBox[5];
-
- double blockLength = blocknx[0]*deltax;
-
- //Grid definition
- SPtr<Grid3D> grid(new Grid3D(comm));
- grid->setDeltaX(deltax);
- grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
- grid->setPeriodicX1(false);
- grid->setPeriodicX2(false);
- grid->setPeriodicX3(false);
-
- //boundary conditions definition
- //////////////////////////////////////////////////////////////////////////////
- SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
- //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new ThinWallNoSlipBCAlgorithm()));
-
- mu::Parser fct;
- fct.SetExpr("U");
- fct.DefineConst("U", uLB);
- SPtr<BCAdapter> inflowAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
- //inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
-
- SPtr<BCAdapter> outflowAdapter(new DensityBCAdapter(rhoLB));
- outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new EqDensityBCAlgorithm()));
- //outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
- //outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm()));
-
- //sphere BC
- mu::Parser fct2;
- fct2.SetExpr("U");
- fct2.DefineConst("U", 0.0);
- SPtr<BCAdapter> velocityBcParticleAdapter(new VelocityBCAdapter(true, false, false, fct2, 0, BCFunction::INFCONST));
- velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
-
- //boundary conditions visitor
- SPtr<BoundaryConditionsBlockVisitor> bcVisitor(new BoundaryConditionsBlockVisitor());
- bcVisitor->addBC(noSlipBCAdapter);
- bcVisitor->addBC(inflowAdapter);
- bcVisitor->addBC(outflowAdapter);
- bcVisitor->addBC(velocityBcParticleAdapter);
- //////////////////////////////////////////////////////////////////////////////////
-
- //LBM kernel definition
- SPtr<LBMKernel> kernel;
- kernel = SPtr<LBMKernel>(new IncompressibleCumulantLBMKernel());
- //SPtr<BCProcessor> bcProc(new BCProcessor());
- SPtr<BCProcessor> bcProc(new ThinWallBCProcessor());
- kernel->setBCProcessor(bcProc);
-
- //if (myid==0) UBLOG(logINFO, "Read obstacleGeo1:start");
- //SPtr<GbTriFaceMesh3D> obstacleGeo1geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+obstacleGeo1, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- //if (myid==0) UBLOG(logINFO, "Read obstacleGeo1:end");
- //if (myid==0) GbSystem3D::writeGeoObject(obstacleGeo1geo.get(), pathOut+"/geo/obstacleGeo1", WbWriterVtkXmlBinary::getInstance());
- //g_minX1 = obstacleGeo1geo->getX1Minimum();
- //g_minX2 = obstacleGeo1geo->getX2Minimum();
- //g_minX3 = obstacleGeo1geo->getX3Minimum();
- //g_maxX1 = obstacleGeo1geo->getX1Maximum();
- //g_maxX2 = obstacleGeo1geo->getX2Maximum();
- //g_maxX3 = obstacleGeo1geo->getX3Maximum();
-
-
-
- //blocks generating
- SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
- if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathOut + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
- GenBlocksGridVisitor genBlocks(gridCube);
- grid->accept(genBlocks);
-
-
- //{
- //SPtr<Interactor3D> obstacleGeo1int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(obstacleGeo1geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- //SPtr<Grid3DVisitor> peVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
- //InteractorsHelper intHelper(grid, peVisitor, true);
- //intHelper.addInteractor(obstacleGeo1int);
- //intHelper.selectBlocks();
-
- ////create LBM kernel
- ////SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, 1);
- ////grid->accept(kernelVisitor);
-
- ////SPtr<Interactor3D> obstacleGeo1int = SPtr<D3Q27Interactor>(new D3Q27Interactor(obstacleGeo1geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- ////UBLOG(logINFO, "Obst: start");
- ////std::vector< std::shared_ptr<Block3D> > blockVector;
- ////UbTupleInt3 blockNX=grid->getBlockNX();
- ////SPtr<GbObject3D> geoObject(obstacleGeo1int->getGbObject3D());
- ////double ext = 0.0;
- ////std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
- ////grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
- ////for (std::shared_ptr<Block3D> block : blockVector)
- ////{
- ////if (block->getKernel())
- ////{
- ////obstacleGeo1int->setBCBlock(block);
- ////}
- ////}
- ////UBLOG(logINFO, "Obst: select blocks");
- ////obstacleGeo1int->initInteractor();
- ////UBLOG(logINFO, "Obst: end");
-
- //SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm));
- //ppblocks->process(0);
- //ppblocks.reset();
- //}
-
- //return;
-
-
- /////////////////////////////////////////////////////
- ////PE domain test
- //std::array<double, 6> simulationDomain ={ g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3 };
- //Vector3D minOffset(peMinOffset[0], peMinOffset[1], peMinOffset[2]);
- //Vector3D maxOffset(peMaxOffset[0], peMaxOffset[1], peMaxOffset[2]);
- //SPtr<GbObject3D> boxPE(new GbCuboid3D(simulationDomain[0]+minOffset[0], simulationDomain[1]+minOffset[1], simulationDomain[2]+minOffset[2], simulationDomain[3]+maxOffset[0], simulationDomain[4]+maxOffset[1], simulationDomain[5]+maxOffset[2]));
- //GbSystem3D::writeGeoObject(boxPE.get(), pathOut + "/geo/boxPE", WbWriterVtkXmlBinary::getInstance());
- //return;
- //////////////////////////////////////////////////////
-
-
- if (myid == 0)
- {
- UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "* uLB = " << uLB);
- UBLOG(logINFO, "* rhoLB = " << rhoLB);
- UBLOG(logINFO, "* nuLB = " << nuLB);
- UBLOG(logINFO, "* deltaX = " << deltax);
- UBLOG(logINFO, "* radius = " << radiusLB);
- UBLOG(logINFO, "* Re = " << Re);
- UBLOG(logINFO, "* number of threads = "<<numOfThreads);
- UBLOG(logINFO, "* number of processes = "<<comm->getNumberOfProcesses());
- UBLOG(logINFO, "* path = "<<pathOut);
- UBLOG(logINFO, "Preprocess - start");
- }
-
- //GbCuboid3DPtr geoInjector2(new GbCuboid3D(-12, -5, 1210, 64, 105, 1320));
- //if (myid == 0) GbSystem3D::writeGeoObject(geoInjector2.get(), pathOut + "/geo/geoInjector2", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr geoInjector5(new GbCuboid3D(-12, 1415, 205, 64, 1525, 315));
- //if (myid == 0) GbSystem3D::writeGeoObject(geoInjector5.get(), pathOut + "/geo/geoInjector5", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr geoInjector4(new GbCuboid3D(-12, -5, 205, 64, 105, 315));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInjector4.get(), pathOut + "/geo/geoInjector4", WbWriterVtkXmlASCII::getInstance());
-
- //GbCuboid3DPtr geoInjector7(new GbCuboid3D(28, 705, 542, 103, 815, 652));
- //if (myid == 0) GbSystem3D::writeGeoObject(geoInjector7.get(), pathOut + "/geo/geoInjector7", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr testWallGeo(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(testWallGeo.get(), pathOut + "/geo/testWallGeo", WbWriterVtkXmlASCII::getInstance());
-
- if (!restart)
- {
- //box
- SPtr<GbObject3D> box(new GbCuboid3D(g_minX1-blockLength, g_minX2, g_minX3, g_maxX1+blockLength, g_maxX2, g_maxX3));
- GbSystem3D::writeGeoObject(box.get(), pathOut + "/geo/box", WbWriterVtkXmlBinary::getInstance());
-
- //michel
- if (myid==0) UBLOG(logINFO, "Read michelGeo:start");
- SPtr<GbTriFaceMesh3D> michelGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+michel, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read michelGeo:end");
- if (myid==0) GbSystem3D::writeGeoObject(michelGeo.get(), pathOut+"/geo/michelGeo", WbWriterVtkXmlBinary::getInstance());
-
- //plexiglas
- if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:start");
- SPtr<GbTriFaceMesh3D> plexiglasGeo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+plexiglas, "plexiglasGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read plexiglasGeo:end");
- if (myid==0) GbSystem3D::writeGeoObject(plexiglasGeo.get(), pathOut+"/geo/plexiglasGeo", WbWriterVtkXmlBinary::getInstance());
-
- //inflow
- GbCuboid3DPtr geoOutflowMichel(new GbCuboid3D(g_minX1-blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflowMichel.get(), pathOut + "/geo/geoOutflowMichel", WbWriterVtkXmlASCII::getInstance());
-
- //outflow
- GbCuboid3DPtr geoOutflowPlexiglas(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflowPlexiglas.get(), pathOut + "/geo/geoOutflowPlexiglas", WbWriterVtkXmlASCII::getInstance());
-
- //set boundary conditions for blocks and create process decomposition for MPI
- SPtr<D3Q27Interactor> boxInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
-
- //inflow
- //SPtr<D3Q27Interactor> inflowInjector2Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector2, grid, inflowAdapter, Interactor3D::SOLID));
- //SPtr<D3Q27Interactor> inflowInjector5Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector5, grid, inflowAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> inflowInjector4Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector4, grid, inflowAdapter, Interactor3D::SOLID));
- //SPtr<D3Q27Interactor> inflowInjector7Int = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInjector7, grid, inflowAdapter, Interactor3D::SOLID));
-
- SPtr<D3Q27Interactor> outflowMichelInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflowMichel, grid, outflowAdapter, Interactor3D::SOLID));
-
- //outflow
- SPtr<D3Q27Interactor> outflowPlexiglasInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflowPlexiglas, grid, outflowAdapter, Interactor3D::SOLID));
-
- //michel
- SPtr<Interactor3D> michelInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(michelGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- //plexiglas
- SPtr<Interactor3D> plexiglasInt = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(plexiglasGeo, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- SPtr<D3Q27Interactor> testWallInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(testWallGeo, grid, inflowAdapter, Interactor3D::SOLID));
-
- SPtr<Interactor3D> obstacleGeo1int, obstacleGeo2int, obstacleGeo3int;
- if (obstacle)
- {
- //obstacleGeo1
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo1:start");
- SPtr<GbTriFaceMesh3D> obstacleGeo1geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+obstacleGeo1, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo1:end");
- if (myid==0) GbSystem3D::writeGeoObject(obstacleGeo1geo.get(), pathOut+"/geo/obstacleGeo1", WbWriterVtkXmlBinary::getInstance());
- obstacleGeo1int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(obstacleGeo1geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- //obstacleGeo2
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo2:start");
- SPtr<GbTriFaceMesh3D> obstacleGeo2geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+obstacleGeo2, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo2:end");
- if (myid==0) GbSystem3D::writeGeoObject(obstacleGeo2geo.get(), pathOut+"/geo/obstacleGeo2", WbWriterVtkXmlBinary::getInstance());
- obstacleGeo2int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(obstacleGeo2geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- //obstacleGeo3
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo3:start");
- SPtr<GbTriFaceMesh3D> obstacleGeo3geo = SPtr<GbTriFaceMesh3D>(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+obstacleGeo3, "michelGeo", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- if (myid==0) UBLOG(logINFO, "Read obstacleGeo3:end");
- if (myid==0) GbSystem3D::writeGeoObject(obstacleGeo3geo.get(), pathOut+"/geo/obstacleGeo3", WbWriterVtkXmlBinary::getInstance());
- obstacleGeo3int = SPtr<D3Q27TriFaceMeshInteractor>(new D3Q27TriFaceMeshInteractor(obstacleGeo3geo, grid, noSlipBCAdapter, Interactor3D::SOLID));
- }
-
- //////////////////////////////////////////////////////////////////////////
- //SPtr<Grid3DVisitor> peVisitor(new PePartitioningGridVisitor(comm, demCoProcessor));
- SPtr<Grid3DVisitor> peVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
- InteractorsHelper intHelper(grid, peVisitor, true);
-
- //intHelper.addInteractor(obstacleGeo1int);
-
- intHelper.addInteractor(boxInt);
- intHelper.addInteractor(michelInt);
- intHelper.addInteractor(plexiglasInt);
- //intHelper.addInteractor(inflowInjector2Int);
- //intHelper.addInteractor(inflowInjector5Int);
- intHelper.addInteractor(inflowInjector4Int);
- //intHelper.addInteractor(inflowInjector7Int);
- intHelper.addInteractor(outflowPlexiglasInt);
- intHelper.addInteractor(outflowMichelInt);
- intHelper.addInteractor(obstacleGeo1int);
- intHelper.addInteractor(obstacleGeo2int);
- intHelper.addInteractor(obstacleGeo3int);
- //intHelper.addInteractor(testWallInt);
- intHelper.selectBlocks();
-
- //write data for visualization of block grid
- SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathOut, WbWriterVtkXmlBinary::getInstance(), comm));
- ppblocks->process(0);
- ppblocks.reset();
-
- unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
- int ghostLayer = 3;
- unsigned long long numberOfNodesPerBlock = (unsigned long long)(blocknx[0])* (unsigned long long)(blocknx[1])* (unsigned long long)(blocknx[2]);
- unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
- unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
- double needMemAll = double(numberOfNodesPerBlockWithGhostLayer*(27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
- double needMem = needMemAll / double(comm->getNumberOfProcesses());
-
- if (myid == 0)
- {
- UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
- UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
- int minInitLevel = grid->getCoarsestInitializedLevel();
- int maxInitLevel = grid->getFinestInitializedLevel();
- for (int level = minInitLevel; level <= maxInitLevel; level++)
- {
- int nobl = grid->getNumberOfBlocks(level);
- UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
- UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl*numberOfNodesPerBlock);
- }
- UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
- UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
- UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
- }
-
- //create LBM kernel
- SetKernelBlockVisitor kernelVisitor(kernel, nuLB, availMem, needMem);
- grid->accept(kernelVisitor);
-
- addNozzle(grid, comm, noSlipBCAdapter/*,intHelper*/);
-
- intHelper.setBC();
-
-
- ////////////////////////////////////////////////////////////////////////////////////////////////////
- //{
- ////UBLOG(logINFO, "Obst: start, rank="<<myid);
- //std::vector< std::shared_ptr<Block3D> > blockVector;
- //UbTupleInt3 blockNX=grid->getBlockNX();
- //SPtr<GbObject3D> geoObject(obstacleGeo3int->getGbObject3D());
- //double ext = 0.0;
- //std::array<double, 6> AABB ={ geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum() };
- //grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext, AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext, AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
- //for (std::shared_ptr<Block3D> block : blockVector)
- //{
- //if (block->getKernel())
- //{
- //obstacleGeo3int->setBCBlock(block);
- //}
- //}
- //UBLOG(logINFO, "Obst: select blocks, number of blocks="<<blockVector.size()<<", rank="<<myid);
- //obstacleGeo3int->initInteractor();
- //UBLOG(logINFO, "Obst: end, rank="<<myid);
- //}
- //////////////////////////////////////////////////////////////////////////////////////////////////////
- //initialization of distributions
- InitDistributionsBlockVisitor initVisitor;
- //initVisitor.setVx1(uLB);
- grid->accept(initVisitor);
-
- //write data for visualization of boundary conditions
- {
- //SPtr<UbScheduler> geoSch(new UbScheduler(1));
- //WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), comm);
- //ppgeo.process(0);
-
- //WriteMacroscopicQuantitiesCoProcessor ppInit(grid, geoSch, pathOut, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
- //ppInit.process(0);
- }
-
- if (myid == 0) UBLOG(logINFO, "Preprocess - end");
- }
- //restart
- //UBLOG(logINFO, "restart definition - start, rank="<<myid);
- SPtr<UbScheduler> restartSch(new UbScheduler(cpStep, cpStart));
- //SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(grid, restartSch, pathOut, comm));
- SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, restartSch, pathOut, comm));
- restartCoProcessor->setLBMKernel(kernel);
- restartCoProcessor->setBCProcessor(bcProc);
-
- if (restart)
- {
- //restartStep = restartCoProcessor->readCpTimeStep();
- restartCoProcessor->restart(restartStep);
- }
-
- //PE initialization
- double refLengthLb = radiusLB*2.0;
- double refLengthWorld = radiusWorld*2.0;
- const std::shared_ptr<LBMUnitConverter> lbmUnitConverter = std::make_shared<LBMUnitConverter>(refLengthWorld, LBMUnitConverter::WORLD_MATERIAL::AIR_20C, refLengthLb);
- if (myid == 0) std::cout << lbmUnitConverter->toString() << std::endl;
- double rhoSphere = 915 * lbmUnitConverter->getFactorDensityWToLb(); // kg/m^3
- if (myid == 0) UBLOG(logINFO, "rhoSphere = "<<rhoSphere);
- SPtr<PhysicsEngineMaterialAdapter> sphereMaterial(new PePhysicsEngineMaterialAdapter("Polypropylen", rhoSphere, 0, 0.15, 0.1, 0.45, 0.5, 1, 0, 0));
- const int timestep = 2;
- const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
- int maxpeIterations = 10;//endTime/2;
- SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter, maxpeIterations);
- demCoProcessor->setBlockVisitor(bcVisitor);
-
- ////////////////////////////////////////////////////////////////////////////
- ////generating spheres
- //UBLOG(logINFO, "generating spheres - start, rank="<<myid);
- SPtr<UbScheduler> sphereScheduler(new UbScheduler(sphereTime/*10,10,10*/));
- double toleranz = 0.0;//0.05;
- SPtr<CreateDemObjectsCoProcessor> createSphereCoProcessor(new CreateDemObjectsCoProcessor(grid, sphereScheduler, comm, demCoProcessor, sphereMaterial, toleranz));
- //UBLOG(logINFO, "generating spheres - stop, rank="<<myid);
-
- ////restart
- ////UBLOG(logINFO, "restart definition - start, rank="<<myid);
- //SPtr<UbScheduler> restartSch(new UbScheduler(cpStep, cpStart));
- ////SPtr<MPIIORestartCoProcessor> restartCoProcessor(new MPIIORestartCoProcessor(grid, restartSch, pathOut, comm));
- //SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, restartSch, pathOut, comm));
- //restartCoProcessor->setLBMKernel(kernel);
- //restartCoProcessor->setBCProcessor(bcProc);
- SPtr<RestartDemObjectsCoProcessor> restartDemObjectsCoProcessor(new RestartDemObjectsCoProcessor(grid, restartSch, pathOut, demCoProcessor, createSphereCoProcessor, radiusLB, comm));
- //UBLOG(logINFO, "restart definition - stop, rank="<<myid);
-
- if (restart)
- {
- createSphereCoProcessor->setToleranz(0.05);
- restartDemObjectsCoProcessor->restart(restartStep);
- createSphereCoProcessor->setToleranz(toleranz);
- }
-
- //set connectors
- //UBLOG(logINFO, "set connectors - start, rank="<<myid);
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- grid->accept(setConnsVisitor);
- //UBLOG(logINFO, "set connectors - stop, rank="<<myid);
-
- //BC visitor
- //UBLOG(logINFO, "BC visitor - start, rank="<<myid);
- grid->accept(*bcVisitor.get());
- //UBLOG(logINFO, "BC visitor - stop, rank="<<myid);
-
- //sphere prototypes
- //UBLOG(logINFO, "sphere prototypes - start, rank="<<myid);
- double d = 2.0*radiusLB;
- int maxX2 = 5;
- int maxX3 = 5;
- //Vector3D origin1(g_minX1+peMinOffset[0]-1.5*d, geoInjector5->getX2Minimum()+1.4*d-6.0, geoInjector5->getX3Minimum()+1.5*d);
- //createSpheres(radiusLB, origin1, maxX2, maxX3, uLB, createSphereCoProcessor);
- //Vector3D origin2(g_minX1+peMinOffset[0]-1.5*d, geoInjector2->getX2Minimum()+2.2*d, geoInjector2->getX3Minimum()+1.5*d);
- //createSpheres(radiusLB, origin2, maxX2, maxX3, uLB, createSphereCoProcessor);
-
- Vector3D origin2(g_minX1+peMinOffset[0]-1.5*d, geoInjector4->getX2Minimum()+2.4*d, geoInjector4->getX3Minimum()+1.5*d);
- createSpheres(radiusLB,origin2,maxX2,maxX3,uLB,createSphereCoProcessor);
-
- //maxX2 = 7;
- //maxX3 = 7;
- //Vector3D origin3(g_minX1+peMinOffset[0]-1.5*d, geoInjector7->getX2Minimum()+0.5*d, geoInjector7->getX3Minimum()+0.5*d);
- //createSpheres(radiusLB,origin3,maxX2,maxX3,uLB,createSphereCoProcessor);
-
-
- createSphereCoProcessor->process(0);
-
- //write data for visualization of macroscopic quantities
- SPtr<UbScheduler> visSch(new UbScheduler(outTime));
- SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathOut,
- WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
-
- SPtr<WriteBoundaryConditionsCoProcessor> writeBCCoProcessor(new WriteBoundaryConditionsCoProcessor(grid, visSch, pathOut,
- WbWriterVtkXmlBinary::getInstance(), comm));
-
- SPtr<WriteDemObjectsCoProcessor> writeDemObjectsCoProcessor(new WriteDemObjectsCoProcessor(grid, visSch, pathOut, WbWriterVtkXmlBinary::getInstance(), demCoProcessor, comm));
-
- if (!restart)
- {
- writeMQCoProcessor->process(0);
- writeBCCoProcessor->process(0);
- writeDemObjectsCoProcessor->process(0);
- }
- ////performance control
- SPtr<UbScheduler> nupsSch(new UbScheduler(nupsTime[0], nupsTime[1], nupsTime[2]));
- SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
-
- //start simulation
- //omp_set_num_threads(numOfThreads);
- SPtr<UbScheduler> stepGhostLayer(peScheduler);
- SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
-
- calculator->addCoProcessor(npr);
- calculator->addCoProcessor(createSphereCoProcessor);
- calculator->addCoProcessor(demCoProcessor);
- ////calculator->addCoProcessor(writeBCCoProcessor);
- calculator->addCoProcessor(writeDemObjectsCoProcessor);
- calculator->addCoProcessor(writeMQCoProcessor);
- calculator->addCoProcessor(restartDemObjectsCoProcessor);
- calculator->addCoProcessor(restartCoProcessor);
-
- if (myid == 0) UBLOG(logINFO, "Simulation-start");
- calculator->calculate();
- if (myid == 0) UBLOG(logINFO, "Simulation-end");
- if (myid==0) UBLOG(logINFO, "END LOGGING - " << UbSystem::getTimeStamp());
-}
-
-//////////////////////////////////////////////////////////////////////////
-int main(int argc, char* argv[])
-{
- try
- {
- //Sleep(30000);
- walberla::Environment env(argc, argv);
-
- if (argv!=NULL)
- {
- //if (argv[1]!=NULL)
- //{
- //thermoplast(string("thermoplast.cfg"));
- thermoplast(string("d:/Projects/VirtualFluidsGit/source/Applications/Thermoplast/config.txt"));
- //}
- //else
- //{
- //cout<<"Configuration file must be set!: "<<argv[0]<<" <config file>"<<endl<<std::flush;
- //}
- }
- return 0;
- }
- catch (std::exception& e)
- {
- UBLOG(logERROR, e.what());
- }
- catch (std::string& s)
- {
- UBLOG(logERROR, s);
- }
- catch (...)
- {
- UBLOG(logERROR, "unknown exception");
- }
-}
diff --git a/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.cpp b/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.cpp
deleted file mode 100644
index 6f32a053afb4f6f45ae74e32b7f8665ab4fd58db..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.cpp
+++ /dev/null
@@ -1,121 +0,0 @@
-#include "CreateDemObjectsCoProcessor.h"
-#include <mpi/Communicator.h>
-#include "DemCoProcessor.h"
-#include "EquilibriumReconstructor.h"
-#include "ExtrapolationReconstructor.h"
-#include "GbSphere3D.h"
-#include "Grid3D.h"
-#include "LBMReconstructor.h"
-#include "MovableObjectInteractor.h"
-#include "NoSlipBCAlgorithm.h"
-#include "PePhysicsEngineMaterialAdapter.h"
-#include "PhysicsEngineMaterialAdapter.h"
-#include "SetBcBlocksBlockVisitor.h"
-#include "UbScheduler.h"
-#include "VelocityBCAdapter.h"
-#include "VelocityBCAlgorithm.h"
-#include "VelocityBcReconstructor.h"
-#include "VelocityWithDensityBCAlgorithm.h"
-#include "muParser.h"
-
-CreateDemObjectsCoProcessor::CreateDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
- std::shared_ptr<vf::mpi::Communicator> comm,
- SPtr<DemCoProcessor> demCoProcessor,
- SPtr<PhysicsEngineMaterialAdapter> demObjectMaterial,
- double tolerance)
- : CoProcessor(grid, s), comm(comm), demCoProcessor(demCoProcessor), demObjectMaterial(demObjectMaterial),
- tolerance(tolerance)
-{
- mu::Parser fct;
- fct.SetExpr("U");
- fct.DefineConst("U", 0.0);
- velocityBcParticleAdapter =
- SPtr<BCAdapter>(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- velocityBcParticleAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
-
- // const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
- std::shared_ptr<Reconstructor> equilibriumReconstructor(new EquilibriumReconstructor());
- // const std::shared_ptr<Reconstructor> lbmReconstructor(new LBMReconstructor(false));
- extrapolationReconstructor = SPtr<Reconstructor>(new ExtrapolationReconstructor(equilibriumReconstructor));
- demCounter = 0;
-}
-//////////////////////////////////////////////////////////////////////////
-void CreateDemObjectsCoProcessor::process(double step)
-{
- if (scheduler->isDue(step)) {
- int istep = static_cast<int>(step);
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "CreateDemObjectsCoProcessor::process start step: " << istep);
- timer.resetAndStart();
-#endif
-
- createGeoObjects();
-
-#ifdef TIMING
- // if (comm->isRoot()) UBLOG(logINFO, "createGeoObjects() time = "<<timer.stop()<<" s");
- // if (comm->isRoot()) UBLOG(logINFO, "number of objects = "<<(int)(geoObjectPrototypeVector.size()));
- // if (comm->isRoot()) UBLOG(logINFO, "total number of objects = "<<demCounter);
- if (comm->isRoot())
- UBLOG(logINFO, "CreateDemObjectsCoProcessor::process stop step: " << istep);
-#endif
-
- // demCoProcessor->distributeIDs();
-
- //#ifdef TIMING
- // if (comm->isRoot()) UBLOG(logINFO, "demCoProcessor->distributeIDs() time = "<<timer.stop()<<" s");
- //#endif
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void CreateDemObjectsCoProcessor::addGeoObject(SPtr<GbObject3D> geoObjectPrototype, Vector3D initalVelocity)
-{
- geoObjectPrototypeVector.push_back(geoObjectPrototype);
- this->initalVelocity.push_back(initalVelocity);
-}
-
-void CreateDemObjectsCoProcessor::clearGeoObjects()
-{
- geoObjectPrototypeVector.clear();
- initalVelocity.clear();
-}
-
-void CreateDemObjectsCoProcessor::createGeoObjects()
-{
- int size = (int)(geoObjectPrototypeVector.size());
-
- std::vector<std::shared_ptr<Block3D>> blockVector;
-
- for (int i = 0; i < size; i++) {
- SPtr<GbSphere3D> sphere = std::dynamic_pointer_cast<GbSphere3D>(geoObjectPrototypeVector[i]);
- if (demCoProcessor->isSpheresIntersection(sphere->getX1Centroid(), sphere->getX2Centroid(),
- sphere->getX3Centroid(),
- sphere->getRadius() * 2.0 * (1.0 - tolerance))) {
- continue;
- }
-
- SPtr<GbObject3D> geoObject((GbObject3D *)(geoObjectPrototypeVector[i]->clone()));
- SPtr<MovableObjectInteractor> geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(
- geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
- demCoProcessor->addInteractor(geoObjectInt, demObjectMaterial, initalVelocity[i]);
- demCounter++;
- }
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "createGeoObjects() time = " << timer.stop() << " s");
- if (comm->isRoot())
- UBLOG(logINFO, "number of objects = " << (int)(geoObjectPrototypeVector.size()));
- if (comm->isRoot())
- UBLOG(logINFO, "total number of objects = " << demCounter);
- // if (comm->isRoot()) UBLOG(logINFO, "CreateDemObjectsCoProcessor::process stop step: " << istep);
-#endif
-
- demCoProcessor->distributeIDs();
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "demCoProcessor->distributeIDs() time = " << timer.stop() << " s");
-#endif
-}
diff --git a/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.h b/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.h
deleted file mode 100644
index 7da317e67bd932f7d594c68d63ebc117b50c1e85..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/CreateDemObjectsCoProcessor.h
+++ /dev/null
@@ -1,52 +0,0 @@
-#ifndef CreateSphereCoProcessor_h__
-#define CreateSphereCoProcessor_h__
-
-#include "CoProcessor.h"
-#include "Vector3D.h"
-#include <array>
-#include <vector>
-
-//#define TIMING
-
-#ifdef TIMING
-#include "UbTiming.h"
-#endif
-
-class Grid3D;
-class UbScheduler;
-namespace vf::mpi {class Communicator;}
-class DemCoProcessor;
-class GbObject3D;
-class BCAdapter;
-class Reconstructor;
-class PhysicsEngineMaterialAdapter;
-
-class CreateDemObjectsCoProcessor : public CoProcessor
-{
-public:
- CreateDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, std::shared_ptr<vf::mpi::Communicator> comm,
- SPtr<DemCoProcessor> demCoProcessor,
- SPtr<PhysicsEngineMaterialAdapter> geoObjectMaterial, double tolerance = 0);
- void process(double step) override;
- void addGeoObject(SPtr<GbObject3D> geoObjectPrototype, Vector3D initalVelocity);
- void clearGeoObjects();
- void createGeoObjects();
- double getToleranz() const { return tolerance; }
- void setToleranz(double val) { tolerance = val; }
-
-protected:
-private:
- std::shared_ptr<vf::mpi::Communicator> comm;
- SPtr<DemCoProcessor> demCoProcessor;
- std::vector<SPtr<GbObject3D>> geoObjectPrototypeVector;
- SPtr<PhysicsEngineMaterialAdapter> demObjectMaterial;
- std::vector<Vector3D> initalVelocity;
- SPtr<BCAdapter> velocityBcParticleAdapter;
- SPtr<Reconstructor> extrapolationReconstructor;
- int demCounter;
- double tolerance;
-#ifdef TIMING
- UbTimer timer;
-#endif
-};
-#endif // CreateSphereCoProcessor_h__
diff --git a/src/cpu/DemCoupling/DemCoProcessor.cpp b/src/cpu/DemCoupling/DemCoProcessor.cpp
deleted file mode 100644
index 642a942d7d96b73af898690a5737f53d2d88b1a5..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/DemCoProcessor.cpp
+++ /dev/null
@@ -1,515 +0,0 @@
-#include "DemCoProcessor.h"
-
-#include "BCProcessor.h"
-#include <mpi/Communicator.h>
-#include "DataSet3D.h"
-#include "DistributionArray3D.h"
-#include "ForceCalculator.h"
-#include "GbSphere3D.h"
-#include "Grid3D.h"
-#include "ILBMKernel.h"
-#include "MovableObjectInteractor.h"
-#include "SetBcBlocksBlockVisitor.h"
-#include "UbScheduler.h"
-
-#include "PePhysicsEngineGeometryAdapter.h"
-#include "PePhysicsEngineSolverAdapter.h"
-#include "PhysicsEngineGeometryAdapter.h"
-#include "PhysicsEngineMaterialAdapter.h"
-#include "PhysicsEngineSolverAdapter.h"
-
-#include "BCArray3D.h"
-#include "Block3D.h"
-#include "BoundaryConditions.h"
-#include "BoundaryConditionsBlockVisitor.h"
-#include "MPICommunicator.h"
-
-#include "UbLogger.h"
-
-#include <array>
-#include <functional>
-
-DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, std::shared_ptr<vf::mpi::Communicator> comm,
- std::shared_ptr<ForceCalculator> forceCalculator,
- std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver,
- double intermediatePeSteps)
- : CoProcessor(grid, s), comm(comm), forceCalculator(forceCalculator),
- physicsEngineSolver(std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)),
- intermediateDemSteps(intermediatePeSteps)
-{
-#ifdef TIMING
- timer.resetAndStart();
-#endif
-
- std::shared_ptr<walberla::blockforest::BlockForest> forest =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
- std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
-
- for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
- walberla::pe::Storage *storage = blockIt->getData<walberla::pe::Storage>(*storageId.get());
- walberla::pe::BodyStorage *bodyStorage = &(*storage)[0];
- walberla::pe::BodyStorage *bodyStorageShadowCopies = &(*storage)[1];
-
- bodyStorage->registerAddCallback("DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::addPeGeo), this));
- bodyStorage->registerRemoveCallback("DemCoProcessor",
- std::bind1st(std::mem_fun(&DemCoProcessor::removePeGeo), this));
-
- bodyStorageShadowCopies->registerAddCallback("DemCoProcessor",
- std::bind1st(std::mem_fun(&DemCoProcessor::addPeShadowGeo), this));
- bodyStorageShadowCopies->registerRemoveCallback(
- "DemCoProcessor", std::bind1st(std::mem_fun(&DemCoProcessor::removePeShadowGeo), this));
- }
-}
-
-DemCoProcessor::~DemCoProcessor()
-{
- std::shared_ptr<walberla::blockforest::BlockForest> forest =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
- std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
-
- for (auto ¤tBlock : *forest) {
- walberla::pe::Storage *storage = currentBlock.getData<walberla::pe::Storage>(*storageId.get());
- walberla::pe::BodyStorage &localStorage = (*storage)[0];
- walberla::pe::BodyStorage &shadowStorage = (*storage)[1];
-
- localStorage.clearAddCallbacks();
- localStorage.clearRemoveCallbacks();
-
- shadowStorage.clearAddCallbacks();
- shadowStorage.clearRemoveCallbacks();
- }
-}
-
-void DemCoProcessor::addInteractor(std::shared_ptr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial,
- Vector3D initalVelocity)
-{
- interactors.push_back(interactor);
- const int id = static_cast<int>(interactors.size() - 1);
- interactor->setID(id);
- const auto peGeometryAdapter = this->createPhysicsEngineGeometryAdapter(interactor, physicsEngineMaterial);
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)->isActive()) {
- peGeometryAdapter->setLinearVelolocity(initalVelocity);
- geoIdMap.insert(
- std::make_pair(std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)->getSystemID(),
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(peGeometryAdapter)));
- }
- SetBcBlocksBlockVisitor setBcVisitor(interactor);
- grid->accept(setBcVisitor);
-
- // std::vector< std::shared_ptr<Block3D> > blockVector;
- // UbTupleInt3 blockNX=grid->getBlockNX();
- // SPtr<GbObject3D> geoObject(interactor->getGbObject3D());
- // double ext = 0.0;
- // std::array<double, 6> AABB ={
- // geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum()
- // }; grid->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext,
- // AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext,
- // AABB[5]+(double)val<3>(blockNX)*ext, blockVector); for (std::shared_ptr<Block3D> block : blockVector)
- //{
- // if (block->getKernel())
- // {
- // interactor->setBCBlock(block);
- // //UBLOG(logINFO, "DemCoProcessor::addInteractor() rank = "<<comm->getProcessID());
- // }
- //}
-
- interactor->initInteractor();
-
- physicsEngineGeometrieAdapters.push_back(peGeometryAdapter);
-}
-
-std::shared_ptr<PhysicsEngineGeometryAdapter> DemCoProcessor::createPhysicsEngineGeometryAdapter(
- std::shared_ptr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial) const
-{
- const int id = static_cast<int>(interactors.size() - 1);
- SPtr<GbSphere3D> vfSphere = std::static_pointer_cast<GbSphere3D>(interactor->getGbObject3D());
- const Vector3D position(vfSphere->getX1Centroid(), vfSphere->getX2Centroid(), vfSphere->getX3Centroid());
- auto peGeometryAdapter = this->physicsEngineSolver->createPhysicsEngineGeometryAdapter(
- id, position, vfSphere->getRadius(), physicsEngineMaterial);
- interactor->setPhysicsEngineGeometry(peGeometryAdapter);
- return peGeometryAdapter;
-}
-
-void DemCoProcessor::process(double actualTimeStep)
-{
-#ifdef TIMING
- timer.resetAndStart();
-#endif
-
- this->applyForcesOnGeometries();
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::process start step: " << actualTimeStep);
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::applyForcesOnGeometries() time = " << timer.stop() << " s");
-#endif
-
- if (scheduler->isDue(actualTimeStep)) {
- // UBLOG(logINFO, "DemCoProcessor::update - START - timestep = " << actualTimeStep);
- const double demTimeStepsPerIteration = scheduler->getMinStep();
-
- if (demTimeStepsPerIteration != 1)
- this->scaleForcesAndTorques(1.0 / demTimeStepsPerIteration);
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::scaleForcesAndTorques() time = " << timer.stop() << " s");
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::calculateDemTimeStep():");
-#endif
-
- if (this->intermediateDemSteps == 1)
- this->calculateDemTimeStep(demTimeStepsPerIteration);
-
- //#ifdef TIMING
- // if (comm->isRoot()) UBLOG(logINFO, "DemCoProcessor::calculateDemTimeStep() time = "<<timer.stop()<<"
- // s");
- //#endif
- // if ((int)actualTimeStep % 100 == 0)
- //{
- // if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometries[0])->isActive())
- // {
- // //UBLOG(logINFO, "v: (x,y,z) " << physicsEngineGeometries[0]->getLinearVelocity() << " actualTimeStep
- // = " << UbSystem::toString(actualTimeStep));
- // }
- //}
-
- // during the intermediate time steps of the collision response, the currently acting forces
- // (interaction forces, gravitational force, ...) have to remain constant.
- // Since they are reset after the call to collision response, they have to be stored explicitly before.
- // Then they are set again after each intermediate step.
-
- this->moveVfGeoObjects();
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::moveVfGeoObject() time = " << timer.stop() << " s");
-#endif
-
- grid->accept(*boundaryConditionsBlockVisitor.get());
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "grid->accept(*boundaryConditionsBlockVisitor.get()) time = " << timer.stop() << " s");
-#endif
-
- // UBLOG(logINFO, "DemCoProcessor::update - END - timestep = " << actualTimeStep);
- }
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, "DemCoProcessor::process stop step: " << actualTimeStep);
-#endif
-}
-//////////////////////////////////////////////////////////////////////////
-std::shared_ptr<PhysicsEngineSolverAdapter> DemCoProcessor::getPhysicsEngineSolver() { return physicsEngineSolver; }
-
-void DemCoProcessor::applyForcesOnGeometries()
-{
- for (int i = 0; i < physicsEngineGeometrieAdapters.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- this->setForcesToObject(grid, interactors[i], physicsEngineGeometrieAdapters[i]);
-
- // physicsEngineGeometries[i]->setLinearVelolocity(Vector3D(-0.001, 0.0, 0.0));
- // physicsEngineGeometries[i]->setAngularVelocity(Vector3D(0.01, 0.01, 0.01));
- // UBLOG(logINFO, "v: (x,y,z) " << physicsEngineGeometries[i]->getLinearVelocity());
- }
- }
-}
-
-void DemCoProcessor::setForcesToObject(SPtr<Grid3D> grid, SPtr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry)
-{
- for (BcNodeIndicesMap::value_type t : interactor->getBcNodeIndicesMap()) {
- SPtr<Block3D> block = t.first;
- SPtr<ILBMKernel> kernel = block->getKernel();
- SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
- distributions->swap();
-
- std::set<std::vector<int>> &transNodeIndicesSet = t.second;
- for (std::vector<int> node : transNodeIndicesSet) {
- int x1 = node[0];
- int x2 = node[1];
- int x3 = node[2];
-
- if (kernel->isInsideOfDomain(x1, x2, x3) && bcArray->isFluid(x1, x2, x3)) {
- // TODO: calculate assumed boundary position
-
- const Vector3D worldCoordinates = grid->getNodeCoordinates(block, x1, x2, x3);
- const auto boundaryVelocity = physicsEngineGeometry->getVelocityAtPosition(worldCoordinates);
-
- SPtr<BoundaryConditions> bc = bcArray->getBC(x1, x2, x3);
- const Vector3D force = forceCalculator->getForces(x1, x2, x3, distributions, bc, boundaryVelocity);
- physicsEngineGeometry->addForceAtPosition(force, worldCoordinates);
- }
- }
- distributions->swap();
- }
-}
-
-void DemCoProcessor::scaleForcesAndTorques(double scalingFactor)
-{
- for (int i = 0; i < physicsEngineGeometrieAdapters.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- const Vector3D force = physicsEngineGeometrieAdapters[i]->getForce() * scalingFactor;
- const Vector3D torque = physicsEngineGeometrieAdapters[i]->getTorque() * scalingFactor;
-
- physicsEngineGeometrieAdapters[i]->resetForceAndTorque();
-
- physicsEngineGeometrieAdapters[i]->setForce(force);
- physicsEngineGeometrieAdapters[i]->setTorque(torque);
-
- // UBLOG(logINFO, "F: (x,y,z) " << force);
- // UBLOG(logINFO, "T: (x,y,z) " << torque);
- }
- }
-}
-
-void DemCoProcessor::calculateDemTimeStep(double step)
-{
- physicsEngineSolver->runTimestep(step);
-
-#ifdef TIMING
- if (comm->isRoot())
- UBLOG(logINFO, " physicsEngineSolver->runTimestep() time = " << timer.stop() << " s");
-#endif
-}
-
-void DemCoProcessor::moveVfGeoObjects()
-{
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->getSemiactive()) {
- walberla::pe::RigidBody *peGeoObject = getPeGeoObject(
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->getSystemID());
- if (peGeoObject != nullptr) {
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->setGeometry(peGeoObject);
- interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->setSemiactive(false);
- } else {
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->setInactive();
- }
- } else {
- interactors[i]->moveGbObjectTo(physicsEngineGeometrieAdapters[i]->getPosition());
- }
- }
- }
-}
-
-bool DemCoProcessor::isDemObjectInAABB(std::array<double, 6> AABB)
-{
- bool result = false;
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- SPtr<GbObject3D> geoObject = interactors[i]->getGbObject3D();
- std::array<double, 2> minMax1;
- std::array<double, 2> minMax2;
- std::array<double, 2> minMax3;
- minMax1[0] = geoObject->getX1Minimum();
- minMax2[0] = geoObject->getX2Minimum();
- minMax3[0] = geoObject->getX3Minimum();
- minMax1[1] = geoObject->getX1Maximum();
- minMax2[1] = geoObject->getX2Maximum();
- minMax3[1] = geoObject->getX3Maximum();
-
- for (int x3 = 0; x3 < 2; x3++)
- for (int x2 = 0; x2 < 2; x2++)
- for (int x1 = 0; x1 < 2; x1++) {
- result =
- result || (minMax1[x1] >= AABB[0] && minMax2[x2] >= AABB[1] && minMax3[x3] >= AABB[2] &&
- minMax1[x1] <= AABB[3] && minMax2[x2] <= AABB[4] && minMax3[x3] <= AABB[5]);
- }
- }
- }
-
- std::vector<int> values;
- values.push_back((int)result);
- std::vector<int> rvalues = comm->gather(values);
-
- if (comm->isRoot()) {
- for (int i = 0; i < (int)rvalues.size(); i++) {
- result = result || (bool)rvalues[i];
- }
- }
- int iresult = (int)result;
- comm->broadcast(iresult);
- result = (bool)iresult;
-
- return result;
-}
-
-int DemCoProcessor::addSurfaceTriangleSet(std::vector<UbTupleFloat3> &nodes, std::vector<UbTupleInt3> &triangles)
-{
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- interactors[i]->getGbObject3D()->addSurfaceTriangleSet(nodes, triangles);
- }
- }
- return (int)interactors.size();
-}
-
-void DemCoProcessor::getObjectsPropertiesVector(std::vector<double> &p)
-{
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- p.push_back(i);
- p.push_back(interactors[i]->getGbObject3D()->getX1Centroid());
- p.push_back(interactors[i]->getGbObject3D()->getX2Centroid());
- p.push_back(interactors[i]->getGbObject3D()->getX3Centroid());
- Vector3D v = physicsEngineGeometrieAdapters[i]->getLinearVelocity();
- p.push_back(v[0]);
- p.push_back(v[1]);
- p.push_back(v[2]);
- }
- }
-}
-
-void DemCoProcessor::addPeGeo(walberla::pe::RigidBody *peGeo)
-{
- auto geometry = getPeGeoAdapter(peGeo->getSystemID());
- if (geometry != nullptr) {
- geometry->setActive();
- geometry->setGeometry(peGeo);
- return;
- } else
- return;
-}
-
-void DemCoProcessor::removePeGeo(walberla::pe::RigidBody *peGeo)
-{
- auto geometry = getPeGeoAdapter(peGeo->getSystemID());
- if (geometry != nullptr) {
- geometry->setSemiactive(true);
- } else
- throw UbException(UB_EXARGS, "PeGeo SystemId=" + UbSystem::toString(peGeo->getSystemID()) +
- " is not matching geometry ID");
-}
-
-void DemCoProcessor::addPeShadowGeo(walberla::pe::RigidBody *peGeo)
-{
- auto geometry = getPeGeoAdapter(peGeo->getSystemID());
- if (geometry != nullptr) {
- geometry->setActive();
- geometry->setGeometry(peGeo);
- return;
- } else
- throw UbException(UB_EXARGS,
- "PeGeo ID=" + UbSystem::toString(peGeo->getSystemID()) + " is not matching geometry ID");
-}
-
-void DemCoProcessor::removePeShadowGeo(walberla::pe::RigidBody *peGeo)
-{
- auto geometry = getPeGeoAdapter(peGeo->getSystemID());
-
- if (geometry != nullptr) {
- geometry->setSemiactive(true);
- } else
- throw UbException(UB_EXARGS,
- "PeGeo ID=" + UbSystem::toString(peGeo->getSystemID()) + " is not matching geometry ID");
-}
-
-bool DemCoProcessor::isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d)
-{
- bool result = false;
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- SPtr<GbObject3D> sphere = interactors[i]->getGbObject3D();
- result = result ||
- (sqrt(pow(sphere->getX1Centroid() - centerX1, 2.0) + pow(sphere->getX2Centroid() - centerX2, 2.0) +
- pow(sphere->getX3Centroid() - centerX3, 2.0)) <= d);
- }
- }
- std::vector<int> values;
- values.push_back((int)result);
- std::vector<int> rvalues = comm->gather(values);
-
- if (comm->isRoot()) {
- for (int i = 0; i < (int)rvalues.size(); i++) {
- result = result || (bool)rvalues[i];
- }
- }
- int iresult = (int)result;
- comm->broadcast(iresult);
- result = (bool)iresult;
-
- return result;
-}
-
-void DemCoProcessor::distributeIDs()
-{
- std::vector<unsigned long long> peIDsSend;
- std::vector<int> vfIDsSend;
-
- for (int i = 0; i < interactors.size(); i++) {
- if (std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])->isActive()) {
- peIDsSend.push_back(
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->getSystemID());
- vfIDsSend.push_back(interactors[i]->getID());
- }
- }
-
- std::vector<unsigned long long> peIDsRecv;
- std::vector<int> vfIDsRecv;
-
- comm->allGather(peIDsSend, peIDsRecv);
- comm->allGather(vfIDsSend, vfIDsRecv);
-
- std::map<int, unsigned long long> idMap;
-
- for (int i = 0; i < peIDsRecv.size(); i++) {
- idMap.insert(std::make_pair(vfIDsRecv[i], peIDsRecv[i]));
- }
-
- for (int i = 0; i < interactors.size(); i++) {
- std::map<int, unsigned long long>::const_iterator it;
- if ((it = idMap.find(interactors[i]->getID())) == idMap.end()) {
- throw UbException(UB_EXARGS, "Interactor ID = " + UbSystem::toString(interactors[i]->getID()) +
- " is invalid! The DEM object may be not in PE domain!");
- }
-
- std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])
- ->setSystemID(it->second);
-
- geoIdMap.insert(std::make_pair(
- it->second, std::dynamic_pointer_cast<PePhysicsEngineGeometryAdapter>(physicsEngineGeometrieAdapters[i])));
- }
-}
-//////////////////////////////////////////////////////////////////////////
-void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
-{
- this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
-}
-//////////////////////////////////////////////////////////////////////////
-walberla::pe::RigidBody *DemCoProcessor::getPeGeoObject(walberla::id_t id)
-{
- std::shared_ptr<walberla::blockforest::BlockForest> forest =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getBlockForest();
- std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getStorageId();
- std::shared_ptr<walberla::pe::BodyStorage> globalBodyStorage =
- std::dynamic_pointer_cast<PePhysicsEngineSolverAdapter>(physicsEngineSolver)->getGlobalBodyStorage();
-
- return walberla::pe::getBody(*globalBodyStorage, *forest, *storageId, id,
- walberla::pe::StorageSelect::LOCAL | walberla::pe::StorageSelect::SHADOW);
-}
-////////////////////////////////////////////////////////////////////////////
-std::shared_ptr<PePhysicsEngineGeometryAdapter> DemCoProcessor::getPeGeoAdapter(unsigned long long systemId)
-{
- std::map<unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter>>::const_iterator it;
- if ((it = geoIdMap.find(systemId)) == geoIdMap.end()) {
- return nullptr;
- } else
- return it->second;
-}
diff --git a/src/cpu/DemCoupling/DemCoProcessor.h b/src/cpu/DemCoupling/DemCoProcessor.h
deleted file mode 100644
index d2946f1e93fcaedc69d44a83a68dc2079910e48f..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/DemCoProcessor.h
+++ /dev/null
@@ -1,94 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef DEM_CO_PROCESSOR_H
-#define DEM_CO_PROCESSOR_H
-
-#include <map>
-#include <memory>
-#include <vector>
-
-#include "Vector3D.h"
-
-#include "CoProcessor.h"
-#include "UbTuple.h"
-
-#include <pe/basic.h>
-
-//#define TIMING
-
-#ifdef TIMING
-#include "UbTiming.h"
-#endif
-
-class PhysicsEngineGeometryAdapter;
-class PhysicsEngineSolverAdapter;
-class PePhysicsEngineSolverAdapter;
-class PhysicsEngineMaterialAdapter;
-class PePhysicsEngineGeometryAdapter;
-
-class UbScheduler;
-class Grid3D;
-class ForceCalculator;
-namespace vf::mpi {class Communicator;}
-class MovableObjectInteractor;
-class BoundaryConditionsBlockVisitor;
-
-class DemCoProcessor : public CoProcessor
-{
-public:
- DemCoProcessor(std::shared_ptr<Grid3D> grid, std::shared_ptr<UbScheduler> s, std::shared_ptr<vf::mpi::Communicator> comm,
- std::shared_ptr<ForceCalculator> forceCalculator,
- std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver, double intermediatePeSteps = 1.0);
- virtual ~DemCoProcessor();
-
- void addInteractor(std::shared_ptr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial,
- Vector3D initalVelocity = Vector3D(0.0, 0.0, 0.0));
- void process(double step) override;
- std::shared_ptr<PhysicsEngineSolverAdapter> getPhysicsEngineSolver();
- void distributeIDs();
- void setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> blockVisitor);
- bool isDemObjectInAABB(std::array<double, 6> AABB);
- int addSurfaceTriangleSet(std::vector<UbTupleFloat3> &nodes, std::vector<UbTupleInt3> &triangles);
- void getObjectsPropertiesVector(std::vector<double> &p);
- void addPeGeo(walberla::pe::RigidBody *peGeo);
- void removePeGeo(walberla::pe::RigidBody *peGeo);
- void addPeShadowGeo(walberla::pe::RigidBody *peGeo);
- void removePeShadowGeo(walberla::pe::RigidBody *peGeo);
- bool isSpheresIntersection(double centerX1, double centerX2, double centerX3, double d);
-
-private:
- std::shared_ptr<PhysicsEngineGeometryAdapter>
- createPhysicsEngineGeometryAdapter(std::shared_ptr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial) const;
- void applyForcesOnGeometries();
- void setForcesToObject(SPtr<Grid3D> grid, std::shared_ptr<MovableObjectInteractor> interactor,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry);
- void scaleForcesAndTorques(double scalingFactor);
- void calculateDemTimeStep(double step);
- void moveVfGeoObjects();
- walberla::pe::RigidBody *getPeGeoObject(walberla::id_t id);
- std::shared_ptr<PePhysicsEngineGeometryAdapter> getPeGeoAdapter(unsigned long long systemId);
-
-private:
- std::shared_ptr<vf::mpi::Communicator> comm;
- std::vector<std::shared_ptr<MovableObjectInteractor>> interactors;
- std::shared_ptr<ForceCalculator> forceCalculator;
- std::shared_ptr<PePhysicsEngineSolverAdapter> physicsEngineSolver;
- std::vector<std::shared_ptr<PhysicsEngineGeometryAdapter>> physicsEngineGeometrieAdapters;
- double intermediateDemSteps;
- SPtr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
- // walberla::pe::BodyStorage* bodyStorage; //!< Reference to the central body storage.
- // walberla::pe::BodyStorage* bodyStorageShadowCopies; //!< Reference to the body storage containing body shadow
- // copies.
-
- std::map<unsigned long long, std::shared_ptr<PePhysicsEngineGeometryAdapter>> geoIdMap;
-
-#ifdef TIMING
- UbTimer timer;
-#endif
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/DemCoupling.cmake b/src/cpu/DemCoupling/DemCoupling.cmake
deleted file mode 100644
index 927c08b6dadae76d2ed023253503f8e7bd804601..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/DemCoupling.cmake
+++ /dev/null
@@ -1,31 +0,0 @@
-INCLUDE(${SOURCE_ROOT}/DemCoupling/CMakePackage.txt)
-INCLUDE(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter/CMakePackage.txt)
-INCLUDE(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter/dummy/CMakePackage.txt)
-INCLUDE(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter/pe/CMakePackage.txt)
-INCLUDE(${SOURCE_ROOT}/DemCoupling/reconstructor/CMakePackage.txt)
-
-INCLUDE(${SOURCE_ROOT}/DemCoupling/IncludsList.cmake)
-
-SET(LINK_LIBRARY optimized ${PE_RELEASE_LIBRARY} debug ${PE_DEBUG_LIBRARY})
-SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
-SET(LINK_LIBRARY optimized ${BLOCKFOREST_RELEASE_LIBRARY} debug ${BLOCKFOREST_DEBUG_LIBRARY})
-SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
-SET(LINK_LIBRARY optimized ${DOMAIN_DECOMPOSITION_RELEASE_LIBRARY} debug ${DOMAIN_DECOMPOSITION_DEBUG_LIBRARY})
-SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
-SET(LINK_LIBRARY optimized ${GEOMETRY_RELEASE_LIBRARY} debug ${GEOMETRY_DEBUG_LIBRARY})
-SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
-SET(LINK_LIBRARY optimized ${CORE_RELEASE_LIBRARY} debug ${CORE_DEBUG_LIBRARY})
-SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-
-IF(${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU")
- SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} "stdc++fs")
-ENDIF()
-
-IF(${USE_METIS})
- SET(LINK_LIBRARY optimized ${METIS_RELEASE_LIBRARY} debug ${METIS_DEBUG_LIBRARY})
- SET(CAB_ADDITIONAL_LINK_LIBRARIES ${CAB_ADDITIONAL_LINK_LIBRARIES} ${LINK_LIBRARY})
-ENDIF()
\ No newline at end of file
diff --git a/src/cpu/DemCoupling/IncludsList.cmake b/src/cpu/DemCoupling/IncludsList.cmake
deleted file mode 100644
index 7ebf198e6082131956d5c1e146031394f39e37d5..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/IncludsList.cmake
+++ /dev/null
@@ -1,8 +0,0 @@
-INCLUDE_DIRECTORIES(${SOURCE_ROOT}/DemCoupling)
-INCLUDE_DIRECTORIES(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter)
-INCLUDE_DIRECTORIES(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter/dummy)
-INCLUDE_DIRECTORIES(${SOURCE_ROOT}/DemCoupling/physicsEngineAdapter/pe)
-INCLUDE_DIRECTORIES(${SOURCE_ROOT}/DemCoupling/reconstructor)
-
-INCLUDE_DIRECTORIES(${PE_ROOT}/src)
-INCLUDE_DIRECTORIES(${PE_BINARY_DIR}/src)
\ No newline at end of file
diff --git a/src/cpu/DemCoupling/MovableObjectInteractor.cpp b/src/cpu/DemCoupling/MovableObjectInteractor.cpp
deleted file mode 100644
index 17185c8bb1cdfedfca1d76fa799cc0810a3fb43d..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/MovableObjectInteractor.cpp
+++ /dev/null
@@ -1,254 +0,0 @@
-#include "MovableObjectInteractor.h"
-
-#include "GbObject3D.h"
-#include "UbLogger.h"
-#include "Vector3D.h"
-
-#include "BCAdapter.h"
-#include "BCArray3D.h"
-#include "BCProcessor.h"
-#include "Block3D.h"
-#include "CoordinateTransformation3D.h"
-#include "Grid3D.h"
-#include "ILBMKernel.h"
-
-#include "BoundaryConditionsBlockVisitor.h"
-#include "SetBcBlocksBlockVisitor.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-#include "Reconstructor.h"
-
-#include <array>
-
-//#define TIMING
-
-#ifdef TIMING
-#include "UbTiming.h"
-#endif
-
-MovableObjectInteractor::MovableObjectInteractor(std::shared_ptr<GbObject3D> geoObject3D, std::shared_ptr<Grid3D> grid,
- std::shared_ptr<BCAdapter> bcAdapter, int type,
- std::shared_ptr<Reconstructor> reconstructor, State state)
- : D3Q27Interactor(geoObject3D, grid, bcAdapter, type), reconstructor(reconstructor), state(state)
-{
- // grid->getBlocks(0, grid->getRank(), true, blockVector);
-}
-
-MovableObjectInteractor::~MovableObjectInteractor() {}
-
-void MovableObjectInteractor::setPhysicsEngineGeometry(
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry)
-{
- this->physicsEngineGeometry = physicsEngineGeometry;
- physicsEngineGeometry->changeState(this->state);
-}
-
-void MovableObjectInteractor::moveGbObjectTo(const Vector3D &position)
-{
- // UBLOG(logINFO, "new position: (x,y,z) " << val<1>(position) << ", " << val<2>(position) << ", " <<
- // val<3>(position));
-
- this->getGbObject3D()->setCenterCoordinates(UbTupleDouble3(position[0], position[1], position[2]));
- this->rearrangeGrid();
-}
-
-void MovableObjectInteractor::rearrangeGrid()
-{
-#ifdef TIMING
- UbTimer timer;
- timer.resetAndStart();
-#endif
-
-#ifdef TIMING
- UBLOG(logINFO, "MovableObjectInteractor::rearrangeGrid():start");
-#endif
-
- this->reconstructDistributionOnSolidNodes();
-
-#ifdef TIMING
- UBLOG(logINFO, "reconstructDistributionOnSolidNodes() time = " << timer.stop() << " s");
-#endif
-
- this->setSolidNodesToFluid();
-
-#ifdef TIMING
- UBLOG(logINFO, "setSolidNodesToFluid() time = " << timer.stop() << " s");
-#endif
-
- this->setBcNodesToFluid();
-
-#ifdef TIMING
- UBLOG(logINFO, "setBcNodesToFluid() time = " << timer.stop() << " s");
-#endif
-
- this->removeSolidBlocks();
-
-#ifdef TIMING
- UBLOG(logINFO, "removeSolidBlocks() time = " << timer.stop() << " s");
-#endif
-
- this->removeBcBlocks();
-
-#ifdef TIMING
- UBLOG(logINFO, "removeBcBlocks() time = " << timer.stop() << " s");
-#endif
-
- this->setBcBlocks();
-
-#ifdef TIMING
- UBLOG(logINFO, "setBcBlocks() time = " << timer.stop() << " s");
-#endif
-
- this->initInteractor();
-
-#ifdef TIMING
- UBLOG(logINFO, "initInteractor() time = " << timer.stop() << " s");
-#endif
-
- this->updateVelocityBc();
-
-#ifdef TIMING
- UBLOG(logINFO, "updateVelocityBc() time = " << timer.stop() << " s");
-#endif
-}
-
-void MovableObjectInteractor::updateNodeLists()
-{
- // for (BcNodeIndicesMap::value_type t : bcNodeIndicesMap)
- //{
- // SPtr<Block3D> block = t.first;
- // std::set< UbTupleInt3 >& bcNodeIndices = t.second;
-
- // SPtr<ILBMKernel> kernel = block->getKernel();
-
- // for (UbTupleInt3 node : bcNodeIndices)
- // {
-
- // }
- //}
-}
-
-void MovableObjectInteractor::reconstructDistributionOnSolidNodes()
-{
- for (SolidNodeIndicesMap::value_type t : solidNodeIndicesMap) {
- SPtr<Block3D> block = t.first;
- std::set<UbTupleInt3> &solidNodeIndices = t.second;
-
- SPtr<ILBMKernel> kernel = block->getKernel();
-
- for (UbTupleInt3 node : solidNodeIndices) {
- const int x1 = val<1>(node);
- const int x2 = val<2>(node);
- const int x3 = val<3>(node);
-
- const Vector3D worldCoordinates = this->grid.lock()->getNodeCoordinates(block, x1, x2, x3);
-
- if (kernel->isInsideOfDomain(x1, x2, x3))
- reconstructor->reconstructNode(x1, x2, x3, worldCoordinates, physicsEngineGeometry, kernel);
- }
- }
-}
-
-void MovableObjectInteractor::setSolidNodesToFluid()
-{
- for (SolidNodeIndicesMap::value_type t : solidNodeIndicesMap) {
- SPtr<Block3D> block = t.first;
- std::set<UbTupleInt3> &solidNodeIndices = t.second;
-
- SPtr<ILBMKernel> kernel = block->getKernel();
- SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
-
- for (UbTupleInt3 node : solidNodeIndices)
- bcArray->setFluid(val<1>(node), val<2>(node), val<3>(node));
- }
-}
-
-void MovableObjectInteractor::setBcNodesToFluid()
-{
- for (BcNodeIndicesMap::value_type t : bcNodeIndicesMap) {
- SPtr<Block3D> block = t.first;
- std::set<std::vector<int>> &bcNodeIndices = t.second;
-
- SPtr<ILBMKernel> kernel = block->getKernel();
- SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
-
- for (std::vector<int> node : bcNodeIndices)
- bcArray->setFluid(node[0], node[1], node[2]);
- }
-}
-
-void MovableObjectInteractor::setBcBlocks()
-{
- SetBcBlocksBlockVisitor v(shared_from_this());
- this->grid.lock()->accept(v);
-
- //////////////////////////////////////////////////////////////////////////
- // SPtr<GbObject3D> geoObject = this->getGbObject3D();
- // std::array<double, 6> AABB ={
- // geoObject->getX1Minimum(),geoObject->getX2Minimum(),geoObject->getX3Minimum(),geoObject->getX1Maximum(),geoObject->getX2Maximum(),geoObject->getX3Maximum()
- // }; blockVector.clear(); UbTupleInt3 blockNX=grid.lock()->getBlockNX(); double ext = 0.0;
- // grid.lock()->getBlocksByCuboid(AABB[0]-(double)val<1>(blockNX)*ext, AABB[1]-(double)val<2>(blockNX)*ext,
- // AABB[2]-(double)val<3>(blockNX)*ext, AABB[3]+(double)val<1>(blockNX)*ext, AABB[4]+(double)val<2>(blockNX)*ext,
- // AABB[5]+(double)val<3>(blockNX)*ext, blockVector);
-
- // for(std::shared_ptr<Block3D> block : this->blockVector)
- //{
- // if (block->getKernel())
- // {
- // setBCBlock(block);
- // }
- //}
- //////////////////////////////////////////////////////////////////////////
- // SPtr<GbObject3D> geoObject = this->getGbObject3D();
- // std::array <double, 2> minMax1;
- // std::array <double, 2> minMax2;
- // std::array <double, 2> minMax3;
- // minMax1[0] = geoObject->getX1Minimum();
- // minMax2[0] = geoObject->getX2Minimum();
- // minMax3[0] = geoObject->getX3Minimum();
- // minMax1[1] = geoObject->getX1Maximum();
- // minMax2[1] = geoObject->getX2Maximum();
- // minMax3[1] = geoObject->getX3Maximum();
-
- // SPtr<CoordinateTransformation3D> trafo = grid.lock()->getCoordinateTransformator();
-
- // for (int x3 = 0; x3 < 2; x3++)
- // for (int x2 = 0; x2 < 2; x2++)
- // for (int x1 = 0; x1 < 2; x1++)
- // {
- // int ix1 = (int)trafo->transformForwardToX1Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
- // int ix2 = (int)trafo->transformForwardToX2Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
- // int ix3 = (int)trafo->transformForwardToX3Coordinate(minMax1[x1], minMax2[x2], minMax3[x3]);
- // blockVector.push_back(grid.lock()->getBlock(ix1, ix2, ix3, 0));
- // }
- // for(std::shared_ptr<Block3D> block : this->blockVector)
- //{
- // if (block->getKernel())
- // {
- // setBCBlock(block);
- // }
- //}
-}
-
-void MovableObjectInteractor::updateVelocityBc()
-{
- for (BcNodeIndicesMap::value_type t : this->getBcNodeIndicesMap()) {
- SPtr<Block3D> block = t.first;
- std::set<std::vector<int>> &bcNodeIndices = t.second;
-
- SPtr<BCArray3D> bcArray = block->getKernel()->getBCProcessor()->getBCArray();
-
- for (std::vector<int> node : bcNodeIndices)
- setGeometryVelocityToBoundaryCondition(node, block, bcArray);
- }
-}
-
-void MovableObjectInteractor::setGeometryVelocityToBoundaryCondition(std::vector<int> node, SPtr<Block3D> block,
- SPtr<BCArray3D> bcArray) const
-{
- const SPtr<BoundaryConditions> bc = bcArray->getBC(node[0], node[1], node[2]);
- const Vector3D worldCoordinates = this->grid.lock()->getNodeCoordinates(block, node[0], node[1], node[2]);
- const Vector3D velocity = this->physicsEngineGeometry->getVelocityAtPosition(worldCoordinates);
-
- bc->setBoundaryVelocity(velocity);
-}
diff --git a/src/cpu/DemCoupling/MovableObjectInteractor.h b/src/cpu/DemCoupling/MovableObjectInteractor.h
deleted file mode 100644
index e0e4343a066d5fd69ddb3cf68d2106337a51d031..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/MovableObjectInteractor.h
+++ /dev/null
@@ -1,60 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef D3Q27_MOVABLE_OBJECT_INTERACTOR_H
-#define D3Q27_MOVABLE_OBJECT_INTERACTOR_H
-
-#include <memory>
-#include <vector>
-
-#include "D3Q27Interactor.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-#include "Vector3D.h"
-
-class Grid3D;
-class Block3D;
-class BCArray3D;
-class BCAdapter;
-class GbObject3D;
-
-class PhysicsEngineGeometryAdapter;
-class Reconstructor;
-
-class MovableObjectInteractor : public D3Q27Interactor
-{
-public:
- typedef std::map<SPtr<Block3D>, std::set<std::array<int, 3>>> InBcNodeIndicesMap;
- typedef std::map<SPtr<Block3D>, std::set<std::array<int, 3>>> OutBcNodeIndicesMap;
-
-public:
- MovableObjectInteractor(std::shared_ptr<GbObject3D> geoObject3D, std::shared_ptr<Grid3D> grid,
- std::shared_ptr<BCAdapter> bcAdapter, int type,
- std::shared_ptr<Reconstructor> reconstructor, State isPinned);
- virtual ~MovableObjectInteractor();
-
- void setPhysicsEngineGeometry(std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry);
-
- void moveGbObjectTo(const Vector3D &position);
-
-private:
- void rearrangeGrid();
- void updateNodeLists();
- void setSolidNodesToFluid();
- void setBcNodesToFluid();
- void reconstructDistributionOnSolidNodes();
- void setBcBlocks();
-
- void updateVelocityBc();
- void setGeometryVelocityToBoundaryCondition(std::vector<int> node, std::shared_ptr<Block3D> block,
- std::shared_ptr<BCArray3D> bcArray) const;
-
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry;
-
- std::shared_ptr<Reconstructor> reconstructor;
- State state;
- std::vector<std::shared_ptr<Block3D>> blockVector;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/PePartitioningGridVisitor.cpp b/src/cpu/DemCoupling/PePartitioningGridVisitor.cpp
deleted file mode 100644
index 429eaeb8be0d3a601b64199e5e86279f7d05ce8f..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/PePartitioningGridVisitor.cpp
+++ /dev/null
@@ -1,148 +0,0 @@
-#if defined VF_METIS && defined VF_MPI
-
-#include "PePartitioningGridVisitor.h"
-#include "Block3D.h"
-#include <mpi/Communicator.h>
-#include "CoordinateTransformation3D.h"
-#include "Grid3D.h"
-#include "UbLogger.h"
-#include <math.h>
-#include <shared_mutex>
-
-#include "DemCoProcessor.h"
-
-using namespace std;
-
-PePartitioningGridVisitor::PePartitioningGridVisitor(std::shared_ptr<vf::mpi::Communicator> comm, std::shared_ptr<DemCoProcessor> dem)
- : Grid3DVisitor(), comm(comm), dem(dem)
-{
- forest = dynamicPointerCast<PePhysicsEngineSolverAdapter>(dem->getPhysicsEngineSolver())->getForest();
-}
-//////////////////////////////////////////////////////////////////////////
-PePartitioningGridVisitor::~PePartitioningGridVisitor() {}
-//////////////////////////////////////////////////////////////////////////
-void PePartitioningGridVisitor::visit(SPtr<Grid3D> grid)
-{
- UBLOG(logDEBUG1, "PePartitioningGridVisitor::visit() - start");
-
- collectData(grid);
- distributePartitionData(grid);
-
- UBLOG(logDEBUG1, "PePartitioningGridVisitor::visit() - end");
-}
-//////////////////////////////////////////////////////////////////////////
-void PePartitioningGridVisitor::collectData(SPtr<Grid3D> grid)
-{
- // int minInitLevel = grid->getCoarsestInitializedLevel();
- // int maxInitLevel = grid->getFinestInitializedLevel();
-
- walberla::uint_t peRank;
-
- for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
- forest->getProcessRank(peRank, blockIt->getId());
- vector<SPtr<Block3D>> blocks;
- walberla::AABB aabb = blockIt->getAABB();
-
- // getBlocksByCuboid((double)aabb.xMin(), (double)aabb.yMin(), (double)aabb.zMin(), (double)aabb.xMax(),
- // (double)aabb.yMax(), (double)aabb.zMax(), blocks, grid); for (SPtr<Block3D> block : blocks)
- //{
- // ids.push_back(block->getGlobalID());
- // ranks.push_back((int)peRank);
- //}
- SPtr<Block3D> block = getBlockByMinUniform((double)aabb.xMin(), (double)aabb.yMin(), (double)aabb.zMin(), grid);
- if (block) {
- ids.push_back(block->getGlobalID());
- ranks.push_back((int)peRank);
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-// void PePartitioningGridVisitor::getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double
-// maxX2, double maxX3, std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid)
-//{
-// int coarsestLevel = grid->getCoarsestInitializedLevel();
-// int finestLevel = grid->getFinestInitializedLevel();
-//
-// SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
-//
-// //////////////////////////////////////////////////////////////////////////
-// //MINIMALE BLOCK-INDIZES BESTIMMEN
-// //
-// //min:
-// double dMinX1 = trafo->transformForwardToX1Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
-// double dMinX2 = trafo->transformForwardToX2Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
-// double dMinX3 = trafo->transformForwardToX3Coordinate(minX1, minX2, minX3)*(1<<finestLevel);
-//
-// //Achtung, wenn minX1 genau auf grenze zwischen zwei bloecken -> der "kleinere" muss genommen werden,
-// //da beim Transformieren der "groessere" Index rauskommt
-// int iMinX1 = (int)dMinX1; //if (UbMath::zero(dMinX1-iMinX1)) iMinX1-=1;
-// int iMinX2 = (int)dMinX2; //if (UbMath::zero(dMinX2-iMinX2)) iMinX2-=1;
-// int iMinX3 = (int)dMinX3; //if (UbMath::zero(dMinX3-iMinX3)) iMinX3-=1;
-//
-// //max (hier kann die Zusatzabfrage vernachlaessigt werden):
-// int iMaxX1 = (int)(trafo->transformForwardToX1Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
-// int iMaxX2 = (int)(trafo->transformForwardToX2Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
-// int iMaxX3 = (int)(trafo->transformForwardToX3Coordinate(maxX1, maxX2, maxX3)*(1<<finestLevel));
-//
-// SPtr<Block3D> block;
-//
-// //set, um doppelte bloecke zu vermeiden, die u.U. bei periodic auftreten koennen
-// std::set<SPtr<Block3D>> blockset;
-// for (int level=coarsestLevel; level<=finestLevel; level++)
-// {
-// //damit bei negativen werten auch der "kleinere" genommen wird -> floor!
-// int minx1 = (int)std::floor((double)iMinX1/(1<<(finestLevel-level)));
-// int minx2 = (int)std::floor((double)iMinX2/(1<<(finestLevel-level)));
-// int minx3 = (int)std::floor((double)iMinX3/(1<<(finestLevel-level)));
-//
-// int maxx1 = iMaxX1/(1<<(finestLevel-level));
-// int maxx2 = iMaxX2/(1<<(finestLevel-level));
-// int maxx3 = iMaxX3/(1<<(finestLevel-level));
-//
-// for (int ix1=minx1; ix1<maxx1; ix1++)
-// for (int ix2=minx2; ix2<maxx2; ix2++)
-// for (int ix3=minx3; ix3<maxx3; ix3++)
-// if ((block=grid->getBlock(ix1, ix2, ix3, level)))
-// {
-// blockset.insert(block);
-// }
-// }
-//
-// blocks.resize(blockset.size());
-// std::copy(blockset.begin(), blockset.end(), blocks.begin());
-//}
-
-SPtr<Block3D> PePartitioningGridVisitor::getBlockByMinUniform(double minX1, double minX2, double minX3,
- SPtr<Grid3D> grid)
-{
- SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
-
- int ix1 = (int)trafo->transformForwardToX1Coordinate(minX1, minX2, minX3);
- int ix2 = (int)trafo->transformForwardToX2Coordinate(minX1, minX2, minX3);
- int ix3 = (int)trafo->transformForwardToX3Coordinate(minX1, minX2, minX3);
-
- return grid->getBlock(ix1, ix2, ix3, 0);
-}
-
-//////////////////////////////////////////////////////////////////////////
-void PePartitioningGridVisitor::distributePartitionData(SPtr<Grid3D> grid)
-{
- std::vector<int> totalIDs;
- std::vector<int> totalRanks;
-
- assert(ids.size() != 0);
- assert(ranks.size() != 0);
-
- comm->allGather(ids, totalIDs);
- comm->allGather(ranks, totalRanks);
-
- assert(totalIDs.size() == totalRanks.size());
- for (int i = 0; i < totalIDs.size(); i++) {
- SPtr<Block3D> block = grid->getBlock(totalIDs[i]);
- if (block)
- block->setRank(totalRanks[i]);
- }
-}
-//////////////////////////////////////////////////////////////////////////
-
-#endif // VF_METIS
diff --git a/src/cpu/DemCoupling/PePartitioningGridVisitor.h b/src/cpu/DemCoupling/PePartitioningGridVisitor.h
deleted file mode 100644
index cad80c0f4d986c45560c6111e8943226df136d24..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/PePartitioningGridVisitor.h
+++ /dev/null
@@ -1,57 +0,0 @@
-#ifndef PePartitioningGridVisitor_h
-#define PePartitioningGridVisitor_h
-
-#if defined VF_MPI
-
-#include <PointerDefinitions.h>
-#include <vector>
-
-#include "Grid3DVisitor.h"
-
-#include "PePhysicsEngineSolverAdapter.h"
-
-#include <array>
-
-////////////////////////////////////////////////////////////////////////
-//! \brief The class implements domain decomposition with PE library
-//! \author Konstantin Kutscher
-//////////////////////////////////////////////////////////////////////////
-namespace vf::mpi {class Communicator;}
-class Grid3D;
-class Block3D;
-class DemCoProcessor;
-// class walberla::blockforest::BlockForest;
-
-class PePartitioningGridVisitor : public Grid3DVisitor
-{
-public:
- //! This describe different types of decomposition
- enum GraphType { LevelIntersected, LevelBased };
-
-public:
- //! Constructor
- //! \param comm - communicator
-
- PePartitioningGridVisitor(std::shared_ptr<vf::mpi::Communicator> comm, std::shared_ptr<DemCoProcessor> dem);
- virtual ~PePartitioningGridVisitor();
- void visit(SPtr<Grid3D> grid) override;
-
-protected:
- void collectData(SPtr<Grid3D> grid);
- void distributePartitionData(SPtr<Grid3D> grid);
- // void getBlocksByCuboid(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3,
- // std::vector<SPtr<Block3D>>& blocks, SPtr<Grid3D> grid);
- SPtr<Block3D> getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid);
-
-private:
- std::shared_ptr<vf::mpi::Communicator> comm;
- std::shared_ptr<DemCoProcessor> dem;
-
- std::vector<int> ids;
- std::vector<int> ranks;
-
- std::shared_ptr<walberla::blockforest::BlockForest> forest;
-};
-
-#endif // VF_MPI
-#endif
diff --git a/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.cpp b/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.cpp
deleted file mode 100644
index ff6cbe7e5a3e394bac18016507a57308d0f1ecbf..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.cpp
+++ /dev/null
@@ -1,119 +0,0 @@
-#include "RestartDemObjectsCoProcessor.h"
-
-#include <mpi/Communicator.h>
-#include "CreateDemObjectsCoProcessor.h"
-#include "DemCoProcessor.h"
-#include "GbSphere3D.h"
-#include "Grid3D.h"
-#include "UbFileInputBinary.h"
-#include "UbFileOutputBinary.h"
-#include "UbScheduler.h"
-#include "UbSystem.h"
-#include "Vector3D.h"
-
-RestartDemObjectsCoProcessor::RestartDemObjectsCoProcessor() {}
-
-RestartDemObjectsCoProcessor::RestartDemObjectsCoProcessor(
- SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path, SPtr<DemCoProcessor> demCoProcessor,
- SPtr<CreateDemObjectsCoProcessor> createDemObjectsCoProcessor, double radius, std::shared_ptr<vf::mpi::Communicator> comm)
- : CoProcessor(grid, s), path(path), demCoProcessor(demCoProcessor),
- createDemObjectsCoProcessor(createDemObjectsCoProcessor), radius(radius), comm(comm)
-{
-}
-
-void RestartDemObjectsCoProcessor::process(double step)
-{
- if (scheduler->isDue(step)) {
- int istep = static_cast<int>(step);
-
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::write step: " << istep);
-
- write(istep);
- }
-}
-
-void RestartDemObjectsCoProcessor::restart(double step)
-{
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::read step: " << (int)step);
-
- read((int)step);
-}
-
-void RestartDemObjectsCoProcessor::write(int step)
-{
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::write start ");
- std::vector<double> p;
-
- demCoProcessor->getObjectsPropertiesVector(p);
-
- // TODO implement getherv
- std::vector<double> rvalues;
- comm->allGather(p, rvalues);
-
- if (comm->isRoot()) {
- std::map<int, std::vector<double>> infMap;
- int size = (int)rvalues.size();
- for (int i = 0; i < size; i += 7) {
- std::vector<double> infVector(6);
- for (int j = 0; j < 6; j++) {
- infVector[j] = rvalues[i + 1 + j];
- }
- infMap.insert(std::make_pair((int)rvalues[i], infVector));
- }
- std::vector<double> wvalues;
- typedef std::map<int, std::vector<double>>::iterator it_type;
- for (it_type iterator = infMap.begin(); iterator != infMap.end(); iterator++) {
- // iterator->first = key
- // iterator->second = value
- std::vector<double>::iterator it = wvalues.end();
- it = wvalues.insert(it, iterator->second.begin(), iterator->second.end());
- }
- std::string subfolder = "dem_cp_" + UbSystem::toString(step);
- std::string filePath = path + "/dem_cp/" + subfolder + "/dem_cp.bin";
- UbFileOutputBinary fo(filePath);
- fo.writeInteger((int)wvalues.size());
- fo.writeVector<double>(wvalues);
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::write number of objects = " << wvalues.size() / 6);
- }
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::write stop ");
-}
-
-void RestartDemObjectsCoProcessor::read(int step)
-{
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::read start ");
- std::vector<double> p;
-
- if (comm->isRoot()) {
- std::string subfolder = "dem_cp_" + UbSystem::toString(step);
- std::string filePath = path + "/dem_cp/" + subfolder + "/dem_cp.bin";
- UbFileInputBinary fi(filePath);
- int size = fi.readInteger();
- p.resize(size);
- fi.readVector<double>(p);
- }
- comm->broadcast(p);
-
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::read number of objects = " << p.size() / 6);
-
- createDemObjectsCoProcessor->clearGeoObjects();
-
- int size = (int)p.size();
-
- for (int i = 0; i < size; i += 6) {
- SPtr<GbObject3D> sphere(new GbSphere3D(p[i], p[i + 1], p[i + 2], radius));
- createDemObjectsCoProcessor->addGeoObject(sphere, Vector3D(p[i + 3], p[i + 4], p[i + 5]));
- }
-
- createDemObjectsCoProcessor->createGeoObjects();
-
- createDemObjectsCoProcessor->clearGeoObjects();
-
- if (comm->isRoot())
- UBLOG(logINFO, "RestartDemObjectsCoProcessor::read stop ");
-}
diff --git a/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.h b/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.h
deleted file mode 100644
index 5123a2d6e51ece8e96d6623d573141a8c272026f..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/RestartDemObjectsCoProcessor.h
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * Author: K. Kutscher
- * mail: kutscher@irmb.tu-bs.de
- */
-#ifndef RestartDemObjectsCoProcessor_H
-#define RestartDemObjectsCoProcessor_H
-
-#include <PointerDefinitions.h>
-#include <string>
-#include <vector>
-
-#include "CoProcessor.h"
-
-namespace vf::mpi {class Communicator;}
-class Grid3D;
-class UbScheduler;
-class DemCoProcessor;
-class CreateDemObjectsCoProcessor;
-
-class RestartDemObjectsCoProcessor : public CoProcessor
-{
-public:
- RestartDemObjectsCoProcessor();
- RestartDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path,
- SPtr<DemCoProcessor> demCoProcessor,
- SPtr<CreateDemObjectsCoProcessor> createDemObjectsCoProcessor, double radius,
- std::shared_ptr<vf::mpi::Communicator> comm);
- ~RestartDemObjectsCoProcessor() {}
- void process(double step) override;
- void restart(double step);
- void write(int step);
- void read(int step);
-
-private:
- std::string path;
- double radius;
- std::shared_ptr<vf::mpi::Communicator> comm;
- SPtr<DemCoProcessor> demCoProcessor;
- SPtr<CreateDemObjectsCoProcessor> createDemObjectsCoProcessor;
-};
-#endif
diff --git a/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.cpp b/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.cpp
deleted file mode 100644
index 3e22c90cf266fa8593b0036d160d79080a3ad31c..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-#include "WriteDemObjectsCoProcessor.h"
-
-#include "basics/writer/WbWriterVtkXmlASCII.h"
-#include "basics/writer/WbWriterVtkXmlBinary.h"
-
-#include <mpi/Communicator.h>
-#include "DemCoProcessor.h"
-#include "Grid3D.h"
-#include "UbScheduler.h"
-#include "UbSystem.h"
-
-WriteDemObjectsCoProcessor::WriteDemObjectsCoProcessor() {}
-//////////////////////////////////////////////////////////////////////////
-WriteDemObjectsCoProcessor::WriteDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path,
- WbWriter *const writer, SPtr<DemCoProcessor> demCoProcessor,
- std::shared_ptr<vf::mpi::Communicator> comm)
- : CoProcessor(grid, s), path(path), writer(writer), demCoProcessor(demCoProcessor), comm(comm)
-{
-}
-//////////////////////////////////////////////////////////////////////////
-void WriteDemObjectsCoProcessor::process(double step)
-{
- if (scheduler->isDue(step)) {
- std::vector<UbTupleFloat3> nodes;
- std::vector<UbTupleInt3> triangles;
-
- int numObjcts = demCoProcessor->addSurfaceTriangleSet(nodes, triangles);
-
- int istep = static_cast<int>(step);
-
- std::string pfilePath, partPath, subfolder, cfilePath;
-
- subfolder = "dem" + UbSystem::toString(istep);
- pfilePath = path + "/dem/" + subfolder;
- cfilePath = path + "/dem/dem_collection";
- partPath = pfilePath + "/dem" + UbSystem::toString(comm->getProcessID()) + "_" + UbSystem::toString(istep);
-
- std::string partName = writer->writeTriangles(partPath, nodes, triangles);
- size_t found = partName.find_last_of("/");
- std::string piece = partName.substr(found + 1);
- piece = subfolder + "/" + piece;
-
- std::vector<std::string> datanames;
- std::vector<std::string> cellDataNames;
- std::vector<std::string> pieces = comm->gather(piece);
- if (comm->isRoot()) {
- std::string pname =
- WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
- found = pname.find_last_of("/");
- piece = pname.substr(found + 1);
-
- std::vector<std::string> filenames;
- filenames.push_back(piece);
- if (step == CoProcessor::scheduler->getMinBegin()) {
- WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
- } else {
- WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
- }
- UBLOG(logINFO, "WriteDemObjectsCoProcessor number of objects: " << numObjcts);
- UBLOG(logINFO, "WriteDemObjectsCoProcessor step: " << istep);
- }
- }
-}
diff --git a/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.h b/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.h
deleted file mode 100644
index 7fb3b045ccd439d772ef565c2013af32c75a7a2d..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/WriteDemObjectsCoProcessor.h
+++ /dev/null
@@ -1,35 +0,0 @@
-/*
- * Author: K. Kutscher
- * mail: kutscher@irmb.tu-bs.de
- */
-#ifndef WriteDemObjectsCoProcessor_H
-#define WriteDemObjectsCoProcessor_H
-
-#include <PointerDefinitions.h>
-#include <string>
-#include <vector>
-
-#include "CoProcessor.h"
-
-namespace vf::mpi {class Communicator;}
-class Grid3D;
-class UbScheduler;
-class DemCoProcessor;
-class WbWriter;
-
-class WriteDemObjectsCoProcessor : public CoProcessor
-{
-public:
- WriteDemObjectsCoProcessor();
- WriteDemObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path, WbWriter *const writer,
- SPtr<DemCoProcessor> demCoProcessor, std::shared_ptr<vf::mpi::Communicator> comm);
- ~WriteDemObjectsCoProcessor() {}
- void process(double step) override;
-
-private:
- std::string path;
- WbWriter *writer;
- std::shared_ptr<vf::mpi::Communicator> comm;
- SPtr<DemCoProcessor> demCoProcessor;
-};
-#endif
diff --git a/src/cpu/DemCoupling/WritePeBlocksCoProcessor.cpp b/src/cpu/DemCoupling/WritePeBlocksCoProcessor.cpp
deleted file mode 100644
index 401ea91bc7225eea7f871cbc2e92be44d1a5c9d7..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/WritePeBlocksCoProcessor.cpp
+++ /dev/null
@@ -1,78 +0,0 @@
-#include "WritePeBlocksCoProcessor.h"
-
-#include "basics/writer/WbWriterVtkXmlASCII.h"
-
-#include "Block3D.h"
-#include <mpi/Communicator.h>
-#include "D3Q27System.h"
-#include "Grid3D.h"
-#include "UbScheduler.h"
-
-WritePeBlocksCoProcessor::WritePeBlocksCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path,
- WbWriter *const writer, std::shared_ptr<vf::mpi::Communicator> comm,
- SPtr<walberla::blockforest::BlockForest> forest)
- : CoProcessor(grid, s), path(path), writer(writer), comm(comm), forest(forest)
-{
-}
-
-WritePeBlocksCoProcessor::~WritePeBlocksCoProcessor() {}
-
-void WritePeBlocksCoProcessor::process(double step)
-{
- if (scheduler->isDue(step))
- collectData(step);
-}
-
-void WritePeBlocksCoProcessor::collectData(double step)
-{
- if (comm->getProcessID() == comm->getRoot()) {
- int istep = int(step);
- std::vector<std::string> filenames;
- std::vector<UbTupleFloat3> nodes;
- std::vector<UbTupleInt8> cells;
- std::vector<std::string> celldatanames;
-
- celldatanames.push_back("ID");
- celldatanames.push_back("rank");
-
- walberla::uint_t rank = 0;
-
- std::vector<std::vector<double>> celldata(celldatanames.size());
-
- int nr = 0;
-
- for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
- walberla::AABB aabb = blockIt->getAABB();
-
- nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMin(), (float)aabb.zMin()));
- nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMin(), (float)aabb.zMin()));
- nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMax(), (float)aabb.zMin()));
- nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMax(), (float)aabb.zMin()));
- nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMin(), (float)aabb.zMax()));
- nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMin(), (float)aabb.zMax()));
- nodes.push_back(makeUbTuple((float)aabb.xMax(), (float)aabb.yMax(), (float)aabb.zMax()));
- nodes.push_back(makeUbTuple((float)aabb.xMin(), (float)aabb.yMax(), (float)aabb.zMax()));
- cells.push_back(makeUbTuple(nr, nr + 1, nr + 2, nr + 3, nr + 4, nr + 5, nr + 6, nr + 7));
- nr += 8;
-
- // data
- celldata[0].push_back((double)blockIt->getId().getID());
- forest->getProcessRank(rank, blockIt->getId());
- celldata[1].push_back((double)rank);
- }
-
- filenames.push_back(writer->writeOctsWithCellData(
- path + "/peBlocks/peBlocks_" + UbSystem::toString(grid->getRank()) + "_" + UbSystem::toString(istep), nodes,
- cells, celldatanames, celldata));
-
- if (istep == CoProcessor::scheduler->getMinBegin()) {
- WbWriterVtkXmlASCII::getInstance()->writeCollection(path + "/peBlocks/peBlocks_collection", filenames,
- istep, false);
- } else {
- WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(path + "/peBlocks/peBlocks_collection", filenames,
- istep, false);
- }
-
- UBLOG(logINFO, "WritePeBlocksCoProcessor step: " << istep);
- }
-}
\ No newline at end of file
diff --git a/src/cpu/DemCoupling/WritePeBlocksCoProcessor.h b/src/cpu/DemCoupling/WritePeBlocksCoProcessor.h
deleted file mode 100644
index ae27d50b3f0bba867db7ad8cce79f2e5d8fd5681..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/WritePeBlocksCoProcessor.h
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * WritePeBlocksCoProcessor.h
- *
- * Created on: 07.09.2018
- * Author: K. Kutscher
- */
-
-#ifndef WritePeBlocksCoProcessor_H_
-#define WritePeBlocksCoProcessor_H_
-
-#include <PointerDefinitions.h>
-#include <string>
-
-#include "CoProcessor.h"
-
-#include <pe/basic.h>
-
-namespace vf::mpi {class Communicator;}
-class Grid3D;
-class UbScheduler;
-class WbWriter;
-
-class WritePeBlocksCoProcessor : public CoProcessor
-{
-public:
- WritePeBlocksCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path, WbWriter *const writer,
- std::shared_ptr<vf::mpi::Communicator> comm, SPtr<walberla::blockforest::BlockForest> forest);
- virtual ~WritePeBlocksCoProcessor();
-
- void process(double step) override;
-
-protected:
- void collectData(double step);
-
- std::string path;
- WbWriter *writer;
- std::shared_ptr<vf::mpi::Communicator> comm;
- SPtr<walberla::blockforest::BlockForest> forest;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineGeometryAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineGeometryAdapter.h
deleted file mode 100644
index 490c1f979eaba285f1c39f834396acf9646584d9..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineGeometryAdapter.h
+++ /dev/null
@@ -1,40 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-#define PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-
-#include "Vector3D.h"
-
-enum class State { PIN, UNPIN };
-
-class PhysicsEngineGeometryAdapter
-{
-public:
- virtual ~PhysicsEngineGeometryAdapter() {}
-
- virtual void addForce(const Vector3D &force) = 0;
- virtual void addTorque(const Vector3D &torque) = 0;
-
- virtual void setForce(const Vector3D &force) = 0;
- virtual void setTorque(const Vector3D &torque) = 0;
-
- virtual void addForceAtPosition(const Vector3D &force, const Vector3D &position) = 0;
- virtual void setLinearVelolocity(const Vector3D &velocity) = 0;
- virtual void setAngularVelocity(const Vector3D &velocity) = 0;
-
- virtual void resetForceAndTorque() = 0;
-
- virtual Vector3D getPosition() const = 0;
- virtual Vector3D getVelocityAtPosition(const Vector3D &position) const = 0;
- virtual Vector3D getLinearVelocity() const = 0;
- virtual Vector3D getAngularVelocity() const = 0;
-
- virtual Vector3D getForce() const = 0;
- virtual Vector3D getTorque() const = 0;
-
- virtual void changeState(State state) = 0;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineMaterialAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineMaterialAdapter.h
deleted file mode 100644
index 30504bee98580f25f23e37030d5ec8180fce3ebc..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineMaterialAdapter.h
+++ /dev/null
@@ -1,39 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PHYSICS_ENGINE_MATERIAL_ADAPTER_H
-#define PHYSICS_ENGINE_MATERIAL_ADAPTER_H
-
-#include <string>
-
-class PhysicsEngineMaterialAdapter
-{
-public:
- PhysicsEngineMaterialAdapter(std::string name, double density, double restitution, double staticFriction,
- double dynamicFriction, double poissonRatio, double youngModul,
- double stiffnessInNormalDirection, double dampingoefficientNormalDirection,
- double dampingTangentialDirection)
- : name(name), density(density), restitution(restitution), staticFriction(staticFriction),
- dynamicFriction(dynamicFriction), poissonRatio(poissonRatio), youngModul(youngModul),
- stiffnessInNormalDirection(stiffnessInNormalDirection),
- dampingoefficientNormalDirection(dampingoefficientNormalDirection),
- dampingTangentialDirection(dampingTangentialDirection)
- {
- }
- virtual ~PhysicsEngineMaterialAdapter() {}
-
-protected:
- std::string name;
- double density;
- double restitution;
- double staticFriction; // Note: pe doubles the input coefficient of friction for material-material contacts.
- double dynamicFriction; // Similar to static friction for low speed friction.
- double poissonRatio;
- double youngModul;
- double stiffnessInNormalDirection;
- double dampingoefficientNormalDirection;
- double dampingTangentialDirection;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineSolverAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineSolverAdapter.h
deleted file mode 100644
index 8e03bf1d6e651f993012acc1e8297e309e993cc8..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/PhysicsEngineSolverAdapter.h
+++ /dev/null
@@ -1,24 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PHYSICS_ENGINE_SOLVER_ADAPTER_H
-#define PHYSICS_ENGINE_SOLVER_ADAPTER_H
-
-#include "Vector3D.h"
-
-class PhysicsEngineGeometryAdapter;
-class PhysicsEngineMaterialAdapter;
-
-class PhysicsEngineSolverAdapter
-{
-public:
- virtual ~PhysicsEngineSolverAdapter() {}
-
- virtual std::shared_ptr<PhysicsEngineGeometryAdapter>
- createPhysicsEngineGeometryAdapter(int id, const Vector3D &position, double radius,
- std::shared_ptr<PhysicsEngineMaterialAdapter> material) const = 0;
- virtual void runTimestep(double step) = 0;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.cpp
deleted file mode 100644
index b18a57532880491a2419ff9d78bc0c64aac108f8..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#include "DummyPhysicsEngineGeometryAdapter.h"
-
-void DummyPhysicsEngineGeometryAdapter::addForce(const Vector3D &force) {}
-
-void DummyPhysicsEngineGeometryAdapter::addTorque(const Vector3D &torque) {}
-
-void DummyPhysicsEngineGeometryAdapter::setForce(const Vector3D &force) {}
-
-void DummyPhysicsEngineGeometryAdapter::setTorque(const Vector3D &torque) {}
-
-void DummyPhysicsEngineGeometryAdapter::addForceAtPosition(const Vector3D &force, const Vector3D &position) {}
-
-void DummyPhysicsEngineGeometryAdapter::setLinearVelolocity(const Vector3D &velocity) { this->velocity = velocity; }
-
-void DummyPhysicsEngineGeometryAdapter::setAngularVelocity(const Vector3D &velocity) {}
-
-void DummyPhysicsEngineGeometryAdapter::resetForceAndTorque() {}
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getVelocityAtPosition(const Vector3D &position) const { return velocity; }
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getLinearVelocity() const { return Vector3D(); }
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getAngularVelocity() const { return Vector3D(); }
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getPosition() const { return Vector3D(); }
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getForce() const { return Vector3D(); }
-
-Vector3D DummyPhysicsEngineGeometryAdapter::getTorque() const { return Vector3D(); }
-
-void DummyPhysicsEngineGeometryAdapter::changeState(State state) {}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.h
deleted file mode 100644
index 70620d3c0681a2d2dc702c4a74f5c8d5a94141ea..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineGeometryAdapter.h
+++ /dev/null
@@ -1,43 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef DUMMY_PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-#define DUMMY_PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-
-#include "UbTuple.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-
-class DummyPhysicsEngineGeometryAdapter : public PhysicsEngineGeometryAdapter
-{
-public:
- DummyPhysicsEngineGeometryAdapter() {}
- virtual ~DummyPhysicsEngineGeometryAdapter() {}
-
- void addForce(const Vector3D &force) override;
- void addTorque(const Vector3D &torque) override;
-
- void setForce(const Vector3D &force) override;
- void setTorque(const Vector3D &torque) override;
-
- void addForceAtPosition(const Vector3D &force, const Vector3D &position) override;
- void setLinearVelolocity(const Vector3D &velocity) override;
- void setAngularVelocity(const Vector3D &velocity) override;
-
- void resetForceAndTorque() override;
-
- Vector3D getVelocityAtPosition(const Vector3D &position) const override;
- Vector3D getLinearVelocity() const override;
- Vector3D getAngularVelocity() const override;
- Vector3D getPosition() const override;
- Vector3D getForce() const override;
- Vector3D getTorque() const override;
-
- void changeState(State state) override;
-
-private:
- Vector3D velocity;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.cpp
deleted file mode 100644
index 7890f966872d5433efa5b40b6358ca0b5a25a40d..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.cpp
+++ /dev/null
@@ -1 +0,0 @@
-#include "DummyPhysicsEngineMaterialAdapter.h"
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.h
deleted file mode 100644
index e84e0f1089a017ed9135383995c96e55cdf0cee6..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineMaterialAdapter.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef DUMMY_PHYSICS_ENGINE_MATERIAL_ADAPTER
-#define DUMMY_PHYSICS_ENGINE_MATERIAL_ADAPTER
-
-#include "PhysicsEngineMaterialAdapter.h"
-
-class DummyPhysicsEngineMaterialAdapter : public PhysicsEngineMaterialAdapter
-{
-public:
- DummyPhysicsEngineMaterialAdapter(std::string name, double density, double restitution, double staticFriction,
- double dynamicFriction, double poissonRatio, double youngModul,
- double stiffnessInNormalDirection, double dampingoefficientNormalDirection,
- double dampingTangentialDirection)
- : PhysicsEngineMaterialAdapter(name, density, restitution, staticFriction, dynamicFriction, poissonRatio,
- youngModul, stiffnessInNormalDirection, dampingoefficientNormalDirection,
- dampingTangentialDirection)
- {
- }
- virtual ~DummyPhysicsEngineMaterialAdapter() {}
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.cpp
deleted file mode 100644
index 321b523ead35a6a8d16d7b48b1fbb7cb66b0fdc8..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.cpp
+++ /dev/null
@@ -1,12 +0,0 @@
-#include "DummyPhysicsEngineSolverAdapter.h"
-
-#include "DummyPhysicsEngineGeometryAdapter.h"
-
-std::shared_ptr<PhysicsEngineGeometryAdapter> DummyPhysicsEngineSolverAdapter::createPhysicsEngineGeometryAdapter(
- int id, const Vector3D &position, double radius, std::shared_ptr<PhysicsEngineMaterialAdapter> material) const
-{
- return std::static_pointer_cast<PhysicsEngineGeometryAdapter>(
- std::make_shared<DummyPhysicsEngineGeometryAdapter>());
-}
-
-void DummyPhysicsEngineSolverAdapter::runTimestep(double step) {}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.h
deleted file mode 100644
index 38e9e7f055b415266cfd35c055951707ba4cda44..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/dummy/DummyPhysicsEngineSolverAdapter.h
+++ /dev/null
@@ -1,26 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef DUMMY_PHYSICS_ENGINE_SOLVER_ADAPTER_H
-#define DUMMY_PHYSICS_ENGINE_SOLVER_ADAPTER_H
-
-#include <memory>
-
-#include "UbTuple.h"
-
-#include "PhysicsEngineSolverAdapter.h"
-
-class DummyPhysicsEngineSolverAdapter : public PhysicsEngineSolverAdapter
-{
-public:
- DummyPhysicsEngineSolverAdapter(){};
- virtual ~DummyPhysicsEngineSolverAdapter() {}
-
- std::shared_ptr<PhysicsEngineGeometryAdapter>
- createPhysicsEngineGeometryAdapter(int id, const Vector3D &position, double radius,
- std::shared_ptr<PhysicsEngineMaterialAdapter> material) const override;
- void runTimestep(double step) override;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapter.h
deleted file mode 100644
index 0373281d60b6618555de3d48190816d79b9ade5b..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapter.h
+++ /dev/null
@@ -1,19 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PE_ADAPTER_H
-#define PE_ADAPTER_H
-
-#include "Vector3D.h"
-#include <pe/basic.h>
-
-class PeConverter
-{
-public:
- static Vector3D convert(walberla::pe::Vec3 vec3) { return Vector3D(vec3[0], vec3[1], vec3[2]); }
-
- static walberla::pe::Vec3 convert(const Vector3D &vec3) { return walberla::pe::Vec3(vec3[0], vec3[1], vec3[2]); }
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapterTest.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapterTest.cpp
deleted file mode 100644
index a92127c3fd0e2b10ad2d73ee88e7fbd00b8ace85..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeAdapterTest.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-//#include "gmock/gmock.h"
-//
-//#include "PeAdapter.h"
-//#include <pe/basic.h>
-//
-//#include "UbTuple.h"
-//
-//
-// TEST(PeAdapterTest, convert_WalberlaVec3_to_Vector3D)
-//{
-// walberla::pe::Vec3 walberlaVec(1.0, -2.0, 3.4);
-// Vector3D ubTuple = PeConverter::convert(walberlaVec);
-//
-// EXPECT_THAT(ubTuple[0], testing::DoubleEq(walberlaVec[0]));
-// EXPECT_THAT(ubTuple[1], testing::DoubleEq(walberlaVec[1]));
-// EXPECT_THAT(ubTuple[2], testing::DoubleEq(walberlaVec[2]));
-//}
-//
-// TEST(PeAdapterTest, convert_Vector3D_to_WalberlaVec3)
-//{
-// Vector3D ubTuple(1.0, -2.0, 3.4);
-// walberla::pe::Vec3 walberlaVec = PeConverter::convert(ubTuple);
-//
-// EXPECT_THAT(ubTuple[0], testing::DoubleEq(walberlaVec[0]));
-// EXPECT_THAT(ubTuple[1], testing::DoubleEq(walberlaVec[1]));
-// EXPECT_THAT(ubTuple[2], testing::DoubleEq(walberlaVec[2]));
-//}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
deleted file mode 100644
index 6c597698d608a287cc3a8bd5db84fbf061539234..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.cpp
+++ /dev/null
@@ -1,47 +0,0 @@
-#include "PeLoadBalancerAdapter.h"
-#include "Block3D.h"
-#include "CoordinateTransformation3D.h"
-#include "Grid3D.h"
-#include "UbLogger.h"
-
-#include "core/debug/CheckFunctions.h"
-
-PeLoadBalancerAdapter::PeLoadBalancerAdapter(SPtr<Grid3D> grid, unsigned numberOfProcesses, int rank)
- : grid(grid), numberOfProcesses(numberOfProcesses), rank(rank)
-{
-}
-
-walberla::uint_t PeLoadBalancerAdapter::operator()(walberla::SetupBlockForest &forest,
- const walberla::uint_t numberOfProcesses,
- const walberla::memory_t perProcessMemoryLimit)
-{
- std::vector<walberla::SetupBlock *> peBlocks;
- forest.getBlocks(peBlocks);
-
- for (auto peBlock = peBlocks.begin(); peBlock != peBlocks.end(); ++peBlock) {
- walberla::AABB aabb = (*peBlock)->getAABB();
- SPtr<Block3D> block = getBlockByMinUniform(aabb.xMin() + 0.5 * (aabb.xMax() - aabb.xMin()),
- aabb.yMin() + 0.5 * (aabb.yMax() - aabb.yMin()),
- aabb.zMin() + 0.5 * (aabb.zMax() - aabb.zMin()), grid);
- if (block) {
- (*peBlock)->assignTargetProcess((walberla::uint_t)block->getRank());
- } else {
- // TODO: the rank of pe blocks is not consistent with VF blocks
- (*peBlock)->assignTargetProcess(0);
- // UBLOG(logINFO, "PeLoadBalancerAdapter::operator() peBlockId="<<(*peBlock)->getId());
- }
- }
-
- return numberOfProcesses;
-}
-
-SPtr<Block3D> PeLoadBalancerAdapter::getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid)
-{
- SPtr<CoordinateTransformation3D> trafo = grid->getCoordinateTransformator();
-
- int ix1 = (int)trafo->transformForwardToX1Coordinate(minX1, minX2, minX3);
- int ix2 = (int)trafo->transformForwardToX2Coordinate(minX1, minX2, minX3);
- int ix3 = (int)trafo->transformForwardToX3Coordinate(minX1, minX2, minX3);
-
- return grid->getBlock(ix1, ix2, ix3, 0);
-}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h
deleted file mode 100644
index 9e1c64dd330cd9b1aa06857d4d441b736c8a39a1..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PeLoadBalancerAdapter.h
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef PeLoadBalancerAdapter_h__
-#define PeLoadBalancerAdapter_h__
-
-#include "PointerDefinitions.h"
-#include "blockforest/SetupBlockForest.h"
-
-class Grid3D;
-class Block3D;
-
-class PeLoadBalancerAdapter
-{
-public:
- PeLoadBalancerAdapter(SPtr<Grid3D> grid, unsigned numberOfProcesses, int rank);
- walberla::uint_t operator()(walberla::SetupBlockForest &forest, const walberla::uint_t numberOfProcesses,
- const walberla::memory_t perProcessMemoryLimit);
- unsigned getNumberOfProcesses() const { return numberOfProcesses; }
- int getRank() const { return rank; }
-
-protected:
- SPtr<Block3D> getBlockByMinUniform(double minX1, double minX2, double minX3, SPtr<Grid3D> grid);
-
-private:
- SPtr<Grid3D> grid;
- unsigned numberOfProcesses;
- int rank;
-};
-
-#endif // PeLoadBalancerAdapter_h__
\ No newline at end of file
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
deleted file mode 100644
index 9800d75b18a78b7eaf0b46c4193b93a55f3ff91b..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-#include "PePhysicsEngineGeometryAdapter.h"
-
-#include <pe/basic.h>
-
-#include "PeAdapter.h"
-
-// PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject) :
-// peGeoObject(peGeoObject)
-//{
-// this->id = peGeoObject->getID();
-// this->active = true;
-//}
-
-PePhysicsEngineGeometryAdapter::PePhysicsEngineGeometryAdapter()
-{
- this->id = -999;
- this->systemID = -999;
- this->active = false;
- this->semiactive = false;
- shadowCounter = 0;
- counter = 0;
-}
-
-void PePhysicsEngineGeometryAdapter::addForce(const Vector3D &force)
-{
- peGeoObject->addForce(PeConverter::convert(force));
-}
-
-void PePhysicsEngineGeometryAdapter::addTorque(const Vector3D &torque)
-{
- peGeoObject->addTorque(PeConverter::convert(torque));
-}
-
-void PePhysicsEngineGeometryAdapter::setForce(const Vector3D &force)
-{
- peGeoObject->setForce(PeConverter::convert(force));
-}
-
-void PePhysicsEngineGeometryAdapter::setTorque(const Vector3D &torque)
-{
- peGeoObject->setTorque(PeConverter::convert(torque));
-}
-
-void PePhysicsEngineGeometryAdapter::addForceAtPosition(const Vector3D &force, const Vector3D &position)
-{
- peGeoObject->addForceAtPos(PeConverter::convert(force), PeConverter::convert(position));
-}
-
-void PePhysicsEngineGeometryAdapter::setLinearVelolocity(const Vector3D &velocity)
-{
- peGeoObject->setLinearVel(PeConverter::convert(velocity));
-}
-
-void PePhysicsEngineGeometryAdapter::setAngularVelocity(const Vector3D &velocity)
-{
- peGeoObject->setAngularVel(PeConverter::convert(velocity));
-}
-
-void PePhysicsEngineGeometryAdapter::resetForceAndTorque() { peGeoObject->resetForceAndTorque(); }
-
-Vector3D PePhysicsEngineGeometryAdapter::getVelocityAtPosition(const Vector3D &position) const
-{
- return PeConverter::convert(peGeoObject->velFromWF(PeConverter::convert(position)));
-}
-
-Vector3D PePhysicsEngineGeometryAdapter::getLinearVelocity() const
-{
- return PeConverter::convert(peGeoObject->getLinearVel());
-}
-
-Vector3D PePhysicsEngineGeometryAdapter::getAngularVelocity() const
-{
- return PeConverter::convert(peGeoObject->getAngularVel());
-}
-
-Vector3D PePhysicsEngineGeometryAdapter::getPosition() const
-{
- return PeConverter::convert(peGeoObject->getPosition());
-}
-
-Vector3D PePhysicsEngineGeometryAdapter::getForce() const { return PeConverter::convert(peGeoObject->getForce()); }
-
-Vector3D PePhysicsEngineGeometryAdapter::getTorque() const { return PeConverter::convert(peGeoObject->getTorque()); }
-
-void PePhysicsEngineGeometryAdapter::changeState(State state)
-{
- if (state == State::PIN)
- peGeoObject->setMassAndInertiaToInfinity();
-}
-
-int PePhysicsEngineGeometryAdapter::getId() const { return id; }
-
-void PePhysicsEngineGeometryAdapter::setId(int id) { this->id = id; }
-
-void PePhysicsEngineGeometryAdapter::setGeometry(walberla::pe::RigidBody *peGeoObject)
-{
- this->peGeoObject = peGeoObject;
-}
-
-//////////////////////////////////////////////////////////////////////////
-void PePhysicsEngineGeometryAdapter::setActive() { active = true; }
-//////////////////////////////////////////////////////////////////////////
-void PePhysicsEngineGeometryAdapter::setInactive() { active = false; }
-//////////////////////////////////////////////////////////////////////////
-bool PePhysicsEngineGeometryAdapter::isActive() { return active; }
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
deleted file mode 100644
index a8eaa7d33ede840f4d76ae90ffb44bef30139f99..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineGeometryAdapter.h
+++ /dev/null
@@ -1,74 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PE_PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-#define PE_PHYSICS_ENGINE_GEOMETRY_ADAPTER_H
-
-#include "PhysicsEngineGeometryAdapter.h"
-#include <core/DataTypes.h>
-
-namespace walberla
-{
-namespace pe
-{
-class RigidBody;
-}
-} // namespace walberla
-
-class PePhysicsEngineGeometryAdapter : public PhysicsEngineGeometryAdapter
-{
-public:
- PePhysicsEngineGeometryAdapter();
- // PePhysicsEngineGeometryAdapter(walberla::pe::RigidBody* peGeoObject);
- virtual ~PePhysicsEngineGeometryAdapter() {}
-
- void addForce(const Vector3D &force) override;
- void addTorque(const Vector3D &torque) override;
-
- void setForce(const Vector3D &force) override;
- void setTorque(const Vector3D &torque) override;
-
- void addForceAtPosition(const Vector3D &force, const Vector3D &position) override;
- void setLinearVelolocity(const Vector3D &velocity) override;
- void setAngularVelocity(const Vector3D &velocity) override;
-
- void resetForceAndTorque() override;
-
- Vector3D getVelocityAtPosition(const Vector3D &position) const override;
- Vector3D getLinearVelocity() const override;
- Vector3D getAngularVelocity() const override;
- Vector3D getPosition() const override;
- Vector3D getForce() const override;
- Vector3D getTorque() const override;
-
- void changeState(State state) override;
-
- int getId() const;
- void setId(int id);
- void setGeometry(walberla::pe::RigidBody *peGeoObject);
-
- void setActive();
- void setInactive();
- bool isActive();
- // void increaseShadowCounter();
- // void decreaseShad
- int shadowCounter;
- int counter;
-
- unsigned long long getSystemID() const { return systemID; }
- void setSystemID(unsigned long long val) { systemID = val; }
- bool getSemiactive() const { return semiactive; }
- void setSemiactive(bool val) { semiactive = val; }
-
-private:
- walberla::pe::RigidBody *peGeoObject;
- // unsigned long long id;
- int id;
- // walberla::id_t systemId;
- unsigned long long systemID;
- bool active;
- bool semiactive;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.cpp
deleted file mode 100644
index 6a36fa1a3b6415e255e0083ec621f1bcea03e3de..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#include "PePhysicsEngineMaterialAdapter.h"
-
-walberla::pe::MaterialID PePhysicsEngineMaterialAdapter::getPeMaterial() const
-{
- if (walberla::pe::Material::find(name) != -1)
- return walberla::pe::Material::find(name);
-
- return walberla::pe::createMaterial(name, density, restitution, staticFriction, dynamicFriction, poissonRatio,
- youngModul, stiffnessInNormalDirection, dampingoefficientNormalDirection,
- dampingTangentialDirection);
-}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.h
deleted file mode 100644
index 6ebfa8d1d9ca67760bffb9e06cafec256b19cf4f..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineMaterialAdapter.h
+++ /dev/null
@@ -1,28 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PE_PHYSICS_ENGINE_MATERIAL_ADAPTER
-#define PE_PHYSICS_ENGINE_MATERIAL_ADAPTER
-
-#include "../PhysicsEngineMaterialAdapter.h"
-#include <pe/basic.h>
-
-class PePhysicsEngineMaterialAdapter : public PhysicsEngineMaterialAdapter
-{
-public:
- PePhysicsEngineMaterialAdapter(std::string name, double density, double restitution, double staticFriction,
- double dynamicFriction, double poissonRatio, double youngModul,
- double stiffnessInNormalDirection, double dampingoefficientNormalDirection,
- double dampingTangentialDirection)
- : PhysicsEngineMaterialAdapter(name, density, restitution, staticFriction, dynamicFriction, poissonRatio,
- youngModul, stiffnessInNormalDirection, dampingoefficientNormalDirection,
- dampingTangentialDirection)
- {
- }
- virtual ~PePhysicsEngineMaterialAdapter() {}
-
- virtual walberla::pe::MaterialID getPeMaterial() const;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
deleted file mode 100644
index 14cef406392fbfbd9862a71b0c054df85a8608ec..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.cpp
+++ /dev/null
@@ -1,235 +0,0 @@
-#include "PePhysicsEngineSolverAdapter.h"
-
-#include <exception>
-
-#include "pe/rigidbody/BoxFactory.h"
-#include "pe/rigidbody/PlaneFactory.h"
-#include "pe/rigidbody/SphereFactory.h"
-#include <pe/basic.h>
-#include <pe/rigidbody/UnionFactory.h>
-//#include "geometry/GeometricalFunctions.h"
-#include <mpi/Communicator.h>
-#include "PeAdapter.h"
-#include "PeLoadBalancerAdapter.h"
-#include "PePhysicsEngineGeometryAdapter.h"
-#include "PePhysicsEngineMaterialAdapter.h"
-#include "UbException.h"
-#include "UbLogger.h"
-#include "UbSystem.h"
-#include <boost/tuple/tuple.hpp>
-#include <memory>
-
-using namespace walberla;
-using namespace walberla::pe;
-
-typedef boost::tuple<walberla::pe::Box, walberla::pe::Sphere, walberla::pe::Plane> BodyTypeTuple;
-
-PePhysicsEngineSolverAdapter::PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter,
- std::shared_ptr<PeLoadBalancerAdapter> loadBalancer)
- : peParameter(peParameter), loadBalancer(loadBalancer)
-{
- this->initalizePeEnvironment();
-}
-
-void PePhysicsEngineSolverAdapter::initalizePeEnvironment()
-{
- this->initialPeBodyStorage();
- this->initialPeBlockForest();
- this->initalBlockData();
- this->initalPeIntegrator();
- this->executePeBodyTypeTuple();
- this->initialPeChannel();
-}
-
-std::shared_ptr<PhysicsEngineGeometryAdapter> PePhysicsEngineSolverAdapter::createPhysicsEngineGeometryAdapter(
- int id, const Vector3D &position, double radius, std::shared_ptr<PhysicsEngineMaterialAdapter> material) const
-{
- const std::shared_ptr<PePhysicsEngineMaterialAdapter> peMaterial =
- std::dynamic_pointer_cast<PePhysicsEngineMaterialAdapter>(material);
- std::shared_ptr<PePhysicsEngineGeometryAdapter> peGeometryAdapter(new PePhysicsEngineGeometryAdapter());
-
- // UBLOG(logINFO, "PePhysicsEngineSolverAdapter::createSphere():start");
- walberla::pe::GeomID peGeometry = createSphere(*globalBodyStorage, *forest, *storageId, id,
- PeConverter::convert(position), radius, peMaterial->getPeMaterial());
- // UBLOG(logINFO, "PePhysicsEngineSolverAdapter::createSphere():end");
-
- if (peGeometry) {
- peGeometryAdapter->setId(id);
- peGeometryAdapter->setSystemID(peGeometry->getSystemID());
- peGeometryAdapter->setActive();
- peGeometryAdapter->setGeometry(peGeometry);
- return peGeometryAdapter;
- } else {
- peGeometryAdapter->setId(id);
- peGeometryAdapter->setInactive();
- return peGeometryAdapter;
- }
-
- walberla::pe::syncNextNeighbors<BodyTypeTuple>(*forest, *storageId);
-}
-
-void PePhysicsEngineSolverAdapter::runTimestep(double step)
-{
- cr->timestep(walberla::real_c(step));
- walberla::pe::syncNextNeighbors<BodyTypeTuple>(*forest, *storageId);
-}
-
-void PePhysicsEngineSolverAdapter::initialPeBodyStorage()
-{
- globalBodyStorage = std::make_shared<walberla::pe::BodyStorage>();
-}
-
-void PePhysicsEngineSolverAdapter::initialPeBlockForest()
-{
-
- // walberla::SetupBlockForest sforest =
- // walberla::blockforest::createUniformBlockGrid(walberla::AABB(peParameter->simulationDomain[0],
- // peParameter->simulationDomain[1], peParameter->simulationDomain[2],
- // peParameter->simulationDomain[3], peParameter->simulationDomain[4], peParameter->simulationDomain[5]), //
- // simulationDomain walberla::uint_t(val<1>(peParameter->numberOfBlocks)),
- // walberla::uint_t(val<2>(peParameter->numberOfBlocks)),
- // walberla::uint_t(val<3>(peParameter->numberOfBlocks)),walberla::uint_t(10),walberla::uint_t(10),walberla::uint_t(10),
- // 5.0,false);
- walberla::SetupBlockForest sforest;
- // sforest.addWorkloadMemorySUIDAssignmentFunction( uniformWorkloadAndMemoryAssignment );
- sforest.init(walberla::AABB(peParameter->simulationDomain[0], peParameter->simulationDomain[1],
- peParameter->simulationDomain[2], peParameter->simulationDomain[3],
- peParameter->simulationDomain[4], peParameter->simulationDomain[5]), // simulationDomain
- walberla::uint_t(val<1>(peParameter->numberOfBlocks)),
- walberla::uint_t(val<2>(peParameter->numberOfBlocks)),
- walberla::uint_t(val<3>(peParameter->numberOfBlocks)), // blocks in each direction
- val<1>(peParameter->isPeriodic), val<2>(peParameter->isPeriodic), val<3>(peParameter->isPeriodic));
- sforest.balanceLoad(*loadBalancer.get(), loadBalancer->getNumberOfProcesses());
- forest = std::shared_ptr<walberla::blockforest::BlockForest>(
- new walberla::blockforest::BlockForest(walberla::uint_c(loadBalancer->getRank()), sforest));
-
- auto mpiManager = walberla::MPIManager::instance();
- mpiManager->useWorldComm();
- if (!forest)
- throw std::runtime_error("No PE BlockForest created ... ");
-}
-
-void PePhysicsEngineSolverAdapter::initalBlockData()
-{
- storageId = std::make_shared<walberla::domain_decomposition::BlockDataID>(
- forest->addBlockData(walberla::pe::createStorageDataHandling<BodyTypeTuple>(), "Storage"));
-}
-
-void PePhysicsEngineSolverAdapter::initalPeIntegrator()
-{
- auto ccdID =
- forest->addBlockData(walberla::pe::ccd::createHashGridsDataHandling(globalBodyStorage, *storageId), "CCD");
- auto fcdID = forest->addBlockData(
- walberla::pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, walberla::pe::fcd::AnalyticCollideFunctor>(),
- "FCD");
-
- cr = std::make_shared<walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers>(globalBodyStorage, forest,
- *storageId, ccdID, fcdID);
- cr->setMaxIterations(peParameter->maxPeIterations);
- cr->setRelaxationModel(
- walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers::ApproximateInelasticCoulombContactByDecoupling);
- cr->setRelaxationParameter(walberla::real_t(peParameter->relaxationParameter));
- cr->setGlobalLinearAcceleration(PeConverter::convert(peParameter->globalLinearAcceleration));
-}
-
-void PePhysicsEngineSolverAdapter::executePeBodyTypeTuple() { walberla::pe::SetBodyTypeIDs<BodyTypeTuple>::execute(); }
-
-void PePhysicsEngineSolverAdapter::initialPeChannel() const
-{
- const walberla::pe::MaterialID material = peParameter->planes->getPeMaterial();
-
- auto simulationDomain = forest->getDomain();
-
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(1, 0, 0), simulationDomain.minCorner(), material);
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(-1, 0, 0), simulationDomain.maxCorner(), material);
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 1, 0), simulationDomain.minCorner(), material);
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, -1, 0), simulationDomain.maxCorner(), material);
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 0, 1), simulationDomain.minCorner(), material);
- // createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 0, -1), simulationDomain.maxCorner(), material);
-
- Vector3D minOffset = peParameter->minOffset;
- Vector3D maxOffset = peParameter->maxOffset;
-
- walberla::pe::Vec3 minX1_Offset(minOffset.X1(), 0, 0);
- walberla::pe::Vec3 maxX1_Offset(maxOffset.X1(), 0, 0);
- walberla::pe::Vec3 minX2_Offset(0, minOffset.X2(), 0);
- walberla::pe::Vec3 maxX2_Offset(0, maxOffset.X2(), 0);
- walberla::pe::Vec3 minX3_Offset(0, 0, minOffset.X3());
- walberla::pe::Vec3 maxX3_Offset(0, 0, maxOffset.X3());
-
- walberla::pe::Vec3 minCorner = simulationDomain.minCorner();
- walberla::pe::Vec3 maxCorner = simulationDomain.maxCorner();
-
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(1, 0, 0), minCorner + minX1_Offset, material);
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(-1, 0, 0), maxCorner + maxX1_Offset, material);
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 1, 0), minCorner + minX2_Offset, material);
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, -1, 0), maxCorner + maxX2_Offset, material);
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 0, 1), minCorner + minX3_Offset, material);
- createPlane(*globalBodyStorage, 0, walberla::pe::Vec3(0, 0, -1), maxCorner + maxX3_Offset, material);
-}
-
-std::shared_ptr<walberla::blockforest::BlockForest> PePhysicsEngineSolverAdapter::getForest() { return forest; }
-
-void PePhysicsEngineSolverAdapter::saveToFile(const std::string &path)
-{
- forest->saveToFile(path + "SerializeDeserialize.sbf");
- forest->saveBlockData("SerializeDeserialize.dump", *storageId.get());
-}
-
-void PePhysicsEngineSolverAdapter::loadFromFile(const std::string &path)
-{
- // forest = std::make_shared< walberla::blockforest::BlockForest >( walberla::uint_c(
- // walberla::MPIManager::instance()->rank() ), path+"SerializeDeserialize.sbf", true, false );
- std::string file = path + "SerializeDeserialize.sbf";
- forest = std::shared_ptr<walberla::blockforest::BlockForest>(new walberla::blockforest::BlockForest(
- walberla::uint_c(walberla::MPIManager::instance()->rank()), file.c_str(), true, false));
- storageId = std::make_shared<walberla::domain_decomposition::BlockDataID>(forest->loadBlockData(
- path + "SerializeDeserialize.dump", walberla::pe::createStorageDataHandling<BodyTypeTuple>(), "Storage"));
-
- this->initalPeIntegrator();
-
- auto ccdID =
- forest->addBlockData(walberla::pe::ccd::createHashGridsDataHandling(globalBodyStorage, *storageId), "CCD");
- auto fcdID = forest->addBlockData(
- walberla::pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, walberla::pe::fcd::AnalyticCollideFunctor>(),
- "FCD");
-
- cr = std::make_shared<walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers>(globalBodyStorage, forest,
- *storageId, ccdID, fcdID);
- cr->setMaxIterations(peParameter->maxPeIterations);
- cr->setRelaxationModel(
- walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers::ApproximateInelasticCoulombContactByDecoupling);
- cr->setRelaxationParameter(walberla::real_t(peParameter->relaxationParameter));
- cr->setGlobalLinearAcceleration(PeConverter::convert(peParameter->globalLinearAcceleration));
-
- this->executePeBodyTypeTuple();
- this->initialPeChannel();
-
- for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
- walberla::pe::ccd::ICCD *ccd = blockIt->getData<walberla::pe::ccd::ICCD>(ccdID);
- ccd->reloadBodies();
- }
-}
-
-std::shared_ptr<walberla::blockforest::BlockForest> PePhysicsEngineSolverAdapter::getBlockForest() { return forest; }
-
-std::shared_ptr<walberla::domain_decomposition::BlockDataID> PePhysicsEngineSolverAdapter::getStorageId()
-{
- return storageId;
-}
-
-std::shared_ptr<walberla::pe::BodyStorage> PePhysicsEngineSolverAdapter::getGlobalBodyStorage()
-{
- return globalBodyStorage;
-}
-
-void PePhysicsEngineSolverAdapter::createObstacle(const Vector3D ¢er, const Vector3D &lengths)
-{
- const walberla::pe::MaterialID material = peParameter->planes->getPeMaterial();
- bool global = true;
- bool communicating = false;
- bool infiniteMass = true;
-
- walberla::pe::createBox(*globalBodyStorage, *forest, *storageId, 0, PeConverter::convert(center),
- PeConverter::convert(lengths), material, global, communicating, infiniteMass);
-}
diff --git a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h b/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
deleted file mode 100644
index 5b2ef94b8a59b5073aebc7eb999db6e9eb860e3d..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/physicsEngineAdapter/pe/PePhysicsEngineSolverAdapter.h
+++ /dev/null
@@ -1,107 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef PE_PHYSICS_ENGINE_SOLVER_ADAPTER_H
-#define PE_PHYSICS_ENGINE_SOLVER_ADAPTER_H
-
-#include <memory>
-#include <shared_mutex>
-
-#include "UbTuple.h"
-#include <pe/basic.h>
-
-#include "PePhysicsEngineSolverAdapter.h"
-#include "PhysicsEngineSolverAdapter.h"
-
-class PePhysicsEngineMaterialAdapter;
-class PhysicsEngineGeometryAdapter;
-class PePhysicsEngineGeometryAdapter;
-class PeLoadBalancerAdapter;
-
-namespace walberla
-{
-namespace domain_decomposition
-{
-class BlockDataID;
-}
-namespace blockforest
-{
-class BlockForest;
-}
-namespace pe
-{
-class BodyStorage;
-class RigidBody;
-namespace cr
-{
-class HardContactSemiImplicitTimesteppingSolvers;
-}
-} // namespace pe
-} // namespace walberla
-
-struct PeParameter {
- PeParameter(double relaxationParameter, int maxPeIterations, Vector3D globalLinearAcceleration,
- std::shared_ptr<PePhysicsEngineMaterialAdapter> planes, std::array<double, 6> simulationDomain,
- UbTupleInt3 numberOfBlocks, UbTupleBool3 isPeriodic, Vector3D minOffset, Vector3D maxOffset)
- : relaxationParameter(relaxationParameter), maxPeIterations(maxPeIterations),
- globalLinearAcceleration(globalLinearAcceleration), simulationDomain(simulationDomain),
- numberOfBlocks(numberOfBlocks), isPeriodic(isPeriodic), planes(planes), minOffset(minOffset),
- maxOffset(maxOffset)
- {
- }
-
- double relaxationParameter;
- int maxPeIterations;
- Vector3D globalLinearAcceleration;
-
- std::array<double, 6> simulationDomain;
- UbTupleInt3 numberOfBlocks;
- UbTupleBool3 isPeriodic;
-
- std::shared_ptr<PePhysicsEngineMaterialAdapter> planes;
-
- Vector3D minOffset;
- Vector3D maxOffset;
-};
-
-class PePhysicsEngineSolverAdapter : public PhysicsEngineSolverAdapter
-{
-public:
- PePhysicsEngineSolverAdapter(std::shared_ptr<PeParameter> peParameter,
- std::shared_ptr<PeLoadBalancerAdapter> loadBalancer);
- virtual ~PePhysicsEngineSolverAdapter() {}
-
- std::shared_ptr<PhysicsEngineGeometryAdapter>
- createPhysicsEngineGeometryAdapter(int id, const Vector3D &position, double radius,
- std::shared_ptr<PhysicsEngineMaterialAdapter> material) const override;
- void runTimestep(double step) override;
- std::shared_ptr<walberla::blockforest::BlockForest> getForest();
- void saveToFile(const std::string &path);
- void loadFromFile(const std::string &path);
- std::shared_ptr<walberla::blockforest::BlockForest> getBlockForest();
- std::shared_ptr<walberla::domain_decomposition::BlockDataID> getStorageId();
- std::shared_ptr<walberla::pe::BodyStorage> getGlobalBodyStorage();
- void createObstacle(const Vector3D ¢er, const Vector3D &lengths);
-
-private:
- void initalizePeEnvironment();
- void initialPeBodyStorage();
- void initialPeBlockForest();
- void initalBlockData();
-
- void initalPeIntegrator();
- static void executePeBodyTypeTuple();
- void initialPeChannel() const;
-
-private:
- std::shared_ptr<PeParameter> peParameter;
- std::shared_ptr<PeLoadBalancerAdapter> loadBalancer;
-
- std::shared_ptr<walberla::pe::BodyStorage> globalBodyStorage;
- std::shared_ptr<walberla::blockforest::BlockForest> forest;
- std::shared_ptr<walberla::domain_decomposition::BlockDataID> storageId;
- std::shared_ptr<walberla::pe::cr::HardContactSemiImplicitTimesteppingSolvers> cr;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.cpp b/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.cpp
deleted file mode 100644
index c5486db0e23d6df29a04e724471bac74c01ceec8..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.cpp
+++ /dev/null
@@ -1,54 +0,0 @@
-#include "EquilibriumReconstructor.h"
-
-#include "BCArray3D.h"
-#include "BCProcessor.h"
-#include "D3Q27System.h"
-#include "DataSet3D.h"
-#include "ILBMKernel.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-
-void EquilibriumReconstructor::reconstructNode(const int &x1, const int &x2, const int &x3,
- const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- const double averageDensity = this->getLocalAverageDensity(x1, x2, x3, kernel);
- LBMReal feq[27];
- const Vector3D boundaryVelocity = physicsEngineGeometry->getVelocityAtPosition(worldCoordinates);
-
- if (kernel->getCompressible())
- D3Q27System::calcCompFeq(feq, averageDensity, boundaryVelocity[0], boundaryVelocity[1], boundaryVelocity[2]);
- else
- D3Q27System::calcIncompFeq(feq, averageDensity, boundaryVelocity[0], boundaryVelocity[1], boundaryVelocity[2]);
-
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
- // distributions->setDistribution(feq, x1, x2, x3);
- distributions->setDistributionInv(feq, x1, x2, x3);
-}
-
-double EquilibriumReconstructor::getLocalAverageDensity(const int &x1, const int &x2, const int &x3,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- int nAverage = 0;
- double averageDensity = 0.0;
-
- SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
-
- LBMReal f[D3Q27System::ENDF + 1];
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
-
- int neighborX1, neighborX2, neighborX3;
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++) {
- neighborX1 = x1 + D3Q27System::DX1[fDir];
- neighborX2 = x2 + D3Q27System::DX2[fDir];
- neighborX3 = x3 + D3Q27System::DX3[fDir];
-
- if (bcArray->isFluid(neighborX1, neighborX2, neighborX3)) {
- distributions->getDistribution(f, neighborX1, neighborX2, neighborX3);
- averageDensity += D3Q27System::getDensity(f);
- ++nAverage;
- }
- }
- return (nAverage > 0) ? averageDensity / nAverage : 0.0;
-}
diff --git a/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.h b/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.h
deleted file mode 100644
index 4a5c0071922645f66baeb6dcf8577f7003ca9078..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/EquilibriumReconstructor.h
+++ /dev/null
@@ -1,29 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef EQUILIBRIUM_RECONSTRUCTOR_H
-#define EQUILIBRIUM_RECONSTRUCTOR_H
-
-#include "UbTuple.h"
-
-#include "Reconstructor.h"
-
-class ILBMKernel;
-class PhysicsEngineGeometryAdapter;
-
-class EquilibriumReconstructor : public Reconstructor
-{
-public:
- virtual ~EquilibriumReconstructor() {}
-
- void reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const override;
-
-private:
- double getLocalAverageDensity(const int &x1, const int &x2, const int &x3,
- std::shared_ptr<ILBMKernel> kernel) const;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.cpp b/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.cpp
deleted file mode 100644
index 8bd915b6c130fd8a1b70cf9e54cbfdfc60b6efcc..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.cpp
+++ /dev/null
@@ -1,127 +0,0 @@
-#include "ExtrapolationReconstructor.h"
-
-#include "BCArray3D.h"
-#include "BCProcessor.h"
-#include "D3Q27System.h"
-#include "DataSet3D.h"
-#include "ILBMKernel.h"
-
-#include "DistributionArray3D.h"
-#include "PhysicsEngineGeometryAdapter.h"
-
-void ExtrapolationReconstructor::setAlternativeReconstructor(std::shared_ptr<Reconstructor> alternativeReconstructor)
-{
- this->alternativeReconstructor = alternativeReconstructor;
-}
-
-ExtrapolationReconstructor::ExtrapolationReconstructor(std::shared_ptr<Reconstructor> alternativeReconstructor)
- : alternativeReconstructor(alternativeReconstructor)
-{
-}
-
-void ExtrapolationReconstructor::reconstructNode(const int &x1, const int &x2, const int &x3,
- const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- const UbTupleInt3 extrapolationDirection = getSphereDirection(worldCoordinates, physicsEngineGeometry);
- const int numberOfCellsForExtrapolation = getNumberOfExtrapolationCells(x1, x2, x3, extrapolationDirection, kernel);
-
- // if (numberOfCellsForExtrapolation < 2)
- alternativeReconstructor->reconstructNode(x1, x2, x3, worldCoordinates, physicsEngineGeometry, kernel);
- // else
- //{
- // //UBLOG(logINFO, "point (x,y,z) " << val<1>(worldCoordinates) << ", " << val<2>(worldCoordinates) << ", " <<
- // val<3>(worldCoordinates));
- // //UBLOG(logINFO, "extradir (x,y,z) " << val<1>(extrapolationDirection) << ", " <<
- // val<2>(extrapolationDirection) << ", " << val<3>(extrapolationDirection));
- // //UBLOG(logINFO, "numberOfCellsForExtrapolation: " << numberOfCellsForExtrapolation );
-
- // this->extrapolatePdFs(x1, x2, x3, extrapolationDirection, numberOfCellsForExtrapolation, kernel);
- //}
-}
-
-UbTupleInt3 ExtrapolationReconstructor::getSphereDirection(
- const Vector3D &worldCoordinates, std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry) const
-{
- const Vector3D spherePosition = physicsEngineGeometry->getPosition();
- const Vector3D bodyNormal = worldCoordinates - spherePosition;
- return this->getCorrespondingLatticeDirection(bodyNormal);
-}
-
-UbTupleInt3 ExtrapolationReconstructor::getCorrespondingLatticeDirection(const Vector3D &direction) const
-{
- int correspondingDirection = 0;
- double innerProduct = 0.0;
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++) {
- // compute inner product <dir,c_i>
- const double temporaryInnerProduct = direction[0] * D3Q27System::cNorm[0][fDir] +
- direction[1] * D3Q27System::cNorm[1][fDir] +
- direction[2] * D3Q27System::cNorm[2][fDir];
- if (temporaryInnerProduct > innerProduct) {
- innerProduct = temporaryInnerProduct;
- correspondingDirection = fDir;
- }
- }
-
- return UbTupleInt3(D3Q27System::DX1[correspondingDirection], D3Q27System::DX2[correspondingDirection],
- D3Q27System::DX3[correspondingDirection]);
-}
-
-int ExtrapolationReconstructor::getNumberOfExtrapolationCells(const int x1, const int x2, const int x3,
- const UbTupleInt3 &extrapolationDirection,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- if (extrapolationDirection == UbTupleInt3(0, 0, 0))
- return 0;
-
- const int desiredCellsInExtrapolationDirection = 3;
-
- for (int numCells = 1; numCells <= desiredCellsInExtrapolationDirection; ++numCells) {
- UbTupleInt3 neighbor(x1 + numCells * val<1>(extrapolationDirection),
- x2 + numCells * val<2>(extrapolationDirection),
- x3 + numCells * val<3>(extrapolationDirection));
-
- if (!kernel->isInsideOfDomain(val<1>(neighbor), val<2>(neighbor), val<3>(neighbor)))
- return numCells - 1;
-
- if (!kernel->getBCProcessor()->getBCArray()->isFluid(val<1>(neighbor), val<2>(neighbor), val<3>(neighbor)))
- return numCells - 1;
- }
- return desiredCellsInExtrapolationDirection;
-}
-
-void ExtrapolationReconstructor::extrapolatePdFs(const int x1, const int x2, const int x3,
- const UbTupleInt3 &extrapolationDirection,
- int numberOfCellsForExtrapolation,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
-
- const int nx1 = val<1>(extrapolationDirection);
- const int nx2 = val<2>(extrapolationDirection);
- const int nx3 = val<3>(extrapolationDirection);
-
- LBMReal pdf[D3Q27System::ENDF + 1];
- LBMReal pdfNeighbor1[D3Q27System::ENDF + 1];
- LBMReal pdfNeighbor2[D3Q27System::ENDF + 1];
-
- distributions->getDistribution(pdf, x1, x2, x3);
- distributions->getDistribution(pdfNeighbor1, x1 + nx1, x2 + nx2, x3 + nx3);
- distributions->getDistribution(pdfNeighbor2, x1 + 2 * nx1, x2 + 2 * nx2, x3 + 2 * nx3);
-
- if (numberOfCellsForExtrapolation == 3) // quadratic normal extrapolation
- {
- LBMReal pdfNeighbor3[D3Q27System::ENDF + 1];
- distributions->getDistribution(pdfNeighbor3, x1 + 3 * nx1, x2 + 3 * nx2, x3 + 3 * nx3);
-
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++)
- pdf[fDir] = 3 * pdfNeighbor1[fDir] - 3 * pdfNeighbor2[fDir] + pdfNeighbor3[fDir];
- } else // numberOfCellsForExtrapolation == 2 // linear normal extrapolation
- {
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++)
- pdf[fDir] = 2 * pdfNeighbor1[fDir] - pdfNeighbor2[fDir];
- }
-
- distributions->setDistribution(pdf, x1, x2, x3);
-}
diff --git a/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.h b/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.h
deleted file mode 100644
index 2844fdf2490a43e6db31ea9bf32dbcdea67d6f34..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/ExtrapolationReconstructor.h
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef EXTRAPOLATION_RECONSTRUCTOR_H
-#define EXTRAPOLATION_RECONSTRUCTOR_H
-
-#include <memory>
-
-#include "UbTuple.h"
-
-#include "Reconstructor.h"
-
-class ILBMKernel;
-class PhysicsEngineGeometryAdapter;
-
-class ExtrapolationReconstructor : public Reconstructor
-{
-public:
- ExtrapolationReconstructor(std::shared_ptr<Reconstructor> alternativeReconstructor);
- virtual ~ExtrapolationReconstructor() {}
-
- void reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const override;
-
- void setAlternativeReconstructor(std::shared_ptr<Reconstructor> alternativeReconstructor);
-
-private:
- int getNumberOfExtrapolationCells(const int x1, const int x2, const int x3, const UbTupleInt3 &ubTuple,
- std::shared_ptr<ILBMKernel> kernel) const;
- UbTupleInt3 getSphereDirection(const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry) const;
- UbTupleInt3 getCorrespondingLatticeDirection(const Vector3D &direction) const;
- void extrapolatePdFs(const int x1, const int x2, const int x3, const UbTupleInt3 &ubTuple,
- int numberOfCellsForExtrapolation, std::shared_ptr<ILBMKernel> kernel) const;
-
- std::shared_ptr<Reconstructor> alternativeReconstructor;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/reconstructor/LBMReconstructor.cpp b/src/cpu/DemCoupling/reconstructor/LBMReconstructor.cpp
deleted file mode 100644
index c6dfdc6dfbe064918fc50faf3222ed14d2b3a8d3..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/LBMReconstructor.cpp
+++ /dev/null
@@ -1,133 +0,0 @@
-#include "LBMReconstructor.h"
-
-#include "BCArray3D.h"
-#include "BCProcessor.h"
-#include "D3Q27System.h"
-#include "DataSet3D.h"
-#include "ILBMKernel.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-
-using namespace D3Q27System;
-
-LBMReconstructor::LBMReconstructor(bool compressible)
-{
- if (compressible) {
- calcMacrosFct = &D3Q27System::calcCompMacroscopicValues;
- } else {
- calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues;
- }
-}
-
-void LBMReconstructor::reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- LBMReal pdf[D3Q27System::ENDF + 1];
-
- LBMReal rho, vx1, vx2, vx3;
- calcMacrosFct(pdf, rho, vx1, vx2, vx3);
-
- LBMReal rho_dif = 1;
-
- while (rho_dif > 1e-5) {
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++) {
-
- UbTupleInt3 neighbor(x1 + D3Q27System::DX1[fDir], x2 + D3Q27System::DX2[fDir], x3 + D3Q27System::DX3[fDir]);
-
- if (!kernel->getBCProcessor()->getBCArray()->isFluid(val<1>(neighbor), val<2>(neighbor),
- val<3>(neighbor))) {
- LBMReal pdfNeighbor[D3Q27System::ENDF + 1];
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
- const int invDir = D3Q27System::INVDIR[fDir];
- distributions->getDistributionForDirection(pdfNeighbor[invDir], val<1>(neighbor), val<2>(neighbor),
- val<3>(neighbor));
- distributions->setDistributionInvForDirection(pdf[invDir], x1, x2, x3, invDir);
- }
- }
- }
-
- LBMReal collFactor = kernel->getCollisionFactor();
- collide(pdf, collFactor);
-}
-
-void LBMReconstructor::collide(LBMReal *f, LBMReal collFactor)
-{
-
- LBMReal drho, vx1, vx2, vx3;
- LBMReal feq[D3Q27System::ENDF + 1];
-
- drho = ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
- (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
- ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
- ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[ZERO];
-
- vx1 = ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[BSE] - f[TNW]) + (f[BNE] - f[TSW]))) +
- (((f[BE] - f[TW]) + (f[TE] - f[BW])) + ((f[SE] - f[NW]) + (f[NE] - f[SW]))) + (f[E] - f[W]));
-
- vx2 = ((((f[TNE] - f[BSW]) + (f[BNW] - f[TSE])) + ((f[TNW] - f[BSE]) + (f[BNE] - f[TSW]))) +
- (((f[BN] - f[TS]) + (f[TN] - f[BS])) + ((f[NW] - f[SE]) + (f[NE] - f[SW]))) + (f[N] - f[S]));
-
- vx3 = ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[TNW] - f[BSE]) + (f[TSW] - f[BNE]))) +
- (((f[TS] - f[BN]) + (f[TN] - f[BS])) + ((f[TW] - f[BE]) + (f[TE] - f[BW]))) + (f[T] - f[B]));
-
- LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
-
- feq[ZERO] = c8o27 * (drho - cu_sq);
- feq[E] = c2o27 * (drho + 3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq);
- feq[W] = c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
- feq[N] = c2o27 * (drho + 3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq);
- feq[S] = c2o27 * (drho + 3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq);
- feq[T] = c2o27 * (drho + 3.0 * (vx3) + c9o2 * (vx3) * (vx3)-cu_sq);
- feq[B] = c2o27 * (drho + 3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq);
- feq[NE] = c1o54 * (drho + 3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq);
- feq[SW] = c1o54 * (drho + 3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq);
- feq[SE] = c1o54 * (drho + 3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq);
- feq[NW] = c1o54 * (drho + 3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq);
- feq[TE] = c1o54 * (drho + 3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq);
- feq[BW] = c1o54 * (drho + 3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq);
- feq[BE] = c1o54 * (drho + 3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq);
- feq[TW] = c1o54 * (drho + 3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq);
- feq[TN] = c1o54 * (drho + 3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq);
- feq[BS] = c1o54 * (drho + 3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq);
- feq[BN] = c1o54 * (drho + 3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq);
- feq[TS] = c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq);
- feq[TNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq);
- feq[BSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
- feq[BNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq);
- feq[TSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
- feq[TSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq);
- feq[BNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
- feq[BSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
- feq[TNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
-
- // Relaxation
- f[ZERO] += (feq[ZERO] - f[ZERO]) * collFactor;
- f[E] += (feq[E] - f[E]) * collFactor;
- f[W] += (feq[W] - f[W]) * collFactor;
- f[N] += (feq[N] - f[N]) * collFactor;
- f[S] += (feq[S] - f[S]) * collFactor;
- f[T] += (feq[T] - f[T]) * collFactor;
- f[B] += (feq[B] - f[B]) * collFactor;
- f[NE] += (feq[NE] - f[NE]) * collFactor;
- f[SW] += (feq[SW] - f[SW]) * collFactor;
- f[SE] += (feq[SE] - f[SE]) * collFactor;
- f[NW] += (feq[NW] - f[NW]) * collFactor;
- f[TE] += (feq[TE] - f[TE]) * collFactor;
- f[BW] += (feq[BW] - f[BW]) * collFactor;
- f[BE] += (feq[BE] - f[BE]) * collFactor;
- f[TW] += (feq[TW] - f[TW]) * collFactor;
- f[TN] += (feq[TN] - f[TN]) * collFactor;
- f[BS] += (feq[BS] - f[BS]) * collFactor;
- f[BN] += (feq[BN] - f[BN]) * collFactor;
- f[TS] += (feq[TS] - f[TS]) * collFactor;
-
- f[TNE] += (feq[TNE] - f[TNE]) * collFactor;
- f[BSW] += (feq[BSW] - f[BSW]) * collFactor;
- f[BNE] += (feq[BNE] - f[BNE]) * collFactor;
- f[TSW] += (feq[TSW] - f[TSW]) * collFactor;
- f[TSE] += (feq[TSE] - f[TSE]) * collFactor;
- f[BNW] += (feq[BNW] - f[BNW]) * collFactor;
- f[BSE] += (feq[BSE] - f[BSE]) * collFactor;
- f[TNW] += (feq[TNW] - f[TNW]) * collFactor;
-}
\ No newline at end of file
diff --git a/src/cpu/DemCoupling/reconstructor/LBMReconstructor.h b/src/cpu/DemCoupling/reconstructor/LBMReconstructor.h
deleted file mode 100644
index 173f008a30bc9f1edf4b79eecce87fc0941d9f62..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/LBMReconstructor.h
+++ /dev/null
@@ -1,35 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef LBM_RECONSTRUCTOR_H
-#define LBM_RECONSTRUCTOR_H
-
-#include "UbTuple.h"
-
-#include "Reconstructor.h"
-
-#include "LBMSystem.h"
-
-class ILBMKernel;
-class PhysicsEngineGeometryAdapter;
-
-class LBMReconstructor : public Reconstructor
-{
-public:
- LBMReconstructor(bool compressible);
- virtual ~LBMReconstructor() {}
-
- void reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const override;
-
-private:
- static void collide(LBMReal *f, LBMReal collFactor);
-
- typedef void (*CalcMacrosFct)(const LBMReal *const & /*f[27]*/, LBMReal & /*rho*/, LBMReal & /*vx1*/,
- LBMReal & /*vx2*/, LBMReal & /*vx3*/);
- CalcMacrosFct calcMacrosFct;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/reconstructor/Reconstructor.h b/src/cpu/DemCoupling/reconstructor/Reconstructor.h
deleted file mode 100644
index 15355d515dc332521583955f225c3e8758bb5fb7..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/Reconstructor.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef RECONSTRCUTOR_H
-#define RECONSTRCUTOR_H
-
-#include <PointerDefinitions.h>
-
-#include "Vector3D.h"
-
-class ILBMKernel;
-class PhysicsEngineGeometryAdapter;
-
-class Reconstructor
-{
-public:
- virtual ~Reconstructor() {}
-
- virtual void reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- SPtr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const = 0;
-};
-
-#endif
diff --git a/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.cpp b/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.cpp
deleted file mode 100644
index c48586ca4ed170d7129990590a7e7ec32154c5f7..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.cpp
+++ /dev/null
@@ -1,96 +0,0 @@
-#include "VelocityBcReconstructor.h"
-
-#include <exception>
-
-#include "BCArray3D.h"
-#include "BCProcessor.h"
-#include "D3Q27System.h"
-#include "DataSet3D.h"
-#include "EsoTwist3D.h"
-#include "ILBMKernel.h"
-
-#include "PhysicsEngineGeometryAdapter.h"
-
-void VelocityBcReconstructor::reconstructNode(const int &x1, const int &x2, const int &x3,
- const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const
-{
- if (kernel->getCompressible())
- throw std::runtime_error("not implemented yet!");
-
- const Vector3D boundaryVelocity = physicsEngineGeometry->getVelocityAtPosition(worldCoordinates);
- // TODO: move to D3Q27 system
- LBMReal wijk[D3Q27System::ENDF + 1];
- D3Q27System::calcIncompFeq(wijk, 1, 0, 0, 0);
-
- SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
-
- SPtr<BoundaryConditions> bc = SPtr<BoundaryConditions>(new BoundaryConditions());
- bc->setBoundaryVelocityX1((float)boundaryVelocity[0]);
- bc->setBoundaryVelocityX2((float)boundaryVelocity[1]);
- bc->setBoundaryVelocityX3((float)boundaryVelocity[2]);
-
- LBMReal feqNullRho[D3Q27System::ENDF + 1];
- D3Q27System::calcIncompFeq(feqNullRho, 0, boundaryVelocity[0], boundaryVelocity[1], boundaryVelocity[2]);
-
- LBMReal fpre[D3Q27System::ENDF + 1];
- LBMReal fpost[D3Q27System::ENDF + 1];
- SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
-
- distributions->swap();
- distributions->getDistributionInv(fpost, x1, x2, x3);
- distributions->swap();
- distributions->getDistribution(fpre, x1, x2, x3);
-
- int neighborX1, neighborX2, neighborX3;
- int neighborX1Inv, neighborX2Inv, neighborX3Inv;
-
- double sumRho = 0, sumWijk = 0;
- double collFactor = kernel->getCollisionFactor();
-
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++) {
- neighborX1 = x1 + D3Q27System::DX1[fDir];
- neighborX2 = x2 + D3Q27System::DX2[fDir];
- neighborX3 = x3 + D3Q27System::DX3[fDir];
-
- if (bcArray->isFluid(neighborX1, neighborX2, neighborX3)) {
- int invDir = D3Q27System::INVDIR[fDir];
-
- neighborX1Inv = x1 + D3Q27System::DX1[invDir];
- neighborX2Inv = x2 + D3Q27System::DX2[invDir];
- neighborX3Inv = x3 + D3Q27System::DX3[invDir];
- if (!bcArray->isFluid(neighborX1Inv, neighborX2Inv, neighborX3Inv)) {
-
- double velocity = bc->getBoundaryVelocity(invDir);
-
- fpre[fDir] = fpre[invDir] - velocity;
- double Omega = fpost[fDir] - fpre[fDir];
-
- sumRho += Omega / collFactor + fpre[fDir] - feqNullRho[fDir];
- sumWijk += wijk[fDir];
- }
- }
- }
-
- double rho = 0.0;
- if (sumWijk > 0.0)
- rho = sumRho / sumWijk;
-
- for (int fDir = D3Q27System::FSTARTDIR; fDir <= D3Q27System::FENDDIR; fDir++) {
- neighborX1 = x1 + D3Q27System::DX1[fDir];
- neighborX2 = x2 + D3Q27System::DX2[fDir];
- neighborX3 = x3 + D3Q27System::DX3[fDir];
-
- if (!bcArray->isFluid(neighborX1, neighborX2, neighborX3)) {
- int invDir = D3Q27System::INVDIR[fDir];
- neighborX1Inv = x1 + D3Q27System::DX1[invDir];
- neighborX2Inv = x2 + D3Q27System::DX2[invDir];
- neighborX3Inv = x3 + D3Q27System::DX3[invDir];
- if (!bcArray->isFluid(neighborX1Inv, neighborX2Inv, neighborX3Inv)) {
- fpre[fDir] = D3Q27System::getIncompFeqForDirection(
- fDir, rho, bc->getBoundaryVelocityX1(), bc->getBoundaryVelocityX2(), bc->getBoundaryVelocityX3());
- }
- }
- }
-}
diff --git a/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.h b/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.h
deleted file mode 100644
index 9f5b3f0b67be91edbcfcdadbc8f5c637a87827dc..0000000000000000000000000000000000000000
--- a/src/cpu/DemCoupling/reconstructor/VelocityBcReconstructor.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * Author: S. Peters
- * mail: peters@irmb.tu-bs.de
- */
-#ifndef VELOCITY_BC_RECONSTRUCTOR_H
-#define VELOCITY_BC_RECONSTRUCTOR_H
-
-#include "UbTuple.h"
-
-#include "Reconstructor.h"
-
-class ILBMKernel;
-class PhysicsEngineGeometryAdapter;
-
-class VelocityBcReconstructor : public Reconstructor
-{
-public:
- virtual ~VelocityBcReconstructor() {}
-
- void reconstructNode(const int &x1, const int &x2, const int &x3, const Vector3D &worldCoordinates,
- std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry,
- std::shared_ptr<ILBMKernel> kernel) const override;
-};
-
-#endif
diff --git a/src/cpu/VirtualFluidsCore/CMakeLists.txt b/src/cpu/VirtualFluidsCore/CMakeLists.txt
index 0485e00ec407dcfe0e2ac1c2a58da044fb33f42d..5300e898bd17e45b2a0a5c5e2b0d083c975ad1fb 100644
--- a/src/cpu/VirtualFluidsCore/CMakeLists.txt
+++ b/src/cpu/VirtualFluidsCore/CMakeLists.txt
@@ -16,10 +16,6 @@ IF(${USE_CATALYST})
list(APPEND VF_LIBRARIES optimized vtkParallelMPI debug vtkParallelMPI )
ENDIF()
-IF(${USE_DEM_COUPLING})
- INCLUDE(${CMAKE_CURRENT_SOURCE_DIR}/../DemCoupling/DemCoupling.cmake)
-ENDIF()
-
if(BUILD_USE_OPENMP)
list(APPEND VF_LIBRARIES OpenMP::OpenMP_CXX)
endif()