Merge remote-tracking branch 'origin/development/PrimitiveQs'

src/GridGenerator/geometries/Object.cu

@@ -6,3 +6,8 @@ void Object::findInnerNodes(SPtr<GridImp> grid)
 {
     grid->getGridStrategy()->findInnerNodes( grid );
 }
+
+HOST int Object::getIntersection(const Vertex & P, const Vertex & direction, Vertex & pointOnObject, real & qVal)
+{
+    return 1;
+}

src/GridGenerator/geometries/Object.h

@@ -12,6 +12,7 @@
 #include "global.h"
 
 class GridImp;
+struct Vertex;
 
 class VF_PUBLIC Object
 {
@@ -48,6 +49,8 @@ public:
     }
 
     HOST virtual void findInnerNodes(SPtr<GridImp> grid);
+
+    HOST virtual int getIntersection(const Vertex &P, const Vertex &direction, Vertex &pointOnObject, real &qVal);
 };
 
 

src/GridGenerator/geometries/Sphere/Sphere.cu

 #include "Sphere.h"
 
+#include <algorithm>    // std::min
+
+#include "geometries/Vertex/Vertex.h"
+
 Sphere::Sphere(const double& centerX, const double& centerY, const double& centerZ, const double& radius)
     : centerX(centerX), centerY(centerY), centerZ(centerZ), radius(radius)
 {
@@ -85,3 +89,49 @@ void Sphere::scale(double delta)
 {
     this->radius += delta;
 }
+
+HOST int Sphere::getIntersection(const Vertex & point, const Vertex & direction, Vertex & pointOnObject, real & qVal)
+{
+    
+    Vertex relativePoint( point.x - this->centerX, 
+                          point.y - this->centerY, 
+                          point.z - this->centerZ );
+
+    real directionSquare = direction.x * direction.x
+                         + direction.y * direction.y
+                         + direction.z * direction.z;
+
+    real p = 2* ( relativePoint.x * direction.x
+                + relativePoint.y * direction.y
+                + relativePoint.z * direction.z )
+                / directionSquare;
+
+    real q = ( relativePoint.x * relativePoint.x
+             + relativePoint.y * relativePoint.y
+             + relativePoint.z * relativePoint.z
+             - this->radius    * this->radius )
+           / directionSquare;
+
+    real discriminant = 0.25 * p * p - q;
+
+
+    if( discriminant < 0.0 ) return 1;
+
+    real result1 = - 0.5 * p + std::sqrt( discriminant );
+    real result2 = - 0.5 * p - std::sqrt( discriminant );
+
+    if( result1 < 0.0 && result2 < 0.0 ) return 1;
+
+    if( result1 < 0.0 ) result1 = 1e99;
+    if( result2 < 0.0 ) result2 = 1e99;
+
+    real t = std::min( result1, result2 );
+
+    pointOnObject.x = point.x + t * direction.x;
+    pointOnObject.y = point.y + t * direction.y;
+    pointOnObject.z = point.z + t * direction.z;
+
+    qVal = t;
+
+    return 0;
+}

src/GridGenerator/geometries/Sphere/Sphere.h

@@ -35,7 +35,10 @@ public:
 
 
     void scale(double delta) override;
-   
+    
+    HOST int getIntersection(const Vertex &P, const Vertex &direction, Vertex &pointOnObject, real &qVal) override;
+
+
 protected:
     double centerX;
     double centerY;

src/GridGenerator/grid/GridImp.cu

@@ -1225,6 +1225,8 @@ HOST void GridImp::findQs(Object* object) //TODO: enable qs for primitive object
     TriangularMesh* triangularMesh = dynamic_cast<TriangularMesh*>(object);
     if (triangularMesh)
         findQs(*triangularMesh);
+    else
+        findQsPrimitive(object);
 }
 
 HOST void GridImp::findQs(TriangularMesh &triangularMesh)
@@ -1286,6 +1288,114 @@ HOSTDEVICE void GridImp::findQs(Triangle &triangle)
     }
 }
 
