From 6186a067350bdb053cde85273c5e932297123fad Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Thu, 7 Jun 2018 18:26:57 +0200
Subject: [PATCH] - moves BoundaryConditionVisitor update from
MovableObjectInteractor to DemCoProcessor - fixes BC nodes update in
MovableObjectInteractor
---
source/Applications/PoiseuilleFlow/pf1.cpp | 2 +-
.../Applications/Thermoplast/thermoplast.cpp | 23 +-
.../CreateGeoObjectsCoProcessor.cpp | 9 +-
.../DemCoupling/CreateGeoObjectsCoProcessor.h | 4 +-
source/DemCoupling/DemCoProcessor.cpp | 8 +-
source/DemCoupling/DemCoProcessor.h | 5 +
.../DemCoupling/MovableObjectInteractor.cpp | 28 +-
source/DemCoupling/MovableObjectInteractor.h | 5 +-
.../numerics/geometry3d/GbCuboid3D.cpp | 292 +++++++++---------
.../WriteBoundaryConditionsCoProcessor.cpp | 4 +-
.../Grid/BasicCalculator.cpp | 3 -
11 files changed, 191 insertions(+), 192 deletions(-)
diff --git a/source/Applications/PoiseuilleFlow/pf1.cpp b/source/Applications/PoiseuilleFlow/pf1.cpp
index 85582c56f..1df2bffaa 100644
--- a/source/Applications/PoiseuilleFlow/pf1.cpp
+++ b/source/Applications/PoiseuilleFlow/pf1.cpp
@@ -139,7 +139,7 @@ void pf1()
//write data for visualization of boundary conditions
{
SPtr<UbScheduler> geoSch(new UbScheduler(1));
- WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm);
+ WriteBoundaryConditionsCoProcessor ppgeo(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm);
ppgeo.process(0);
}
diff --git a/source/Applications/Thermoplast/thermoplast.cpp b/source/Applications/Thermoplast/thermoplast.cpp
index 18e20f509..0af46a907 100644
--- a/source/Applications/Thermoplast/thermoplast.cpp
+++ b/source/Applications/Thermoplast/thermoplast.cpp
@@ -76,7 +76,7 @@ void pf1()
double deltax = 1;
double rhoLB = 0.0;
double nuLB = 0.005;
- double uLB = 0.01;
+ double uLB = 0.1;
//geometry definition
@@ -94,10 +94,6 @@ void pf1()
SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
- //cylinder
- //SPtr<GbObject3D> cylinder(new GbCylinder3D(g_minX1 - 2.0*deltax, 0.0, 0.0, g_maxX1 + 2.0*deltax, 0.0, 0.0, g_maxX2));
- //GbSystem3D::writeGeoObject(cylinder.get(), pathname + "/geo/cylinder", WbWriterVtkXmlBinary::getInstance());
-
//box
SPtr<GbObject3D> box(new GbCuboid3D(g_minX1-blockLength, g_minX2, g_minX3, g_maxX1+blockLength, g_maxX2, g_maxX3));
GbSystem3D::writeGeoObject(box.get(), pathname + "/geo/box", WbWriterVtkXmlBinary::getInstance());
@@ -138,7 +134,7 @@ void pf1()
fct.SetExpr("U");
fct.DefineConst("U", uLB);
SPtr<BCAdapter> inflowAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+ inflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
SPtr<BCAdapter> outflowAdapter(new DensityBCAdapter(rhoLB));
outflowAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
@@ -167,7 +163,7 @@ void pf1()
//////////////////////////////////////////////////////////////////////////////////
//set boundary conditions for blocks and create process decomposition for MPI
- SPtr<D3Q27Interactor> cylinderInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
+ SPtr<D3Q27Interactor> boxInt(new D3Q27Interactor(box, grid, noSlipBCAdapter, Interactor3D::INVERSESOLID));
//sphere bc object
SPtr<MovableObjectInteractor> sphereInt1;
@@ -175,8 +171,6 @@ void pf1()
const std::shared_ptr<Reconstructor> extrapolationReconstructor = std::make_shared<ExtrapolationReconstructor>(velocityReconstructor);
sphereInt1 = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(sphere1, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
- sphereInt1->setBlockVisitor(bcVisitor);
-
//inflow
SPtr<D3Q27Interactor> inflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, inflowAdapter, Interactor3D::SOLID));
@@ -186,8 +180,8 @@ void pf1()
SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::B));
InteractorsHelper intHelper(grid, metisVisitor);
- intHelper.addInteractor(cylinderInt);
- //intHelper.addInteractor(sphereInt1);
+ intHelper.addInteractor(boxInt);
+ intHelper.addInteractor(sphereInt1);
intHelper.addInteractor(inflowInt);
intHelper.addInteractor(outflowInt);
@@ -289,13 +283,14 @@ void pf1()
const int timestep = 2;
const SPtr<UbScheduler> peScheduler(new UbScheduler(timestep));
SPtr<DemCoProcessor> demCoProcessor = makePeCoProcessor(grid, comm, peScheduler, lbmUnitConverter);
- //demCoProcessor->addInteractor(sphereInt1, sphereMaterial, Vector3D(0.0, 0.0, 0.0));
+ demCoProcessor->addInteractor(sphereInt1, sphereMaterial, Vector3D(0.0, 0.0, 0.0));
demCoProcessor->distributeIDs();
+ demCoProcessor->setBlockVisitor(bcVisitor);
double g = 9.81 * lbmUnitConverter->getFactorAccWToLb();
double rhoFluid = lbmUnitConverter->getFactorDensityWToLb() * 830.0; // 1 // kg/m^3
//////////////////////////////////////////////////////////////////////////
SPtr<UbScheduler> sphereScheduler(new UbScheduler(2.0*radius/uLB,1));
- SPtr<CreateGeoObjectsCoProcessor> createSphereCoProcessor(new CreateGeoObjectsCoProcessor(grid,sphereScheduler,demCoProcessor,sphere1,sphereMaterial,bcVisitor));
+ SPtr<CreateGeoObjectsCoProcessor> createSphereCoProcessor(new CreateGeoObjectsCoProcessor(grid,sphereScheduler,demCoProcessor,sphere1,sphereMaterial));
//start simulation
@@ -304,7 +299,7 @@ void pf1()
SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
//SPtr<Calculator> calculator(new DemCalculator(grid, stepGhostLayer, endTime));
calculator->addCoProcessor(npr);
- calculator->addCoProcessor(createSphereCoProcessor);
+ //calculator->addCoProcessor(createSphereCoProcessor);
calculator->addCoProcessor(demCoProcessor);
calculator->addCoProcessor(writeBCCoProcessor);
calculator->addCoProcessor(writeMQCoProcessor);
diff --git a/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp b/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
index f631b196f..f43abd271 100644
--- a/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
+++ b/source/DemCoupling/CreateGeoObjectsCoProcessor.cpp
@@ -4,7 +4,6 @@
#include "GbSphere3D.h"
#include "VelocityBCAdapter.h"
#include "VelocityWithDensityBCAlgorithm.h"
-#include "BoundaryConditionsBlockVisitor.h"
#include "MovableObjectInteractor.h"
#include "VelocityBcReconstructor.h"
#include "ExtrapolationReconstructor.h"
@@ -15,12 +14,11 @@
#include "Grid3D.h"
CreateGeoObjectsCoProcessor::CreateGeoObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<DemCoProcessor> demCoProcessor,
- SPtr<GbObject3D> geoObjectPrototype, SPtr<PhysicsEngineMaterialAdapter> geoObjectMaterial, SPtr<BoundaryConditionsBlockVisitor> bcVisitor) :
+ SPtr<GbObject3D> geoObjectPrototype, SPtr<PhysicsEngineMaterialAdapter> geoObjectMaterial) :
CoProcessor(grid, s),
demCoProcessor(demCoProcessor),
geoObjectPrototype(geoObjectPrototype),
- geoObjectMaterial(geoObjectMaterial),
- bcVisitor(bcVisitor)
+ geoObjectMaterial(geoObjectMaterial)
{
}
@@ -38,16 +36,13 @@ void CreateGeoObjectsCoProcessor::process(double step)
SPtr<MovableObjectInteractor> geoObjectInt;
const std::shared_ptr<Reconstructor> velocityReconstructor(new VelocityBcReconstructor());
const std::shared_ptr<Reconstructor> extrapolationReconstructor(new ExtrapolationReconstructor(velocityReconstructor));
- //SPtr<GbObject3D> geoObject(dynamicPointerCast<GbSphere3D>(geoObjectPrototype)->clone());
SPtr<GbObject3D> geoObject((GbObject3D*)(geoObjectPrototype->clone()));
geoObjectInt = SPtr<MovableObjectInteractor>(new MovableObjectInteractor(geoObject, grid, velocityBcParticleAdapter, Interactor3D::SOLID, extrapolationReconstructor, State::UNPIN));
- geoObjectInt->setBlockVisitor(bcVisitor);
//SetSolidOrBoundaryBlockVisitor setSolidVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::SOLID);
//grid->accept(setSolidVisitor);
SetSolidOrBoundaryBlockVisitor setBcVisitor(geoObjectInt, SetSolidOrBoundaryBlockVisitor::BC);
grid->accept(setBcVisitor);
geoObjectInt->initInteractor();
- grid->accept(*bcVisitor.get());
demCoProcessor->addInteractor(geoObjectInt, geoObjectMaterial, Vector3D(0.0, 0.0, 0.0));
demCoProcessor->distributeIDs();
diff --git a/source/DemCoupling/CreateGeoObjectsCoProcessor.h b/source/DemCoupling/CreateGeoObjectsCoProcessor.h
index e3c717d5d..ff7a9eef6 100644
--- a/source/DemCoupling/CreateGeoObjectsCoProcessor.h
+++ b/source/DemCoupling/CreateGeoObjectsCoProcessor.h
@@ -7,7 +7,6 @@ class Grid3D;
class UbScheduler;
class DemCoProcessor;
class GbObject3D;
-class BoundaryConditionsBlockVisitor;
class PhysicsEngineMaterialAdapter;
@@ -16,13 +15,12 @@ class CreateGeoObjectsCoProcessor : public CoProcessor
public:
CreateGeoObjectsCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
SPtr<DemCoProcessor> demCoProcessor, SPtr<GbObject3D> spherePrototype,
- SPtr<PhysicsEngineMaterialAdapter> sphereMaterial, SPtr<BoundaryConditionsBlockVisitor> bcVisitor);
+ SPtr<PhysicsEngineMaterialAdapter> sphereMaterial);
void process(double step) override;
protected:
private:
SPtr<DemCoProcessor> demCoProcessor;
SPtr<GbObject3D> geoObjectPrototype;
SPtr<PhysicsEngineMaterialAdapter> geoObjectMaterial;
- SPtr<BoundaryConditionsBlockVisitor> bcVisitor;
};
#endif // CreateSphereCoProcessor_h__
diff --git a/source/DemCoupling/DemCoProcessor.cpp b/source/DemCoupling/DemCoProcessor.cpp
index f162a2b3f..00c8b4aff 100644
--- a/source/DemCoupling/DemCoProcessor.cpp
+++ b/source/DemCoupling/DemCoProcessor.cpp
@@ -20,7 +20,7 @@
#include "Block3D.h"
#include "BCArray3D.h"
#include "MPICommunicator.h"
-
+#include "BoundaryConditionsBlockVisitor.h"
DemCoProcessor::DemCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Communicator> comm, std::shared_ptr<ForceCalculator> forceCalculator, std::shared_ptr<PhysicsEngineSolverAdapter> physicsEngineSolver, double intermediatePeSteps) :
CoProcessor(grid, s), comm(comm), forceCalculator(forceCalculator), physicsEngineSolver(physicsEngineSolver), intermediateDemSteps(intermediatePeSteps)
@@ -100,6 +100,8 @@ void DemCoProcessor::process(double actualTimeStep)
this->moveVfGeoObject();
+ grid->accept(*boundaryConditionsBlockVisitor.get());
+
//UBLOG(logINFO, "DemCoProcessor::update - END" << step);
}
}
@@ -252,3 +254,7 @@ void DemCoProcessor::distributeIDs()
}
}
+void DemCoProcessor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
+{
+ this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
+}
diff --git a/source/DemCoupling/DemCoProcessor.h b/source/DemCoupling/DemCoProcessor.h
index bb4606b03..fc5436ad6 100644
--- a/source/DemCoupling/DemCoProcessor.h
+++ b/source/DemCoupling/DemCoProcessor.h
@@ -23,6 +23,7 @@ class ForceCalculator;
class Communicator;
class MovableObjectInteractor;
class Communicator;
+class BoundaryConditionsBlockVisitor;
class DemCoProcessor : public CoProcessor
@@ -38,6 +39,8 @@ public:
std::shared_ptr<PhysicsEngineSolverAdapter> getPhysicsEngineSolver();
void distributeIDs();
+
+ void setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> blockVisitor);
private:
std::shared_ptr<PhysicsEngineGeometryAdapter> createPhysicsEngineGeometryAdapter(std::shared_ptr<MovableObjectInteractor> interactor, std::shared_ptr<PhysicsEngineMaterialAdapter> physicsEngineMaterial) const;
@@ -61,6 +64,8 @@ private:
double intermediateDemSteps;
+ SPtr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
+
};
diff --git a/source/DemCoupling/MovableObjectInteractor.cpp b/source/DemCoupling/MovableObjectInteractor.cpp
index 2a76405e9..29e340927 100644
--- a/source/DemCoupling/MovableObjectInteractor.cpp
+++ b/source/DemCoupling/MovableObjectInteractor.cpp
@@ -35,11 +35,6 @@ void MovableObjectInteractor::setPhysicsEngineGeometry(std::shared_ptr<PhysicsEn
physicsEngineGeometry->changeState(this->state);
}
-void MovableObjectInteractor::setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor)
-{
- this->boundaryConditionsBlockVisitor = boundaryConditionsBlockVisitor;
-}
-
void MovableObjectInteractor::moveGbObjectTo(const Vector3D& position)
{
//UBLOG(logINFO, "new position: (x,y,z) " << val<1>(position) << ", " << val<2>(position) << ", " << val<3>(position));
@@ -51,7 +46,9 @@ void MovableObjectInteractor::moveGbObjectTo(const Vector3D& position)
void MovableObjectInteractor::rearrangeGrid()
{
this->reconstructDistributionOnSolidNodes();
+
this->setSolidNodesToFluid();
+ this->setBcNodesToFluid();
this->removeSolidBlocks();
this->removeBcBlocks();
@@ -60,7 +57,6 @@ void MovableObjectInteractor::rearrangeGrid()
this->initInteractor();
- this->setBcs();
this->updateVelocityBc();
}
@@ -103,15 +99,25 @@ void MovableObjectInteractor::setSolidNodesToFluid()
}
}
-void MovableObjectInteractor::setBcBlocks()
+void MovableObjectInteractor::setBcNodesToFluid()
{
- SetSolidOrBoundaryBlockVisitor v(shared_from_this(), SetSolidOrBoundaryBlockVisitor::BC);
- this->grid.lock()->accept(v);
+ 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::setBcs() const
+void MovableObjectInteractor::setBcBlocks()
{
- this->grid.lock()->accept(*boundaryConditionsBlockVisitor.get());
+ SetSolidOrBoundaryBlockVisitor v(shared_from_this(), SetSolidOrBoundaryBlockVisitor::BC);
+ this->grid.lock()->accept(v);
}
void MovableObjectInteractor::updateVelocityBc()
diff --git a/source/DemCoupling/MovableObjectInteractor.h b/source/DemCoupling/MovableObjectInteractor.h
index de82d2604..2b43d97dc 100644
--- a/source/DemCoupling/MovableObjectInteractor.h
+++ b/source/DemCoupling/MovableObjectInteractor.h
@@ -19,7 +19,6 @@ class Block3D;
class BCArray3D;
class BCAdapter;
class GbObject3D;
-class BoundaryConditionsBlockVisitor;
class PhysicsEngineGeometryAdapter;
class Reconstructor;
@@ -31,21 +30,19 @@ public:
virtual ~MovableObjectInteractor();
void setPhysicsEngineGeometry(std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry);
- void setBlockVisitor(std::shared_ptr<BoundaryConditionsBlockVisitor> blockVisitor);
void moveGbObjectTo(const Vector3D& position);
private:
void rearrangeGrid();
void setSolidNodesToFluid();
+ void setBcNodesToFluid();
void reconstructDistributionOnSolidNodes();
- void setBcs() const;
void setBcBlocks();
void updateVelocityBc();
void setGeometryVelocityToBoundaryCondition(std::vector<int> node, std::shared_ptr<Block3D> block, std::shared_ptr<BCArray3D> bcArray) const;
- std::shared_ptr<BoundaryConditionsBlockVisitor> boundaryConditionsBlockVisitor;
std::shared_ptr<PhysicsEngineGeometryAdapter> physicsEngineGeometry;
std::shared_ptr<Reconstructor> reconstructor;
diff --git a/source/ThirdParty/Library/numerics/geometry3d/GbCuboid3D.cpp b/source/ThirdParty/Library/numerics/geometry3d/GbCuboid3D.cpp
index eae8f3315..35fff5d87 100644
--- a/source/ThirdParty/Library/numerics/geometry3d/GbCuboid3D.cpp
+++ b/source/ThirdParty/Library/numerics/geometry3d/GbCuboid3D.cpp
@@ -484,150 +484,150 @@ double GbCuboid3D::getCellVolumeInsideGbObject3D(const double& x1a,const double&
return 0.0;
}
/*==========================================================*/
-double GbCuboid3D::getIntersectionRaytraceFactor(const double& x1, const double& x2, const double& x3, const double& rx1, const double& rx2, const double& rx3)
-{
- double minB[3] = { this->getX1Minimum(), this->getX2Minimum(), this->getX3Minimum() };
- double maxB[3] = { this->getX1Maximum(), this->getX2Maximum(), this->getX3Maximum() };
- double origin[3] = { x1, x2, x3 }; //point
- double dir[3] = { rx1, rx2, rx3 }; //ray
-
- bool inside = true;
- char quadrant[3];
- int whichPlane;
- double maxT[3];
- double candidatePlane[3];
-
- /* Find candidate planes; this loop can be avoided if
- rays cast all from the eye(assume perpsective view) */
- for(int i=0; i<3; i++)
- {
- if(origin[i] < minB[i])
- {
- quadrant[i] = 1/*LEFT*/;
- candidatePlane[i] = minB[i];
- inside = false;
- }
- else if(origin[i] > maxB[i])
- {
- quadrant[i] = 0/*RIGHT*/;
- candidatePlane[i] = maxB[i];
- inside = false;
- }
- else
- {
- quadrant[i] = 2/*MIDDLE*/;
- }
- }
- /* Ray origin inside bounding box */
- if(inside)
- {
- //throw UbException(UB_EXARGS,"not done");
- return 0.0;
- }
-
- /* Calculate T distances to candidate planes */
- for(int i=0; i<3; i++)
- {
- if( quadrant[i]!=2/*MIDDLE*/ && fabs(dir[i])>1.E-10 )
- {
- maxT[i] = (candidatePlane[i]-origin[i])/dir[i];
- }
- else maxT[i] = -1.0;
- }
-
- /* Get largest of the maxT's for final choice of intersection */
- whichPlane = 0;
- for(int i=1; i<3; i++)
- if (maxT[whichPlane] < maxT[i])
- whichPlane = i;
-
- /* Check final candidate actually inside box */
- if(maxT[whichPlane]< -1.E-10) return -1.0;
- double dummy;
- for(int i= 0; i<3; i++)
- {
- if( whichPlane!= i)
- {
- dummy = origin[i] + maxT[whichPlane]*dir[i];
- if(dummy < minB[i] || dummy > maxB[i])
- return -1.0;
- }
- }
-
- return maxT[whichPlane] ; /* ray hits box */
-}
+//double GbCuboid3D::getIntersectionRaytraceFactor(const double& x1, const double& x2, const double& x3, const double& rx1, const double& rx2, const double& rx3)
+//{
+// double minB[3] = { this->getX1Minimum(), this->getX2Minimum(), this->getX3Minimum() };
+// double maxB[3] = { this->getX1Maximum(), this->getX2Maximum(), this->getX3Maximum() };
+// double origin[3] = { x1, x2, x3 }; //point
+// double dir[3] = { rx1, rx2, rx3 }; //ray
+//
+// bool inside = true;
+// char quadrant[3];
+// int whichPlane;
+// double maxT[3];
+// double candidatePlane[3];
+//
+// /* Find candidate planes; this loop can be avoided if
+// rays cast all from the eye(assume perpsective view) */
+// for(int i=0; i<3; i++)
+// {
+// if(origin[i] < minB[i])
+// {
+// quadrant[i] = 1/*LEFT*/;
+// candidatePlane[i] = minB[i];
+// inside = false;
+// }
+// else if(origin[i] > maxB[i])
+// {
+// quadrant[i] = 0/*RIGHT*/;
+// candidatePlane[i] = maxB[i];
+// inside = false;
+// }
+// else
+// {
+// quadrant[i] = 2/*MIDDLE*/;
+// }
+// }
+// /* Ray origin inside bounding box */
+// if(inside)
+// {
+// //throw UbException(UB_EXARGS,"not done");
+// return 0.0;
+// }
+//
+// /* Calculate T distances to candidate planes */
+// for(int i=0; i<3; i++)
+// {
+// if( quadrant[i]!=2/*MIDDLE*/ && fabs(dir[i])>1.E-10 )
+// {
+// maxT[i] = (candidatePlane[i]-origin[i])/dir[i];
+// }
+// else maxT[i] = -1.0;
+// }
+//
+// /* Get largest of the maxT's for final choice of intersection */
+// whichPlane = 0;
+// for(int i=1; i<3; i++)
+// if (maxT[whichPlane] < maxT[i])
+// whichPlane = i;
+//
+// /* Check final candidate actually inside box */
+// if(maxT[whichPlane]< -1.E-10) return -1.0;
+// double dummy;
+// for(int i= 0; i<3; i++)
+// {
+// if( whichPlane!= i)
+// {
+// dummy = origin[i] + maxT[whichPlane]*dir[i];
+// if(dummy < minB[i] || dummy > maxB[i])
+// return -1.0;
+// }
+// }
+//
+// return maxT[whichPlane] ; /* ray hits box */
+//}
// /*==========================================================*/
-// double GbCuboid3D::getIntersectionRaytraceFactor(const double& x1, const double& x2, const double& x3, const double& rx1, const double& rx2, const double& rx3)
-// {
-// double absX,absMaxX,absY,absMaxY,absZ,absMaxZ;
-//
-// if(rx1<0.0) absX = this->getX1Maximum() - x1;
-// else absX = this->getX1Minimum() - x1;
-// if(1-(rx1<0.0)) absMaxX = this->getX1Maximum() - x1;
-// else absMaxX = this->getX1Minimum() - x1;
-//
-// if(rx2<0.0) absY = this->getX2Maximum() - x2;
-// else absY = this->getX2Minimum() - x2;
-// if(1-(rx2<0.0)) absMaxY = this->getX2Maximum() - x2;
-// else absMaxY = this->getX2Minimum() - x2;
-//
-// if(rx3<0.0) absZ = this->getX3Maximum() - x3;
-// else absZ = this->getX3Minimum() - x3;
-// if(1-(rx3<0.0)) absMaxZ = this->getX3Maximum() - x3;
-// else absMaxZ = this->getX3Minimum() - x3;
-//
-//
-// //tmin ist die verschneidung des Gerade (Ray) durch die naehere Gerade (MinX oder MaxX)
-// //tmax ist die verschneidung des Gerade (Ray) durch die weiteste Gerade (MinX oder MaxX)
-// //analog fuer tymin und tymax
-// double tmin, tymin, tzmin, tmax, tymax, tzmax;
-//
-// if(!UbMath::zero(rx1)) tmin = tmax = 1.0/rx1;
-// else if(rx1<0.0) tmin = tmax = -UbMath::getPositiveInfinity<double>();
-// else tmin = tmax = UbMath::getPositiveInfinity<double>();
-//
-// if(!UbMath::zero(rx2)) tymin = tymax = 1.0/rx2;
-// else if(rx2<0.0) tymin = tymax = -UbMath::getPositiveInfinity<double>();
-// else tymin = tymax = UbMath::getPositiveInfinity<double>();
-//
-// if(!UbMath::zero(rx3)) tzmin = tzmax = 1.0/rx3;
-// else if(rx1<0.0) tzmin = tzmax = -UbMath::getPositiveInfinity<double>();
-// else tzmin = tzmax = UbMath::getPositiveInfinity<double>();
-//
-// //tmin *= absX;
-// //tmax *= absMaxX;
-// //tymin *= absY;
-// //tymax *= absMaxY;
-// //tzmin *= absZ;
-// //tzmax *= absMaxZ;
-//
-// //0 * 1/0 vermeiden, da es ein Undefined wert produziert
-// if( !UbMath::zero(absX) || !UbMath::zero(rx1) ) tmin *= absX;
-// else tmin = tymin;
-//
-// if( !UbMath::zero(absY) || !UbMath::zero(rx2)) tymin *= absY;
-// else tymin = tmin;
-//
-// if( !UbMath::zero(absZ) || !UbMath::zero(rx3)) tzmin *= absZ;
-// else tzmin = tymin;
-//
-// if( !UbMath::zero(absMaxX) || !UbMath::zero(rx1)) tmax *= absMaxX;
-// else tmax = tymax;
-//
-// if( !UbMath::zero(absMaxY) || !UbMath::zero(rx2)) tymax *= absMaxY;
-// else tymax = tmax;
-//
-// if( !UbMath::zero(absMaxZ) || !UbMath::zero(rx3)) tzmax *= absMaxZ;
-// else tzmax = tymax;
-//
-// //in dieser Fall gibt es keine Verschneidung
-// if( (tmin > tymax) || (tymin > tmax) ) return -1;
-//
-// tmin = UbMath::max3(tmin,tymin,tzmin);
-// tmax = UbMath::min3(tmax,tymax,tzmax);
-//
-// if( (tmin > tzmax) || (tzmin > tmax) ) return -1;
-// if(tmin >= 0.0) return tmin ;
-//
-// return tmax;
-//}
+ double GbCuboid3D::getIntersectionRaytraceFactor(const double& x1, const double& x2, const double& x3, const double& rx1, const double& rx2, const double& rx3)
+ {
+ double absX,absMaxX,absY,absMaxY,absZ,absMaxZ;
+
+ if(rx1<0.0) absX = this->getX1Maximum() - x1;
+ else absX = this->getX1Minimum() - x1;
+ if(1-(rx1<0.0)) absMaxX = this->getX1Maximum() - x1;
+ else absMaxX = this->getX1Minimum() - x1;
+
+ if(rx2<0.0) absY = this->getX2Maximum() - x2;
+ else absY = this->getX2Minimum() - x2;
+ if(1-(rx2<0.0)) absMaxY = this->getX2Maximum() - x2;
+ else absMaxY = this->getX2Minimum() - x2;
+
+ if(rx3<0.0) absZ = this->getX3Maximum() - x3;
+ else absZ = this->getX3Minimum() - x3;
+ if(1-(rx3<0.0)) absMaxZ = this->getX3Maximum() - x3;
+ else absMaxZ = this->getX3Minimum() - x3;
+
+
+ //tmin ist die verschneidung des Gerade (Ray) durch die naehere Gerade (MinX oder MaxX)
+ //tmax ist die verschneidung des Gerade (Ray) durch die weiteste Gerade (MinX oder MaxX)
+ //analog fuer tymin und tymax
+ double tmin, tymin, tzmin, tmax, tymax, tzmax;
+
+ if(!UbMath::zero(rx1)) tmin = tmax = 1.0/rx1;
+ else if(rx1<0.0) tmin = tmax = -UbMath::getPositiveInfinity<double>();
+ else tmin = tmax = UbMath::getPositiveInfinity<double>();
+
+ if(!UbMath::zero(rx2)) tymin = tymax = 1.0/rx2;
+ else if(rx2<0.0) tymin = tymax = -UbMath::getPositiveInfinity<double>();
+ else tymin = tymax = UbMath::getPositiveInfinity<double>();
+
+ if(!UbMath::zero(rx3)) tzmin = tzmax = 1.0/rx3;
+ else if(rx1<0.0) tzmin = tzmax = -UbMath::getPositiveInfinity<double>();
+ else tzmin = tzmax = UbMath::getPositiveInfinity<double>();
+
+ //tmin *= absX;
+ //tmax *= absMaxX;
+ //tymin *= absY;
+ //tymax *= absMaxY;
+ //tzmin *= absZ;
+ //tzmax *= absMaxZ;
+
+ //0 * 1/0 vermeiden, da es ein Undefined wert produziert
+ if( !UbMath::zero(absX) || !UbMath::zero(rx1) ) tmin *= absX;
+ else tmin = tymin;
+
+ if( !UbMath::zero(absY) || !UbMath::zero(rx2)) tymin *= absY;
+ else tymin = tmin;
+
+ if( !UbMath::zero(absZ) || !UbMath::zero(rx3)) tzmin *= absZ;
+ else tzmin = tymin;
+
+ if( !UbMath::zero(absMaxX) || !UbMath::zero(rx1)) tmax *= absMaxX;
+ else tmax = tymax;
+
+ if( !UbMath::zero(absMaxY) || !UbMath::zero(rx2)) tymax *= absMaxY;
+ else tymax = tmax;
+
+ if( !UbMath::zero(absMaxZ) || !UbMath::zero(rx3)) tzmax *= absMaxZ;
+ else tzmax = tymax;
+
+ //in dieser Fall gibt es keine Verschneidung
+ if( (tmin > tymax) || (tymin > tmax) ) return -1;
+
+ tmin = UbMath::max(tmin,tymin,tzmin);
+ tmax = UbMath::min(tmax,tymax,tzmax);
+
+ if( (tmin > tzmax) || (tzmin > tmax) ) return -1;
+ if(tmin >= 0.0) return tmin ;
+
+ return tmax;
+}
diff --git a/source/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp b/source/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
index 4fc7e64b9..d4b4df1d0 100644
--- a/source/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
+++ b/source/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
@@ -181,8 +181,8 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
data[0].push_back(3.0);
else if (bcArray->getBC(ix1, ix2, ix3)->hasSlipBoundary())
data[0].push_back(4.0);
- else
- data[0].push_back(5.0);
+ //else
+ // data[0].push_back(5.0);
if (bcArray->isSolid(ix1, ix2, ix3))
diff --git a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
index 285c9cdbd..2b08c7221 100644
--- a/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
+++ b/source/VirtualFluidsCore/Grid/BasicCalculator.cpp
@@ -46,8 +46,6 @@ void BasicCalculator::calculate()
for (calcStep = startTimeStep; calcStep<=numberOfTimeSteps; calcStep++)
{
- //coProcess((double)(calcStep-1));
-
//////////////////////////////////////////////////////////////////////////
#ifdef TIMING
UBLOG(logINFO, "calcStep = "<<calcStep);
@@ -124,7 +122,6 @@ void BasicCalculator::calculate()
}
}
//exchange data between blocks for visualization
- //if ((int)additionalGhostLayerUpdateScheduler->getNextDueTime()==calcStep)
if (additionalGhostLayerUpdateScheduler->isDue(calcStep))
{
exchangeBlockData(straightStartLevel, maxInitLevel);
--
GitLab