+HOST void GridImp::findQsPrimitive(Object * object)
+{
+
+    if( this->qComputationStage == qComputationStageType::ComputeQs ){
+	    gridStrategy->allocateQs(shared_from_this());
+    }
+
+
+    for( int index = 0; index < this->size; index++ )
+    {
+
+        if( this->qIndices[index] == INVALID_INDEX ) continue;
+
+        real x,y,z;
+
+        this->transIndexToCoords(index,x,y,z);
+        
+        const Vertex point(x, y, z);
+
+        if( this->qComputationStage == qComputationStageType::ComputeQs ){
+            if(this->field.is(index, BC_SOLID))
+            {
+				calculateQs(index, point, object);
+			}
+        }
+        else if( this->qComputationStage == qComputationStageType::FindSolidBoundaryNodes )
+        {
+            if( !this->field.is(index, FLUID) ) continue;
+
+            if( checkIfAtLeastOneValidQ(index, point, object) )
+            {
+                // similar as in void GridImp::findBoundarySolidNode(uint index)
+                this->field.setFieldEntry( index, BC_SOLID );
+                this->qIndices[index] = this->numberOfSolidBoundaryNodes++;
+            }
+        }
+
+    }
+}
+
+HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, Object* object) const
+{
+    Vertex pointOnTriangle, direction;
+
+	real subdistance;
+	int error;
+	for (int i = distribution.dir_start; i <= distribution.dir_end; i++)
+	{
+		direction = Vertex( real(distribution.dirs[i * DIMENSION + 0]), 
+                            real(distribution.dirs[i * DIMENSION + 1]),
+			                real(distribution.dirs[i * DIMENSION + 2]) );
+
+		uint neighborIndex = this->transCoordToIndex(point.x + direction.x * this->delta,
+													    point.y + direction.y * this->delta,
+													    point.z + direction.z * this->delta);
+
+		if (neighborIndex == INVALID_INDEX) continue;
+
+		error = object->getIntersection(point, direction, pointOnTriangle, subdistance);
+
+		subdistance /= this->delta;
+
+		if (error == 0 && vf::Math::lessEqual(subdistance, 1.0) && vf::Math::greaterEqual(subdistance, 0.0))
+		{
+			if ( -0.5        > this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] ||
+                    subdistance < this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] )
+			{
+
+				this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = subdistance;
+                    
+                this->qPatches[ this->qIndices[index] ] = 0;
+
+				//printf("%d %f \n", this->qIndices[index], subdistance);
+			}
+		}
+	}
+}
+
+HOSTDEVICE bool GridImp::checkIfAtLeastOneValidQ(const uint index, const Vertex &point, Object* object) const
+{
+    Vertex pointOnTriangle, direction;
+
+	real subdistance;
+	int error;
+	for (int i = distribution.dir_start; i <= distribution.dir_end; i++)
+	{
+		direction = Vertex( real(distribution.dirs[i * DIMENSION + 0]), 
+                            real(distribution.dirs[i * DIMENSION + 1]),
+			                real(distribution.dirs[i * DIMENSION + 2]) );
+
+		uint neighborIndex = this->transCoordToIndex(point.x + direction.x * this->delta,
+													 point.y + direction.y * this->delta,
+													 point.z + direction.z * this->delta);
+
+		if (neighborIndex == INVALID_INDEX) continue;
+
+		error = object->getIntersection(point, direction, pointOnTriangle, subdistance);
+
+		subdistance /= this->delta;
+
+		if (error == 0 && vf::Math::lessEqual(subdistance, 1.0) && vf::Math::greaterEqual(subdistance, 0.0))
+		{
+			return true;
+		}
+	}
+    return false;
+}
+
 HOSTDEVICE void GridImp::setDebugPoint(uint index, int pointValue)
 {
     if (field.isInvalidCoarseUnderFine(index) && pointValue == INVALID_SOLID)
@@ -1329,7 +1439,6 @@ HOSTDEVICE void GridImp::calculateQs(const Vertex &point, const Triangle &triang
 HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, const Triangle &triangle) const
 {
 	Vertex pointOnTriangle, direction;
-	//VertexInteger solid_node;
 	real subdistance;
 	int error;
 	for (int i = distribution.dir_start; i <= distribution.dir_end; i++)
@@ -1348,14 +1457,6 @@ HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, cons
 
 		if (neighborIndex == INVALID_INDEX) continue;
 
-		//if ( this->field.is(neighborIndex, STOPPER_SOLID) )
-  //      {
-  //          if (this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] < -0.5)
-  //          {
-  //              this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = 1.0;
-  //          }
-  //      }
-
 		error = triangle.getTriangleIntersection(point, direction, pointOnTriangle, subdistance);
 
 		subdistance /= this->delta;
@@ -1378,7 +1479,6 @@ HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, cons
 HOSTDEVICE bool GridImp::checkIfAtLeastOneValidQ(const uint index, const Vertex & point, const Triangle & triangle) const
 {
 	Vertex pointOnTriangle, direction;
-	//VertexInteger solid_node;
 	real subdistance;
 	int error;
 	for (int i = distribution.dir_start; i <= distribution.dir_end; i++)
@@ -1386,8 +1486,9 @@ HOSTDEVICE bool GridImp::checkIfAtLeastOneValidQ(const uint index, const Vertex
 #if defined(__CUDA_ARCH__)
 		direction = Vertex(DIRECTIONS[i][0], DIRECTIONS[i][1], DIRECTIONS[i][2]);
 #else
-		direction = Vertex(real(distribution.dirs[i * DIMENSION + 0]), real(distribution.dirs[i * DIMENSION + 1]),
-			real(distribution.dirs[i * DIMENSION + 2]));
+		direction = Vertex(real(distribution.dirs[i * DIMENSION + 0]), 
+                           real(distribution.dirs[i * DIMENSION + 1]),
+			               real(distribution.dirs[i * DIMENSION + 2]));
 #endif
 
 		uint neighborIndex = this->transCoordToIndex(point.x + direction.x * this->delta,

src/GridGenerator/grid/GridImp.h

@@ -260,6 +260,8 @@ public:
     HOST void findQs(Object* object) override;
     HOST void findQs(TriangularMesh &triangularMesh);
     HOSTDEVICE void findQs(Triangle &triangle);
+
+    HOST void findQsPrimitive(Object* object);
 private:
 
     enum class qComputationStageType{
@@ -275,7 +277,11 @@ private:
     HOSTDEVICE void setDebugPoint(uint index, int pointValue);
 	HOSTDEVICE void calculateQs(const Vertex &point, const Triangle &triangle) const;
 	HOSTDEVICE void calculateQs(const uint index, const Vertex &point, const Triangle &triangle) const;
-    HOSTDEVICE bool checkIfAtLeastOneValidQ(const uint index, const Vertex &point, const Triangle &triangle) const;
+	HOST void calculateQs(const uint index, const Vertex &point, Object* object) const;
+
+    HOST bool checkIfAtLeastOneValidQ(const uint index, const Vertex &point, const Triangle &triangle) const;
+
+    HOST bool checkIfAtLeastOneValidQ(const uint index, const Vertex &point, Object* object) const;
 
 public:
 

targets/apps/LBM/gridGeneratorTest/main.cpp

-
-#define _USE_MATH_DEFINES
-#include <math.h>
-#include <string>
-#include <sstream>
-#include <iostream>
-#include <stdexcept>
-#include <fstream>
-#include <exception>
-#include <memory>
-
-//////////////////////////////////////////////////////////////////////////
-
-#include "Core/DataTypes.h"
-#include "Core/PointerDefinitions.h"
-
-#include "Core/LbmOrGks.h"
-#include "Core/Input/Input.h"
-#include "Core/StringUtilities/StringUtil.h"
-#include "Core/Input/ConfigFileReader/ConfigFileReader.h"
-
-#include "Core/VectorTypes.h"
-#include "Core/Logger/Logger.h"
-
-//////////////////////////////////////////////////////////////////////////
-
-#include "GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
-#include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
-#include "GridGenerator/grid/BoundaryConditions/Side.h"
-#include "GridGenerator/grid/GridFactory.h"
-
-#include "GridGenerator/geometries/Cuboid/Cuboid.h"
-#include "GridGenerator/geometries/TriangularMesh/TriangularMesh.h"
-#include "GridGenerator/geometries/Conglomerate/Conglomerate.h"
-#include "GridGenerator/geometries/TriangularMesh/TriangularMeshStrategy.h"
-
-#include "GridGenerator/io/SimulationFileWriter/SimulationFileWriter.h"
-#include "GridGenerator/io/GridVTKWriter/GridVTKWriter.h"
-#include "GridGenerator/io/STLReaderWriter/STLReader.h"
-#include "GridGenerator/io/STLReaderWriter/STLWriter.h"
-
-//////////////////////////////////////////////////////////////////////////
-
-#include "VirtualFluids_GPU/LBM/Simulation.h"
-#include "VirtualFluids_GPU/Communication/Communicator.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridProvider.h"
-#include "VirtualFluids_GPU/DataStructureInitializer/GridReaderFiles/GridReader.h"
-#include "VirtualFluids_GPU/Parameter/Parameter.h"
-#include "VirtualFluids_GPU/Output/FileWriter.h"
-
-#include "VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.h"
-#include "VirtualFluids_GPU/PreProcessor/PreProcessorFactory/PreProcessorFactoryImp.h"
-
-#include "VirtualFluids_GPU/GPU/CudaMemoryManager.h"
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//
-//          U s e r    s e t t i n g s
-//
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-//LbmOrGks lbmOrGks = GKS;
-//LbmOrGks lbmOrGks = LBM;
-
-const real L  = 2.0;
-
-const real Re = 100000.0;
-
-const real velocity  = 1.0;
-
-const real dt = 0.25e-3;
-
-const uint nx = 32;
-
-std::string path("F:/Work/Computations/out/gridGeneratorTest/"); //LEGOLAS
-//std::string path("E:/DrivenCavity/results/"); //TESLA03
-
-std::string simulationName("DrivenCavity");
-
-const uint timeStepOut = 10;
-const uint timeStepEnd = 1000000;
-
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-void multipleLevel(const std::string& configPath)
-{
-    logging::Logger::addStream(&std::cout);
-    logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
-    logging::Logger::timeStamp(logging::Logger::ENABLE);
-    logging::Logger::enablePrintedRankNumbers(logging::Logger::ENABLE);
-
-    auto gridFactory = GridFactory::make();
-    gridFactory->setGridStrategy(Device::CPU);
-    gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
-
-    auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
-    
-	Communicator* comm = Communicator::getInstanz();
-	SPtr<ConfigFileReader> configReader = ConfigFileReader::getNewInstance();
-	SPtr<ConfigData> configData = configReader->readConfigFile(configPath);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-	real dx = L / real(nx);
-
-	gridBuilder->addCoarseGrid(-L, -L, -L,
-							    L,  L,  L, dx);
-
-    TriangularMesh* cylinderSTL = TriangularMesh::make("F:/Work/Computations/gridGenerator/stl/Cylinder.stl");
-    
-
-    gridBuilder->setNumberOfLayers(6,6);
-    gridBuilder->addGrid(cylinderSTL, 1);
-
-    //Object* cube = new Cuboid(-0.5,-0.5,-0.5,0.5,0.5,0.5);
-    //gridBuilder->addGrid(cube, 1);
-
-    gridBuilder->addGeometry(cylinderSTL);
-
-	gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
-	gridBuilder->buildGrids(LBM, false); // buildGrids() has to be called before setting the BCs!!!!
-
-    gridBuilder->writeGridsToVtk("F:/Work/Computations/out/gridGeneratorTest/grid_");
-
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-	////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    SPtr<Parameter> para = Parameter::make(configData, comm);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    const real velocityLB = velocity * dt / dx; // LB units
-
-	const real vx = velocityLB; // LB units
-
-    const real viscosityLB = nx * velocityLB / Re; // LB units
-
-    *logging::out << logging::Logger::INFO_HIGH << "velocity  [dx/dt] = " << velocityLB << " \n";
-    *logging::out << logging::Logger::INFO_HIGH << "viscosity [dx/dt] = " << viscosityLB << "\n";
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    para->setOutputPath( path );
-    para->setOutputPrefix( simulationName );
-
-    para->setFName(para->getOutputPath() + "/" + para->getOutputPrefix());
-
-    para->setPrintFiles(true);
-    
-    para->setVelocity(velocityLB);
-    para->setViscosity(viscosityLB);
-
-    //para->setVelocityRatio(1.0 / velocityLB);
-    para->setVelocityRatio(1.0);
-
-    para->setTOut( timeStepOut );
-    para->setTEnd( timeStepEnd );
-
-    para->setUseWale(false);
-
-    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-	gridBuilder->setVelocityBoundaryCondition(SideType::MX, vx, 0.0, 0.0);
-	gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
-	//gridBuilder->setVelocityBoundaryCondition(SideType::PX, 0.0, 0.0, 0.0);
-    gridBuilder->setVelocityBoundaryCondition(SideType::PY, vx, 0.0, 0.0);
-    gridBuilder->setVelocityBoundaryCondition(SideType::MY, vx, 0.0, 0.0);
-    gridBuilder->setVelocityBoundaryCondition(SideType::PZ, vx, 0.0, 0.0);
-    gridBuilder->setVelocityBoundaryCondition(SideType::MZ, vx, 0.0, 0.0);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    SimulationFileWriter::write("F:/Work/Computations/out/gridGeneratorTest/grid/", gridBuilder, FILEFORMAT::BINARY);
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-    SPtr<CudaMemoryManager> cudaMemoryManager = CudaMemoryManager::make(para);
-    SPtr<GridProvider> gridGenerator = GridProvider::makeGridGenerator(gridBuilder, para, cudaMemoryManager);
-    //SPtr<GridProvider> gridGenerator = GridProvider::makeGridReader(FILEFORMAT::BINARY, para, cudaMemoryManager);
-
-    Simulation sim;
-    SPtr<FileWriter> fileWriter = SPtr<FileWriter>(new FileWriter());
-    SPtr<KernelFactoryImp> kernelFactory = KernelFactoryImp::getInstance();
-    SPtr<PreProcessorFactoryImp> preProcessorFactory = PreProcessorFactoryImp::getInstance();
-    sim.setFactories(kernelFactory, preProcessorFactory);
-    sim.init(para, gridGenerator, fileWriter, cudaMemoryManager);
-    
-
-    sim.run();
-    sim.free();
-
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-}
-
-int main( int argc, char* argv[])
-{
-    MPI_Init(&argc, &argv);
-    std::string str, str2; 
-    if ( argv != NULL )
-    {
-        //str = static_cast<std::string>(argv[0]);
-        
-        try
-        {
-            //////////////////////////////////////////////////////////////////////////
-
-			std::string targetPath;
-
-			targetPath = __FILE__;
-
-#ifdef _WIN32
-			targetPath = targetPath.substr(0, targetPath.find_last_of('\\') + 1);
-#else
-			targetPath = targetPath.substr(0, targetPath.find_last_of('/') + 1);
-#endif
-
-			std::cout << targetPath << std::endl;
-
-			multipleLevel(targetPath + "config.txt");
-
-            //////////////////////////////////////////////////////////////////////////
-		}
-        catch (const std::exception& e)
-        {
-                
-            *logging::out << logging::Logger::LOGGER_ERROR << e.what() << "\n";
-        }
-        catch (const std::bad_alloc e)
-        {
-                
-            *logging::out << logging::Logger::LOGGER_ERROR << "Bad Alloc:" << e.what() << "\n";
-        }
-        catch (...)
-        {
-            *logging::out << logging::Logger::LOGGER_ERROR << "Unknown exception!\n";
-        }
-    }
-
-   MPI_Finalize();
-   return 0;
-}