From a174cfdd31b106ae5743e03791ef322763386265 Mon Sep 17 00:00:00 2001
From: "LEGOLAS\\lenz" <lenz@irmb.tu-bs.de>
Date: Fri, 16 Nov 2018 10:43:25 +0100
Subject: [PATCH] implements GksMeshAdapter in 3D
---
CMakeLists.txt | 11 +-
src/GksMeshAdapter/GksMeshAdapter.cpp | 1525 +++++++++++++++++
src/GksMeshAdapter/GksMeshAdapter.h | 151 ++
src/GksMeshAdapter/MeshCell.cpp | 32 +
src/GksMeshAdapter/MeshCell.h | 71 +
src/GksMeshAdapter/MeshFace.cpp | 12 +
src/GksMeshAdapter/MeshFace.h | 50 +
src/GksMeshAdapter/package.include | 0
targets/apps/GKS/gksTest/CMakeLists.txt | 6 +-
targets/apps/GKS/gksTest/main.cpp | 1161 +------------
.../libs/GksMeshAdapter/3rdPartyLinking.cmake | 6 +
targets/libs/GksMeshAdapter/CMakeLists.txt | 18 +
.../libs/GksMeshAdapter/CMakePackage.cmake | 8 +
13 files changed, 1939 insertions(+), 1112 deletions(-)
create mode 100644 src/GksMeshAdapter/GksMeshAdapter.cpp
create mode 100644 src/GksMeshAdapter/GksMeshAdapter.h
create mode 100644 src/GksMeshAdapter/MeshCell.cpp
create mode 100644 src/GksMeshAdapter/MeshCell.h
create mode 100644 src/GksMeshAdapter/MeshFace.cpp
create mode 100644 src/GksMeshAdapter/MeshFace.h
create mode 100644 src/GksMeshAdapter/package.include
create mode 100644 targets/libs/GksMeshAdapter/3rdPartyLinking.cmake
create mode 100644 targets/libs/GksMeshAdapter/CMakeLists.txt
create mode 100644 targets/libs/GksMeshAdapter/CMakePackage.cmake
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 02360bd20..c40dec06b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -52,11 +52,16 @@ include_directories(${CMAKE_BINARY_DIR})
set_property(GLOBAL PROPERTY USE_FOLDERS ON)
set_property(GLOBAL PROPERTY PREDEFINED_TARGETS_FOLDER ".cmake")
-set(libraryFolder "libs")
+
+set(libraryFolder "libs")
+set(gksLibraryFolder "libs/GKS")
+
set(testFolder "tests")
-set(appFolder "apps")
+
+set(appFolder "apps")
set(lbmAppFolder "apps/LBM")
set(gksAppFolder "apps/GKS")
+
set(thirdPartyFolder "3rdParty")
IF(MSVC)
@@ -123,6 +128,8 @@ add_subdirectory(3rdParty/metis/metis-5.1.0)
add_subdirectory(targets/libs/VirtualFluidsBasics)
add_subdirectory(targets/libs/Core)
+add_subdirectory(targets/libs/GksMeshAdapter)
+
if(HULC.BUILD_NUMERIC_TESTS)
add_subdirectory(3rdParty/fftw/fftw-3.3.7)
add_subdirectory(targets/tests/TestingHULC)
diff --git a/src/GksMeshAdapter/GksMeshAdapter.cpp b/src/GksMeshAdapter/GksMeshAdapter.cpp
new file mode 100644
index 000000000..e67b3fa13
--- /dev/null
+++ b/src/GksMeshAdapter/GksMeshAdapter.cpp
@@ -0,0 +1,1525 @@
+#include "GksMeshAdapter.h"
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+#include <fstream>
+#include <algorithm>
+#include <numeric>
+#include <functional>
+#include <iostream>
+#include <iomanip>
+
+#include "core/Logger/Logger.h"
+
+#include "GridGenerator/grid/distributions/D3Q27.h"
+#include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
+#include "GridGenerator/grid/NodeValues.h"
+#include "GridGenerator/utilities/math/Math.h"
+
+#include "MeshCell.h"
+#include "MeshFace.h"
+
+void GksMeshAdapter::inputGrid(SPtr<MultipleGridBuilder> gridBuilder)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "inputGrid()" << "\n";
+
+ this->numberOfLevels = gridBuilder->getNumberOfGridLevels();
+
+ std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ this->dxCoarse = grids[0]->getDelta();
+
+ //////////////////////////////////////////////////////////////////////////
+
+ this->gridToMesh.resize( gridBuilder->getNumberOfGridLevels() );
+
+ for( uint level = 0; level < gridBuilder->getNumberOfGridLevels(); level++ ){
+ this->gridToMesh[level].resize( grids[level]->getSize() );
+
+ for( auto& cellIdx : this->gridToMesh[level] ) cellIdx = INVALID_INDEX;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // I d e n t i f y C e l l s i n L B - G r i d
+ //
+ //////////////////////////////////////////////////////////////////////////
+
+ uint numberOfCells = 0;
+
+ for( uint level = 0; level < gridBuilder->getNumberOfGridLevels(); level++ ){
+ for( uint gridIdx = 0; gridIdx < grids[level]->getSize(); gridIdx++ ){
+ if (grids[level]->getFieldEntry(gridIdx) != STOPPER_COARSE_UNDER_FINE &&
+ grids[level]->getFieldEntry(gridIdx) != STOPPER_SOLID &&
+ grids[level]->getFieldEntry(gridIdx) != INVALID_COARSE_UNDER_FINE &&
+ grids[level]->getFieldEntry(gridIdx) != INVALID_OUT_OF_GRID &&
+ grids[level]->getFieldEntry(gridIdx) != INVALID_SOLID )
+ {
+ this->gridToMesh[level][gridIdx] = numberOfCells++;
+ }
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // S e t M e s h t o G r i d i n f o r m a t i o n
+ //
+ //////////////////////////////////////////////////////////////////////////
+
+ this->cells.resize( numberOfCells );
+
+ for( uint level = 0; level < gridBuilder->getNumberOfGridLevels(); level++ ){
+ for( uint gridIdx = 0; gridIdx < grids[level]->getSize(); gridIdx++ ){
+ if ( this->gridToMesh[level][gridIdx] != INVALID_INDEX ){
+
+ uint cellIdx = gridToMesh[level][gridIdx];
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ cell.level = level;
+ cell.gridIdx = gridIdx;
+
+ cell.type = grids[level]->getFieldEntry(gridIdx);
+ }
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+}
+
+void GksMeshAdapter::findQuadtreeConnectivity(SPtr<MultipleGridBuilder> gridBuilder)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "findQuadtreeConnectivity()" << "\n";
+
+ std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ if( cell.type == FLUID_FCC || cell.type == FLUID_CFC ){
+
+ real x, y, z;
+ grids[cell.level]->transIndexToCoords(cell.gridIdx, x, y, z);
+
+ real d = 0.25 * grids[cell.level]->getDelta();
+
+ for( uint idx = 0; idx < 8; idx++ )
+ {
+ Distribution dirs = DistributionHelper::getDistribution27();
+
+ real xSign = dirs.directions[idx + 19][0];
+ real ySign = dirs.directions[idx + 19][1];
+ real zSign = dirs.directions[idx + 19][2];
+
+ cell.children[ idx ] = this->gridToMesh[cell.level+1][ grids[cell.level+1]->transCoordToIndex( x + xSign * d,
+ y + ySign * d,
+ z + zSign * d ) ];
+ }
+
+ // register parent
+ for( uint child = 0; child < 8; child++ )
+ this->cells[ cell.children[child] ].parent = cellIdx;
+
+ // set correct type for CFF cells
+ for( uint child = 0; child < 8; child++ )
+ if( this->cells[ cell.children[child] ].type != FLUID_FCF )
+ this->cells[ cell.children[child] ].type = FLUID_CFF;
+
+ }
+ }
+}
+
+void GksMeshAdapter::findCellToCellConnectivity(SPtr<MultipleGridBuilder> gridBuilder)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "findCellToCellConnectivity()" << "\n";
+
+ std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ real x, y, z;
+ grids[cell.level]->transIndexToCoords(cell.gridIdx, x, y, z);
+
+ real d = grids[cell.level]->getDelta();
+
+ for( uint idx = 0; idx < 27; idx++ )
+ {
+ if( idx == DIR_27_ZERO ) continue;
+
+ Distribution dirs = DistributionHelper::getDistribution27();
+
+ int xSign = dirs.directions[idx][0];
+ int ySign = dirs.directions[idx][1];
+ int zSign = dirs.directions[idx][2];
+
+ uint neighborGridIdx = grids[cell.level]->transCoordToIndex( x + xSign * d,
+ y + ySign * d,
+ z + zSign * d );
+
+ if( neighborGridIdx == INVALID_INDEX || this->gridToMesh[cell.level][neighborGridIdx] == INVALID_INDEX ){
+ if( !cell.isCoarseGhostCell() && cell.type != BC_SOLID )
+ cell.isGhostCell = true;
+
+ continue;
+ }
+
+ cell.cellToCell[ idx ] = this->gridToMesh[cell.level][neighborGridIdx];
+ }
+ }
+}
+
+void GksMeshAdapter::countCells()
+{
+ this->numberOfCellsPerLevel .resize( this->numberOfLevels );
+ this->numberOfBulkCellsPerLevel.resize( this->numberOfLevels );
+ this->startOfCellsPerLevel .resize( this->numberOfLevels );
+
+ for( auto& i : this->numberOfCellsPerLevel ) i = 0;
+ for( auto& i : this->numberOfBulkCellsPerLevel ) i = 0;
+ for( auto& i : this->startOfCellsPerLevel ) i = 0;
+
+ uint level = 0;
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ if( cell.level != level ) level++;
+
+ this->numberOfCellsPerLevel[ level ]++;
+
+ if( ! ( cell.isGhostCell || cell.isCoarseGhostCell() ) )
+ this->numberOfBulkCellsPerLevel[ level ]++;
+ }
+
+ for( uint level = 1; level < this->numberOfLevels; level++ )
+ this->startOfCellsPerLevel[ level ] = this->startOfCellsPerLevel[ level-1 ] + this->numberOfCellsPerLevel[ level-1 ];
+}
+
+void GksMeshAdapter::partitionCells()
+{
+ for( uint level = 0; level < this->numberOfLevels; level++ ){
+
+ std::vector<uint> idxMap( this->cells.size() );
+ std::iota( idxMap.begin(), idxMap.end(), 0 );
+
+ // partition idxMap
+ std::stable_partition( idxMap.begin() + this->startOfCellsPerLevel[level],
+ idxMap.begin() + this->startOfCellsPerLevel[level]
+ + this->numberOfCellsPerLevel[level],
+ [this](int lhs){
+ return ! ( this->cells[ lhs ].isGhostCell || this->cells[ lhs ].isCoarseGhostCell() );
+ }
+ );
+
+ // invert idxMap
+ {
+ std::vector<uint> buffer = idxMap;
+ for( uint idx = 0; idx < idxMap.size(); idx ++ )
+ idxMap[ buffer[idx] ] = idx;
+ }
+
+ // partition cell list
+ std::stable_partition( this->cells.begin() + this->startOfCellsPerLevel[level],
+ this->cells.begin() + this->startOfCellsPerLevel[level]
+ + this->numberOfCellsPerLevel[level],
+ [this](MeshCell lhs){
+ return ! ( lhs.isGhostCell || lhs.isCoarseGhostCell() );
+ }
+ );
+
+ this->refreshCellConnectivity( idxMap );
+ }
+}
+
+void GksMeshAdapter::refreshCellConnectivity(const std::vector<uint>& idxMap)
+{
+ for( auto& cell : this->cells ){
+ for( uint idx = 0; idx < 27; idx++ )
+ if( cell.cellToCell[ idx ] != INVALID_INDEX )
+ cell.cellToCell[ idx ] = idxMap[ cell.cellToCell[ idx ] ];
+
+ if( cell.parent != INVALID_INDEX )
+ cell.parent = idxMap[ cell.parent ];
+
+ for( uint idx = 0; idx < 8; idx++ )
+ if( cell.children[ idx ] != INVALID_INDEX )
+ cell.children[ idx ] = idxMap[ cell.children[ idx ] ];
+ }
+
+ for( auto& grid : this->gridToMesh ){
+ for( auto& cellIdx : grid ){
+ if( cellIdx != INVALID_INDEX )
+ cellIdx = idxMap[ cellIdx ];
+ }
+ }
+}
+
+void GksMeshAdapter::findCornerCells(SPtr<MultipleGridBuilder> gridBuilder, real z0)
+{
+ //SPtr<Grid> grid = gridBuilder->getGrids()[0];
+ //
+ //this->cornerCells[0] = this->gridToMesh[ 0 ][ grid->transCoordToIndex( grid->getStartX(), grid->getStartY(), z0 ) ];
+ //this->cornerCells[1] = this->gridToMesh[ 0 ][ grid->transCoordToIndex( grid->getEndX() , grid->getStartY(), z0 ) ];
+ //this->cornerCells[2] = this->gridToMesh[ 0 ][ grid->transCoordToIndex( grid->getEndX() , grid->getEndY() , z0 ) ];
+ //this->cornerCells[3] = this->gridToMesh[ 0 ][ grid->transCoordToIndex( grid->getStartX(), grid->getEndY() , z0 ) ];
+}
+
+void GksMeshAdapter::generateNodes(SPtr<MultipleGridBuilder> gridBuilder)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "generateNodes()" << "\n";
+
+ std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ nodes.reserve( 2 * this->cells.size() );
+
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ if( cell.type == STOPPER_SOLID ) continue;
+
+ real x, y, z;
+ grids[cell.level]->transIndexToCoords(cell.gridIdx, x, y, z);
+
+ real d = 0.5 * grids[cell.level]->getDelta();
+
+ std::array<Vec3,8> dir;
+
+ for( uint idx = 0; idx < 8; idx++ )
+ {
+ if( cell.cellToNode[idx] == INVALID_INDEX )
+ {
+ Distribution dirs = DistributionHelper::getDistribution27();
+
+ real dx = dirs.directions[idx + 19][0] * d;
+ real dy = dirs.directions[idx + 19][1] * d;
+ real dz = dirs.directions[idx + 19][2] * d;
+
+ nodes.push_back( Vec3( x + dx, y + dy, z + dz ) );
+
+ cell.cellToNode[idx] = nodes.size()-1;
+
+ // register new node at neighbor cells on same level
+ for( uint idx = 0; idx < 8; idx++ )
+ {
+ if( cell.cellToCell[idx] == INVALID_INDEX ) continue;
+
+ Distribution dirs = DistributionHelper::getDistribution27();
+
+ real dxNeighbor = - dirs.directions[idx + 19][0] * d;
+ real dyNeighbor = - dirs.directions[idx + 19][1] * d;
+ real dzNeighbor = - dirs.directions[idx + 19][2] * d;
+
+ real xNeighbor = nodes.back().x + dx + dxNeighbor;
+ real yNeighbor = nodes.back().y + dy + dyNeighbor;
+ real zNeighbor = nodes.back().z + dz + dzNeighbor;
+
+ uint neighborGridIdx = grids[cell.level]->transCoordToIndex( xNeighbor, yNeighbor, zNeighbor );
+
+ if( neighborGridIdx != INVALID_INDEX && gridToMesh[ cell.level ][ neighborGridIdx ] == cell.cellToCell[idx] )
+ {
+ this->cells[ cell.cellToCell[idx] ].cellToNode[ idx ];
+ }
+ }
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ //if( cell.type == FLUID_FCC || cell.type == FLUID_CFC ){
+
+ // // register edge nodes at child cells
+ // this->cells[ cell.children[0] ].cellToNode[0] = cell.cellToNode[0];
+ // this->cells[ cell.children[1] ].cellToNode[1] = cell.cellToNode[1];
+ // this->cells[ cell.children[2] ].cellToNode[2] = cell.cellToNode[2];
+ // this->cells[ cell.children[3] ].cellToNode[3] = cell.cellToNode[3];
+
+ // std::array<arraytypes::Vec2,4> dir;
+
+ // dir[0] = arraytypes::Vec2( - d, 0.0 );
+ // dir[1] = arraytypes::Vec2( 0.0, - d );
+ // dir[2] = arraytypes::Vec2( d, 0.0 );
+ // dir[3] = arraytypes::Vec2( 0.0, d );
+
+ // // introduce new face nodes
+ // if( cell.cellToFaceNode[0] == -1 ){
+ // nodes.push_back( arraytypes::Vec2( x + dir[0].x, y + dir[0].y ) );
+
+ // cell.cellToFaceNode[0] = nodes.size()-1;
+ //
+ // // register face node at adjacent cell
+ // if( cell.cellToCell[0] != -1 ) this->cells[ cell.cellToCell[0] ].cellToFaceNode[2] = nodes.size()-1;
+
+ // }
+ //
+ // if( cell.cellToFaceNode[1] == -1 ){
+ // nodes.push_back( arraytypes::Vec2( x + dir[1].x, y + dir[1].y ) );
+
+ // cell.cellToFaceNode[1] = nodes.size()-1;
+ //
+ // // register face node at adjacent cell
+ // if( cell.cellToCell[1] != -1 ) this->cells[ cell.cellToCell[1] ].cellToFaceNode[3] = nodes.size()-1;
+
+ // }
+ //
+ // if( cell.cellToFaceNode[2] == -1 ){
+ // nodes.push_back( arraytypes::Vec2( x + dir[2].x, y + dir[2].y ) );
+
+ // cell.cellToFaceNode[2] = nodes.size()-1;
+ //
+ // // register face node at adjacent cell
+ // if( cell.cellToCell[2] != -1 ) this->cells[ cell.cellToCell[2] ].cellToFaceNode[0] = nodes.size()-1;
+
+ // }
+ //
+ // if( cell.cellToFaceNode[3] == -1 ){
+ // nodes.push_back( arraytypes::Vec2( x + dir[3].x, y + dir[3].y ) );
+
+ // cell.cellToFaceNode[3] = nodes.size()-1;
+ //
+ // // register face node at adjacent cell
+ // if( cell.cellToCell[3] != -1 ) this->cells[ cell.cellToCell[3] ].cellToFaceNode[1] = nodes.size()-1;
+
+ // }
+
+ // // register face node at child cells
+ // this->cells[ cell.children[3] ].cellToNode[0] = cell.cellToFaceNode[0];
+ // this->cells[ cell.children[0] ].cellToNode[3] = cell.cellToFaceNode[0];
+
+ // this->cells[ cell.children[0] ].cellToNode[1] = cell.cellToFaceNode[1];
+ // this->cells[ cell.children[1] ].cellToNode[0] = cell.cellToFaceNode[1];
+
+ // this->cells[ cell.children[1] ].cellToNode[2] = cell.cellToFaceNode[2];
+ // this->cells[ cell.children[2] ].cellToNode[1] = cell.cellToFaceNode[2];
+
+ // this->cells[ cell.children[2] ].cellToNode[3] = cell.cellToFaceNode[3];
+ // this->cells[ cell.children[3] ].cellToNode[2] = cell.cellToFaceNode[3];
+ //}
+ }
+}
+
+void GksMeshAdapter::inputQs(SPtr<MultipleGridBuilder> gridBuilder)
+{
+
+ //std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( cell.type != BC_SOLID ) continue;
+
+ // cell.Qs[0] = grids[cell.level]->getQValue(cell.gridIdx, 7);
+ // cell.Qs[1] = grids[cell.level]->getQValue(cell.gridIdx, 8);
+ // cell.Qs[2] = grids[cell.level]->getQValue(cell.gridIdx, 6);
+ // cell.Qs[3] = grids[cell.level]->getQValue(cell.gridIdx, 9);
+ //}
+}
+
+void GksMeshAdapter::morphNodes()
+{
+ //this->displacement.resize( this->nodes.size() );
+ //this->nDisplacements.resize( this->nodes.size() );
+
+ //for( uint n : this->nDisplacements ) n = 0;
+
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( cell.type != BC_SOLID ) continue;
+
+ // std::array<arraytypes::Vec2,4> dir;
+ // dir[0] = arraytypes::Vec2( - getDx(cell.level), - getDx(cell.level) );
+ // dir[1] = arraytypes::Vec2( getDx(cell.level), - getDx(cell.level) );
+ // dir[2] = arraytypes::Vec2( getDx(cell.level), getDx(cell.level) );
+ // dir[3] = arraytypes::Vec2( - getDx(cell.level), getDx(cell.level) );
+
+ // for( uint idx = 0; idx < 4; idx ++ ){
+ //
+ // //if( vf::Math::equal( cell.Qs[idx], 0.0 ) ) continue;
+ // if( cell.Qs[idx] < 0.0 || cell.Qs[idx] > 1.0 ) continue;
+
+ // this->displacement[ cell.cellToNode[idx] ].x += ( cell.Qs[idx] - 0.5 ) * dir[idx].x;
+ // this->displacement[ cell.cellToNode[idx] ].y += ( cell.Qs[idx] - 0.5 ) * dir[idx].y;
+ //
+ // this->nDisplacements[ cell.cellToNode[idx] ]++;
+ // }
+ //}
+
+ //for( uint nodeIdx = 0; nodeIdx < this->nodes.size(); nodeIdx++ ){
+ //
+ // if( this->nDisplacements[nodeIdx] == 0 ) continue;
+ //
+ // nodes[ nodeIdx ].x += this->displacement[ nodeIdx ].x / double( this->nDisplacements[ nodeIdx ] );
+ // nodes[ nodeIdx ].y += this->displacement[ nodeIdx ].y / double( this->nDisplacements[ nodeIdx ] );
+ //}
+}
+
+void GksMeshAdapter::repairTriangularCells()
+{
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( cell.type != BC_SOLID ) continue;
+
+ // arraytypes::uint4 cellToNodeTMP = cell.cellToNode;
+
+ // for( uint idx = 0; idx < 4; idx++ ){
+ //
+ // arraytypes::Vec2 v1, v2;
+
+ // v1.x = this->nodes[ cell.cellToNode[ idx ] ].x - this->nodes[ cell.cellToNode[ ( idx - 1 ) % 4 ] ].x;
+ // v1.y = this->nodes[ cell.cellToNode[ idx ] ].y - this->nodes[ cell.cellToNode[ ( idx - 1 ) % 4 ] ].y;
+ // v2.x = - this->nodes[ cell.cellToNode[ idx ] ].x + this->nodes[ cell.cellToNode[ ( idx + 1 ) % 4 ] ].x;
+ // v2.y = - this->nodes[ cell.cellToNode[ idx ] ].y + this->nodes[ cell.cellToNode[ ( idx + 1 ) % 4 ] ].y;
+
+ // double v1Length = sqrt( v1.x*v1.x + v1.y*v1.y );
+ // double v2Length = sqrt( v2.x*v2.x + v2.y*v2.y );
+
+ // double v1DotV2 = v1.x*v2.x + v1.y*v2.y;
+
+ // bool isConcave = v1.x*v2.y - v1.y*v2.x < 0.0;
+
+ // double argument = v1DotV2 / ( v1Length * v2Length );
+
+ // bool isCollinear = false;
+ //
+ // // the acos function does not like the argument 1.0, hence the special treatment
+ // if( fabs( v1DotV2 / ( v1Length * v2Length ) - 1.0 ) < 1.0e-10 )
+ // isCollinear = true;
+ // else
+ // isCollinear = acos( v1DotV2 / ( v1Length * v2Length ) ) < M_PI / 4.0;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // uint neighborIdx1 = cell.cellToCell[ idx ];
+ // uint neighborIdx2 = cell.cellToCell[ ( idx + 1 ) % 4 ];
+
+ // bool neighborCellsExist = ( ( neighborIdx1 != INVALID_IDX && this->cells[ neighborIdx1 ].type != STOPPER_SOLID ) ||
+ // ( neighborIdx2 != INVALID_IDX && this->cells[ neighborIdx2 ].type != STOPPER_SOLID ) );
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // if( isConcave || isCollinear ){
+ // if( !neighborCellsExist ){
+ // cellToNodeTMP[idx] = INVALID_IDX;
+ // }
+ // else{
+ // nodes[ cell.cellToNode[ idx ] ].x -= this->displacement[ cell.cellToNode[ idx ] ].x / double( this->nDisplacements[ cell.cellToNode[ idx ] ] );
+ // nodes[ cell.cellToNode[ idx ] ].y -= this->displacement[ cell.cellToNode[ idx ] ].y / double( this->nDisplacements[ cell.cellToNode[ idx ] ] );
+ // }
+ // }
+ // }
+
+ // cell.cellToNode = cellToNodeTMP;
+ //}
+}
+
+void GksMeshAdapter::computeCellGeometry()
+{
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ Vec3 cellCenter;
+
+ for( uint node = 0; node < 8; node++ ){
+ cellCenter = cellCenter + this->nodes[ cell.cellToNode[node] ];
+ }
+
+ cell.cellCenter.x = cellCenter.x / eight;
+ cell.cellCenter.y = cellCenter.y / eight;
+ cell.cellCenter.z = cellCenter.z / eight;
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // uint nNodes = 0;
+ // for( uint idx : cell.cellToNode ) if( idx != INVALID_IDX ) nNodes++;
+
+ // if( nNodes == 3 ){
+ //
+ // arraytypes::uint3 cellToNode;
+ // for( uint idx = 0, counter = 0; idx < 4; idx++ )
+ // if( cell.cellToNode[idx] != INVALID_IDX )
+ // cellToNode[counter++] = cell.cellToNode[idx];
+ //
+ // cell.cellCenter.x = ( this->nodes[ cellToNode[0] ].x
+ // + this->nodes[ cellToNode[1] ].x
+ // + this->nodes[ cellToNode[2] ].x ) / 3.0;
+ // cell.cellCenter.y = ( this->nodes[ cellToNode[0] ].y
+ // + this->nodes[ cellToNode[1] ].y
+ // + this->nodes[ cellToNode[2] ].y ) / 3.0;
+
+ // cell.cellVolume = 0.5 * ( this->nodes[ cellToNode[0] ].x * ( this->nodes[ cellToNode[1] ].y - this->nodes[ cellToNode[2] ].y )
+ // + this->nodes[ cellToNode[1] ].x * ( this->nodes[ cellToNode[2] ].y - this->nodes[ cellToNode[0] ].y )
+ // + this->nodes[ cellToNode[2] ].x * ( this->nodes[ cellToNode[0] ].y - this->nodes[ cellToNode[1] ].y ) );
+
+ // }
+ // else if( nNodes == 4 ){
+ //
+ // arraytypes::Vec2 triCenter[2];
+ // triCenter[0].x = (this->nodes[ cell.cellToNode[0] ].x + this->nodes[ cell.cellToNode[1] ].x + this->nodes[ cell.cellToNode[3] ].x) / 3.0;
+ // triCenter[0].y = (this->nodes[ cell.cellToNode[0] ].y + this->nodes[ cell.cellToNode[1] ].y + this->nodes[ cell.cellToNode[3] ].y) / 3.0;
+ // triCenter[1].y = ( this->nodes[ cell.cellToNode[1] ].y + this->nodes[ cell.cellToNode[2] ].y + this->nodes[ cell.cellToNode[3] ].y) / 3.0;
+ // triCenter[1].x = ( this->nodes[ cell.cellToNode[1] ].x + this->nodes[ cell.cellToNode[2] ].x + this->nodes[ cell.cellToNode[3] ].x) / 3.0;
+
+ // double triVolume[2];
+ // triVolume[0] = 0.5 * fabs( this->nodes[ cell.cellToNode[0] ].x * ( this->nodes[ cell.cellToNode[1] ].y - this->nodes[ cell.cellToNode[3] ].y )
+ // + this->nodes[ cell.cellToNode[1] ].x * ( this->nodes[ cell.cellToNode[3] ].y - this->nodes[ cell.cellToNode[0] ].y )
+ // + this->nodes[ cell.cellToNode[3] ].x * ( this->nodes[ cell.cellToNode[0] ].y - this->nodes[ cell.cellToNode[1] ].y ) );
+ // triVolume[1] = 0.5 * fabs( this->nodes[ cell.cellToNode[2] ].x * ( this->nodes[ cell.cellToNode[3] ].y - this->nodes[ cell.cellToNode[1] ].y )
+ // + this->nodes[ cell.cellToNode[3] ].x * ( this->nodes[ cell.cellToNode[1] ].y - this->nodes[ cell.cellToNode[2] ].y )
+ // + this->nodes[ cell.cellToNode[1] ].x * ( this->nodes[ cell.cellToNode[2] ].y - this->nodes[ cell.cellToNode[3] ].y ) );
+
+ // cell.cellVolume = triVolume[0] + triVolume[1];
+ // cell.cellCenter.x = ( triCenter[0].x * triVolume[0] + triCenter[1].x * triVolume[1] ) / cell.cellVolume;
+ // cell.cellCenter.y = ( triCenter[0].y * triVolume[0] + triCenter[1].y * triVolume[1] ) / cell.cellVolume;
+ // }
+ }
+}
+
+void GksMeshAdapter::generateSolidGhostCells()
+{
+ //std::vector<MeshCell> newCells;
+
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( cell.type != BC_SOLID ) continue;
+ //
+ // for( uint idx = 0; idx < 4; idx++ ){
+ //
+ // if( cell.cellToCell[ idx ] == INVALID_IDX ){
+ //
+ // cell.cellToCell[ idx ] = this->cells.size() + newCells.size();
+
+ // if( cell.cellToNode[ idx ] == INVALID_IDX )
+ // cell.cellToCell[ ( idx + 1 ) % 4 ] = this->cells.size() + newCells.size();
+
+ // if( cell.cellToNode[ ( idx - 1 ) % 4 ] == INVALID_IDX )
+ // cell.cellToCell[ ( idx - 1 ) % 4 ] = this->cells.size() + newCells.size();
+
+ // uint neighborIdx = cell.cellToCell[ ( idx + 1 ) % 4 ];
+ // if( neighborIdx != INVALID_IDX && neighborIdx < this->cells.size() )
+ // this->cells[ neighborIdx ].cellToCell[ 4 + ( idx - 1 ) % 4 ] = this->cells.size() + newCells.size();
+
+ // neighborIdx = cell.cellToCell[ ( idx - 1 ) % 4 ];
+ // if( neighborIdx != INVALID_IDX && neighborIdx < this->cells.size() )
+ // this->cells[ neighborIdx ].cellToCell[ 4 + ( idx + 0 ) % 4 ] = this->cells.size() + newCells.size();
+
+ // newCells.push_back( MeshCell() );
+ //
+ // newCells.back().isGhostCell = true;
+ // newCells.back().type = STOPPER_SOLID;
+ // newCells.back().level = cell.level;
+ // }
+ // }
+ //}
+
+ //this->cells.insert( this->cells.end(), newCells.begin(), newCells.end() );
+}
+
+void GksMeshAdapter::computeGhostCellGeometry()
+{
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( cell.type != BC_SOLID ) continue;
+
+ // for( uint idx = 0; idx < 4; idx++ ){
+ //
+ // MeshCell& neighborCell = this->cells[ cell.cellToCell[ idx ] ];
+
+ // if( neighborCell.type != STOPPER_SOLID ) continue;
+ //
+ // arraytypes::Vec2 faceCenter;
+
+ // uint node0 = cell.cellToNode[ ( idx - 1 ) % 4 ];
+ // uint node1 = cell.cellToNode[ idx ];
+
+ // if( node0 == INVALID_IDX ) continue;
+
+ // if( node1 == INVALID_IDX ) node1 = cell.cellToNode[ ( idx + 1 ) % 4 ];
+
+ // faceCenter.x = 0.5 * ( this->nodes[ node0 ].x + this->nodes[ node1 ].x );
+ // faceCenter.y = 0.5 * ( this->nodes[ node0 ].y + this->nodes[ node1 ].y );
+
+ // arraytypes::Vec2 newNode;
+
+ // //newNode.x = 2.0 * faceCenter.x - cell.cellCenter.x;
+ // //newNode.y = 2.0 * faceCenter.y - cell.cellCenter.y;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // arraytypes::Vec2 v01 ( this->nodes[ node1 ].x - this->nodes[ node0 ].x,
+ // this->nodes[ node1 ].y - this->nodes[ node0 ].y );
+
+ // arraytypes::Vec2 vD0 ( this->nodes[ node0 ].x - cell.cellCenter.x,
+ // this->nodes[ node0 ].y - cell.cellCenter.y );
+
+ // double alpha = - ( vD0.x * v01.x + vD0.y * v01.y ) / ( v01.x * v01.x + v01.y * v01.y );
+
+ // arraytypes::Vec2 vDP( this->nodes[ node0 ].x + alpha * v01.x - cell.cellCenter.x,
+ // this->nodes[ node0 ].y + alpha * v01.y - cell.cellCenter.y );
+
+ // newNode.x = cell.cellCenter.x + 2.0 * vDP.x;
+ // newNode.y = cell.cellCenter.y + 2.0 * vDP.y;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // neighborCell.cellCenter = newNode;
+ //
+ // neighborCell.cellToNode[0] = node1;
+ // neighborCell.cellToNode[1] = node0;
+ // neighborCell.cellToNode[2] = this->nodes.size();
+
+ // // dummy value to identify it as ghost cell
+ // neighborCell.cellVolume = -1.0;
+
+ // this->nodes.push_back( newNode );
+ // }
+ //}
+}
+
+void GksMeshAdapter::generateFaces(SPtr<MultipleGridBuilder> gridBuilder)
+{
+ std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ this->faces.reserve( 2 * this->cells.size() );
+
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ if( cell.type == BC_SOLID || cell.type == STOPPER_SOLID ) continue;
+
+ // generate faces in positive direction
+ for( uint neighborIdx = 0; neighborIdx < 6; neighborIdx += 2 ){
+
+ if( cell.faceExists[ neighborIdx ] ) continue;
+
+ if( cell.cellToCell[ neighborIdx ] == INVALID_INDEX ) continue;
+
+ uint neighborCellIdx = cell.cellToCell[ neighborIdx ];
+
+ MeshCell& neighborCell = this->cells[ neighborCellIdx ];
+
+ if( cell.isGhostCell && neighborCell.isGhostCell ) continue;
+
+ if( cell.isCoarseGhostCell() || neighborCell.isCoarseGhostCell() ) continue;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ MeshFace newFace;
+
+ newFace.level = cell.level;
+
+ if( neighborIdx == 0 )
+ {
+ newFace.faceToNode[ 0 ] = cell.cellToNode[ 3 ];
+ newFace.faceToNode[ 1 ] = cell.cellToNode[ 1 ];
+ newFace.faceToNode[ 2 ] = cell.cellToNode[ 0 ];
+ newFace.faceToNode[ 3 ] = cell.cellToNode[ 2 ];
+ newFace.orientation = 'x';
+ }
+ if( neighborIdx == 2 )
+ {
+ newFace.faceToNode[ 0 ] = cell.cellToNode[ 5 ];
+ newFace.faceToNode[ 1 ] = cell.cellToNode[ 4 ];
+ newFace.faceToNode[ 2 ] = cell.cellToNode[ 0 ];
+ newFace.faceToNode[ 3 ] = cell.cellToNode[ 1 ];
+ newFace.orientation = 'y';
+ }
+ if( neighborIdx == 4 )
+ {
+ newFace.faceToNode[ 0 ] = cell.cellToNode[ 6 ];
+ newFace.faceToNode[ 1 ] = cell.cellToNode[ 2 ];
+ newFace.faceToNode[ 2 ] = cell.cellToNode[ 0 ];
+ newFace.faceToNode[ 3 ] = cell.cellToNode[ 4 ];
+ newFace.orientation = 'z';
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ cell.faceExists[ neighborIdx ] = true;
+
+ // register face at neighbor
+ for( uint idx = 0; idx < 6; idx++ ){
+ if( neighborCell.cellToCell[ idx ] == cellIdx ){
+ neighborCell.faceExists[ idx ] = true;
+ break;
+ }
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ newFace.negCell = cellIdx;
+ newFace.posCell = neighborCellIdx;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ if ( cell.type == FLUID_CFF && neighborCell.type == FLUID_FCF ) newFace.negCellCoarse = cell.parent;
+ if ( cell.type == FLUID_FCF && neighborCell.type == FLUID_CFF ) newFace.posCellCoarse = neighborCell.parent;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ Vec3 faceCenter;
+
+ for( uint node = 0; node < 4; node++ ){
+ faceCenter = faceCenter + this->nodes[ newFace.faceToNode[node] ];
+ }
+
+ newFace.faceCenter.x = faceCenter.x / four;
+ newFace.faceCenter.y = faceCenter.y / four;
+ newFace.faceCenter.z = faceCenter.z / four;
+
+ this->faces.push_back( newFace );
+ }
+ }
+}
+
+void GksMeshAdapter::generateSolidFaces(SPtr<MultipleGridBuilder> gridBuilder)
+{
+
+ //std::vector< SPtr<Grid> > grids = gridBuilder->getGrids();
+
+ //this->faces.reserve( 2 * this->cells.size() );
+
+ //for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+ //
+ // MeshCell& cell = this->cells[ cellIdx ];
+
+ // if( ! ( cell.type == BC_SOLID ) ) continue;
+
+ // for( uint neighborIdx = 0; neighborIdx < 4; neighborIdx++ ){
+
+ // if( cell.faceExists[ neighborIdx ] ) continue;
+
+ // if( cell.cellToCell[ neighborIdx ] == INVALID_IDX ) continue;
+
+ // uint neighborCellIdx = cell.cellToCell[ neighborIdx ];
+
+ // MeshCell& neighborCell = this->cells[ neighborCellIdx ];
+
+ // if( cell.isGhostCell && neighborCell.isGhostCell ) continue;
+
+ // if( cell.isCoarseGhostCell() || neighborCell.isCoarseGhostCell() ) continue;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // MeshFace newFace;
+
+ // newFace.level = cell.level;
+
+ // uint node0 = cell.cellToNode[ ( neighborIdx - 1 ) % 4 ];
+ // uint node1 = cell.cellToNode[ neighborIdx ];
+
+ // if( node0 == INVALID_IDX ) continue;
+
+ // if( node1 == INVALID_IDX ) node1 = cell.cellToNode[ ( neighborIdx + 1 ) % 4 ];
+
+ // newFace.faceToNode[ 0 ] = node0;
+ // newFace.faceToNode[ 1 ] = node1;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // newFace.orientation = 'a';
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // cell.faceExists[ neighborIdx ] = true;
+
+ // // register face at neighbor
+ // for( uint idx = 0; idx < 4; idx++ ){
+ // if( neighborCell.cellToCell[ idx ] == cellIdx ){
+ // neighborCell.faceExists[ idx ] = true;
+ // break;
+ // }
+ // }
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // newFace.faceToCell[ 0 ] = cellIdx;
+ // newFace.faceToCell[ 1 ] = neighborCellIdx;
+ //
+ // //////////////////////////////////////////////////////////////////////////
+ // //////////////////////////////////////////////////////////////////////////
+
+ // if( false )
+ // {
+ // if (neighborCell.type != STOPPER_SOLID) {
+ // newFace.faceToCell[2] = cell.cellToCell[(neighborIdx + 1) % 4];
+ // newFace.faceToCell[3] = cell.cellToCell[(neighborIdx - 1) % 4];
+
+ // newFace.faceToCell[4] = cell.cellToCell[neighborIdx + 4];
+ // newFace.faceToCell[5] = cell.cellToCell[neighborIdx + 3];
+ // }
+ // else {
+ // newFace.faceToCell[2] = cell.cellToCell[(neighborIdx + 1) % 4];
+ // newFace.faceToCell[3] = cell.cellToCell[(neighborIdx - 1) % 4];
+
+ // if (newFace.faceToCell[2] == newFace.faceToCell[1])
+ // newFace.faceToCell[2] = cell.cellToCell[(neighborIdx + 2) % 4];
+
+ // if (newFace.faceToCell[3] == newFace.faceToCell[1])
+ // newFace.faceToCell[3] = cell.cellToCell[(neighborIdx - 2) % 4];
+ // }
+
+ // }
+ // //////////////////////////////////////////////////////////////////////////
+
+ // if( true )
+ // {
+ // // search for suitable cells
+
+ // std::vector<uint> candidates;
+
+ // std::function<void(uint)> findCandidates = [&](uint candidateIdx)
+ // {
+ // MeshCell candidateCell = this->cells[candidateIdx];
+
+ // bool sharesNode = false;
+
+ // for (auto faceNodeIdx : newFace.faceToNode)
+ // for (auto candidateNodeIdx : candidateCell.cellToNode)
+ // if (faceNodeIdx == candidateNodeIdx)
+ // sharesNode = true;
+
+ // if (!sharesNode) return;
+
+ // candidates.push_back(candidateIdx);
+
+ // for (auto candidateCellIdx : candidateCell.cellToCell)
+ // if (candidateCellIdx != INVALID_IDX)
+ // if (std::find(candidates.begin(), candidates.end(), candidateCellIdx) == candidates.end())
+ // findCandidates(candidateCellIdx);
+ // };
+
+ // findCandidates(cellIdx);
+
+ // // https://stackoverflow.com/questions/3385229/c-erase-vector-element-by-value-rather-than-by-position
+ // candidates.erase(std::remove(candidates.begin(), candidates.end(), cellIdx), candidates.end());
+ // candidates.erase(std::remove(candidates.begin(), candidates.end(), neighborCellIdx), candidates.end());
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // newFace.faceCenter.x = 0.5 * (this->nodes[newFace.faceToNode[1]].x + this->nodes[newFace.faceToNode[0]].x);
+ // newFace.faceCenter.y = 0.5 * (this->nodes[newFace.faceToNode[1]].y + this->nodes[newFace.faceToNode[0]].y);
+
+ // // choose cells
+
+ // if (candidates.size() <= 4) {
+
+ // uint j = 2;
+ // for (uint i = 0; i < candidates.size(); i++) {
+ // newFace.faceToCell[j++] = candidates[i];
+ // }
+ // while (j < 6) {
+ // newFace.faceToCell[j++] = INVALID_IDX;
+ // }
+ // }
+ // else
+ // {
+ // uint numberOfPermutations = pow(candidates.size(), (candidates.size() - 4));
+
+ // uint minimalPermutation = INVALID_IDX;
+ // double minimalDeviation = 1.0e99;
+
+ // for (uint permutation = 0; permutation < numberOfPermutations; permutation++)
+ // {
+ // std::vector<uint> jumpList;
+ // uint tmp = permutation;
+ // do
+ // {
+ // jumpList.push_back(tmp % candidates.size());
+ // tmp = tmp / candidates.size();
+ // } while (tmp > 0);
+
+ // std::sort(jumpList.begin(), jumpList.end());
+ // jumpList.erase(std::unique(jumpList.begin(), jumpList.end()), jumpList.end());
+
+ // if (jumpList.size() != (candidates.size() - 4)) continue;
+
+ // arraytypes::Vec2 centroid;
+
+ // centroid.x += 1.0 / 6.0 * cell.cellCenter.x;
+ // centroid.y += 1.0 / 6.0 * cell.cellCenter.y;
+ // centroid.x += 1.0 / 6.0 * neighborCell.cellCenter.x;
+ // centroid.y += 1.0 / 6.0 * neighborCell.cellCenter.y;
+
+ // for (uint i = 0; i < candidates.size(); i++)
+ // {
+ // if (std::find(jumpList.begin(), jumpList.end(), i) == jumpList.end())
+ // {
+ // centroid.x += 1.0 / 6.0 * this->cells[candidates[i]].cellCenter.x;
+ // centroid.y += 1.0 / 6.0 * this->cells[candidates[i]].cellCenter.y;
+ // }
+ // }
+
+ // double deviation = sqrt((newFace.faceCenter.x - centroid.x)
+ // * (newFace.faceCenter.x - centroid.x)
+ // + (newFace.faceCenter.y - centroid.y)
+ // * (newFace.faceCenter.y - centroid.y));
+
+ // if (deviation < minimalDeviation)
+ // {
+ // minimalPermutation = permutation;
+ // minimalDeviation = deviation;
+ // }
+ // }
+
+ // std::vector<uint> jumpList;
+ // uint tmp = minimalPermutation;
+ // do
+ // {
+ // jumpList.push_back(tmp % candidates.size());
+ // tmp = tmp / candidates.size();
+ // } while (tmp > 0);
+
+ // uint j = 2;
+ // for (uint i = 0; i < candidates.size(); i++)
+ // {
+ // if (std::find(jumpList.begin(), jumpList.end(), i) == jumpList.end())
+ // {
+ // newFace.faceToCell[j++] = candidates[i];
+ // }
+ // }
+ // }
+ // }
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // newFace.negCell = cellIdx;
+ // newFace.posCell = neighborCellIdx;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // if ( cell.type == FLUID_CFF && neighborCell.type == FLUID_FCF ) newFace.negCell = cell.parent;
+ // if ( cell.type == FLUID_FCF && neighborCell.type == FLUID_CFF ) newFace.posCell = neighborCell.parent;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // arraytypes::Vec2 faceVec( this->nodes[ newFace.faceToNode[ 1 ] ].x - this->nodes[ newFace.faceToNode[ 0 ] ].x,
+ // this->nodes[ newFace.faceToNode[ 1 ] ].y - this->nodes[ newFace.faceToNode[ 0 ] ].y );
+
+ // newFace.faceArea = sqrt( faceVec.x * faceVec.x + faceVec.y * faceVec.y );
+
+ // newFace.faceCenter.x = 0.5 * ( this->nodes[ newFace.faceToNode[ 1 ] ].x + this->nodes[ newFace.faceToNode[ 0 ] ].x );
+ // newFace.faceCenter.y = 0.5 * ( this->nodes[ newFace.faceToNode[ 1 ] ].y + this->nodes[ newFace.faceToNode[ 0 ] ].y );
+
+ // newFace.faceNormal.x = faceVec.y / newFace.faceArea;
+ // newFace.faceNormal.y = - faceVec.x / newFace.faceArea;
+
+ // //////////////////////////////////////////////////////////////////////////
+
+ // this->faces.push_back( newFace );
+ // }
+ //}
+}
+
+void GksMeshAdapter::sortFaces()
+{
+ real xMax = ( *std::max_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.x < rhs.faceCenter.x; } ) ).faceCenter.x;
+ real yMax = ( *std::max_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.y < rhs.faceCenter.y; } ) ).faceCenter.y;
+ real zMax = ( *std::max_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.z < rhs.faceCenter.z; } ) ).faceCenter.z;
+
+ real xMin = ( *std::min_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.x < rhs.faceCenter.x; } ) ).faceCenter.x;
+ real yMin = ( *std::min_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.y < rhs.faceCenter.y; } ) ).faceCenter.y;
+ real zMin = ( *std::min_element( this->faces.begin(), this->faces.end(), [this]( MeshFace lhs, MeshFace rhs ){ return lhs.faceCenter.z < rhs.faceCenter.z; } ) ).faceCenter.z;
+
+ real xRange = xMax - xMin;
+ real yRange = yMax - yMin;
+ real zRange = zMax - zMin;
+
+ std::sort( this->faces.begin(), this->faces.end(),
+ [&,this]( MeshFace lhs, MeshFace rhs ){
+
+ if( lhs.level != rhs.level ) return lhs.level < rhs.level;
+
+ if( lhs.orientation != rhs.orientation ){
+ if ( lhs.orientation == 'x' && rhs.orientation == 'y' ) return true;
+ else if ( lhs.orientation == 'y' && rhs.orientation == 'z' ) return true;
+ else if ( lhs.orientation == 'x' && rhs.orientation == 'z' ) return true;
+ else return false;
+ }
+
+ return ( ten * ten * lhs.faceCenter.z / zRange + ten * lhs.faceCenter.y / yRange + lhs.faceCenter.x / xRange ) <
+ ( ten * ten * rhs.faceCenter.z / zRange + ten * rhs.faceCenter.y / yRange + rhs.faceCenter.x / xRange );
+ }
+ );
+}
+
+void GksMeshAdapter::countFaces()
+{
+ this->numberOfFacesPerLevelXYZ.resize( 3 * this->numberOfLevels );
+ this->startOfFacesPerLevelXYZ.resize ( 3 * this->numberOfLevels );
+
+ for( auto& i : this->numberOfFacesPerLevelXYZ ) i = 0;
+ for( auto& i : this->startOfFacesPerLevelXYZ ) i = 0;
+
+ for( auto& face : this->faces ){
+ if ( face.orientation == 'x' ) this->numberOfFacesPerLevelXYZ[ 3 * face.level ]++;
+ else if ( face.orientation == 'y' ) this->numberOfFacesPerLevelXYZ[ 3 * face.level + 1 ]++;
+ else if ( face.orientation == 'z' ) this->numberOfFacesPerLevelXYZ[ 3 * face.level + 2 ]++;
+ }
+
+ this->startOfFacesPerLevelXYZ[0] = 0;
+
+ for( uint level = 1; level < 3 * this->numberOfLevels; level++ ){
+
+ this->startOfFacesPerLevelXYZ[level] = this->startOfFacesPerLevelXYZ [level - 1]
+ + this->numberOfFacesPerLevelXYZ[level - 1];
+ }
+}
+
+void GksMeshAdapter::findSolidFaces()
+{
+ //for( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ // MeshFace& face = this->faces[ faceIdx ];
+
+ // if( this->cells[ face.posCell ].type == STOPPER_SOLID )
+ // this->solidFaces.push_back( faceIdx );
+ //}
+}
+
+void GksMeshAdapter::generateInterfaceConnectivity()
+{
+ this->numberOfFineToCoarsePerLevel.resize( this->numberOfLevels );
+ this->startOfFineToCoarsePerLevel.resize ( this->numberOfLevels );
+ this->numberOfCoarseToFinePerLevel.resize( this->numberOfLevels );
+ this->startOfCoarseToFinePerLevel.resize ( this->numberOfLevels );
+
+ for( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ if( cell.type == FLUID_FCC ){
+
+ uint_9 connectivity;
+
+ connectivity[ 0 ] = cellIdx;
+ connectivity[ 1 ] = cell.children[ 0 ];
+ connectivity[ 2 ] = cell.children[ 1 ];
+ connectivity[ 3 ] = cell.children[ 2 ];
+ connectivity[ 4 ] = cell.children[ 3 ];
+ connectivity[ 5 ] = cell.children[ 4 ];
+ connectivity[ 6 ] = cell.children[ 5 ];
+ connectivity[ 7 ] = cell.children[ 6 ];
+ connectivity[ 8 ] = cell.children[ 7 ];
+
+ this->fineToCoarse.push_back( connectivity );
+
+ this->numberOfFineToCoarsePerLevel[ cell.level ]++;
+ }
+
+ if( cell.type == FLUID_CFC ){
+
+ uint_15 connectivity;
+
+ connectivity[ 0 ] = cellIdx;
+
+ connectivity[ 1 ] = cell.cellToCell[ 0 ];
+ connectivity[ 2 ] = cell.cellToCell[ 1 ];
+ connectivity[ 3 ] = cell.cellToCell[ 2 ];
+ connectivity[ 4 ] = cell.cellToCell[ 3 ];
+ connectivity[ 5 ] = cell.cellToCell[ 4 ];
+ connectivity[ 6 ] = cell.cellToCell[ 5 ];
+
+ connectivity[ 7 ] = cell.children[ 0 ];
+ connectivity[ 8 ] = cell.children[ 1 ];
+ connectivity[ 9 ] = cell.children[ 2 ];
+ connectivity[ 10 ] = cell.children[ 3 ];
+ connectivity[ 11 ] = cell.children[ 4 ];
+ connectivity[ 12 ] = cell.children[ 5 ];
+ connectivity[ 13 ] = cell.children[ 6 ];
+ connectivity[ 14 ] = cell.children[ 7 ];
+
+ this->coarseToFine.push_back( connectivity );
+
+ numberOfCoarseToFinePerLevel[ cell.level ]++;
+ }
+ }
+
+ std::cout << numberOfCoarseToFinePerLevel[ 0 ] << " " << this->numberOfFineToCoarsePerLevel[ 0 ] << std::endl;
+
+ this->startOfFineToCoarsePerLevel[0] = 0;
+ this->startOfCoarseToFinePerLevel[0] = 0;
+
+ for( uint level = 1; level < this->numberOfLevels; level++ ){
+
+ this->startOfFineToCoarsePerLevel[level] = this->startOfFineToCoarsePerLevel [level - 1]
+ + this->numberOfFineToCoarsePerLevel[level - 1];
+
+ this->startOfCoarseToFinePerLevel[level] = this->startOfCoarseToFinePerLevel [level - 1]
+ + this->numberOfCoarseToFinePerLevel[level - 1];
+ }
+}
+
+void GksMeshAdapter::writeMeshVTK(std::string filename)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "writeMeshVTK( " << filename << " )" << "\n";
+
+ std::ofstream file;
+
+ file.open(filename);
+
+ file << "# vtk DataFile Version 3.0\n";
+ file << "by MeshGenerator\n";
+ file << "ASCII\n";
+ file << "DATASET UNSTRUCTURED_GRID\n";
+
+ file << "POINTS " << nodes.size() << " float" << std::endl;
+
+ for (auto node : nodes){
+ file << node.x << " " << node.y << " " << node.z << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "CELLS " << this->cells.size() << " " << this->cells.size() * 9 << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ uint_8 nodes;
+ for( auto& i : nodes ) i = INVALID_INDEX;
+
+ nodes[0] = cell.cellToNode[7];//[ 6 ];
+ nodes[1] = cell.cellToNode[3];//[ 5 ];
+ nodes[2] = cell.cellToNode[1];//[ 2 ];
+ nodes[3] = cell.cellToNode[5];//[ 1 ];
+ nodes[4] = cell.cellToNode[6];//[ 4 ];
+ nodes[5] = cell.cellToNode[2];//[ 7 ];
+ nodes[6] = cell.cellToNode[0];//[ 0 ];
+ nodes[7] = cell.cellToNode[4];//[ 3 ];
+
+ file << 8 << " ";
+
+ for( uint i = 0; i < 8; i++ ) file << nodes[i] << " ";
+
+ file << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "CELL_TYPES " << this->cells.size() << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+ file << 12 << std::endl;
+ }
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "\nCELL_DATA " << this->cells.size() << std::endl;
+
+ file << "FIELD Label " << 4 << std::endl;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "CellIdx 1 " << this->cells.size() << " int" << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ file << cellIdx << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "level 1 " << this->cells.size() << " int" << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ file << cell.level << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "type 1 " << this->cells.size() << " int" << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ file << (uint) cell.type << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "isGhostCell 1 " << this->cells.size() << " int" << std::endl;
+
+ for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ MeshCell& cell = this->cells[ cellIdx ];
+
+ file << (uint) cell.isGhostCell << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file.close();
+}
+
+void GksMeshAdapter::writeMeshFaceVTK(std::string filename)
+{
+ *logging::out << logging::Logger::INFO_INTERMEDIATE << "writeMeshFaceVTK( " << filename << " )" << "\n";
+
+ std::ofstream file;
+
+ file.open(filename);
+
+ file << "# vtk DataFile Version 3.0\n";
+ file << "by MeshGenerator\n";
+ file << "ASCII\n";
+ file << "DATASET UNSTRUCTURED_GRID\n";
+
+ file << "POINTS " << nodes.size() << " float" << std::endl;
+
+ for (auto node : nodes){
+ file << node.x << " " << node.y << " " << node.z << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "CELLS " << this->faces.size() << " " << 5 * this->faces.size() << std::endl;
+
+ for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ file << "4 ";
+
+ file << this->faces[ faceIdx ].faceToNode[0] << " ";
+ file << this->faces[ faceIdx ].faceToNode[1] << " ";
+ file << this->faces[ faceIdx ].faceToNode[2] << " ";
+ file << this->faces[ faceIdx ].faceToNode[3] << " ";
+
+ file << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "CELL_TYPES " << this->faces.size() << std::endl;
+
+ for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+ file << "9" << std::endl;
+ }
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "\nCELL_DATA " << this->faces.size() << std::endl;
+
+ file << "FIELD Label " << 3 << std::endl;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "FaceIdx 1 " << this->faces.size() << " int" << std::endl;
+
+ for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ file << faceIdx << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "level 1 " << this->faces.size() << " int" << std::endl;
+
+ for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ file << this->faces[ faceIdx ].level << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "orientation 1 " << this->faces.size() << " int" << std::endl;
+
+ for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ file << (int)this->faces[ faceIdx ].orientation << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+
+ file << "VECTORS posCell double" << std::endl;
+
+ for ( auto face : this->faces )
+ {
+ MeshCell& cell = this->cells[ face.posCell ];
+
+ Vec3 vec = cell.cellCenter - face.faceCenter;
+
+ file << vec.x << " ";
+ file << vec.y << " ";
+ file << vec.z << std::endl;
+ }
+
+ file << "VECTORS negCell double" << std::endl;
+
+ for ( auto face : this->faces )
+ {
+ MeshCell& cell = this->cells[ face.negCell ];
+
+ Vec3 vec = cell.cellCenter - face.faceCenter;
+
+ file << vec.x << " ";
+ file << vec.y << " ";
+ file << vec.z << std::endl;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+
+ file.close();
+}
+
+void GksMeshAdapter::writeMeshCellToCellVTK(std::string filename)
+{
+ //std::ofstream file;
+
+ //file.open(filename);
+
+ //file << "# vtk DataFile Version 3.0\n";
+ //file << "by MeshGenerator\n";
+ //file << "ASCII\n";
+ //file << "DATASET UNSTRUCTURED_GRID\n";
+
+ //file << "POINTS " << this->cells.size() << " float" << std::endl;
+
+ //for (auto cell : cells){
+ // file << cell.cellCenter.x << " " << cell.cellCenter.y << " 0.0" << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CELLS " << 8 * this->cells.size() << " " << 3 * 8 * this->cells.size() << std::endl;
+
+ //for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // for( uint i = 0; i < 8; i++ )
+ // if( this->cells[ cellIdx ].cellToCell[ i ] != INVALID_INDEX )
+ // file << "2 " << cellIdx << " " << this->cells[ cellIdx ].cellToCell[ i ] << " " << std::endl;
+ // else
+ // file << "2 " << cellIdx << " " << cellIdx << " " << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CELL_TYPES " << 8 * this->cells.size() << std::endl;
+
+ //for ( uint i = 0; i < 8 * this->cells.size(); i++ ){
+ // file << "3" << std::endl;
+ //}
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "\nCELL_DATA " << 8 * this->cells.size() << std::endl;
+
+ //file << "FIELD Label " << 2 << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CellIdx 1 " << 8 * this->cells.size() << " int" << std::endl;
+
+ //for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // for( uint i = 0; i < 8; i++ )
+ // file << cellIdx << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CellToCell 1 " << 8 * this->cells.size() << " int" << std::endl;
+
+ //for ( uint cellIdx = 0; cellIdx < this->cells.size(); cellIdx++ ){
+
+ // for( uint i = 0; i < 8; i++ )
+ // file << i << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file.close();
+}
+
+void GksMeshAdapter::writeMeshFaceToCellVTK(std::string filename)
+{
+ //std::ofstream file;
+
+ //file.open(filename);
+
+ //file << "# vtk DataFile Version 3.0\n";
+ //file << "by MeshGenerator\n";
+ //file << "ASCII\n";
+ //file << "DATASET UNSTRUCTURED_GRID\n";
+
+ //file << "POINTS " << this->cells.size() + this->faces.size() << " float" << std::endl;
+
+ //for (auto cell : cells){
+ // file << cell.cellCenter.x << " " << cell.cellCenter.y << " 0.0" << std::endl;
+ //}
+
+ //for (auto face : faces){
+ // file << face.faceCenter.x << " " << face.faceCenter.y << " 0.0" << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CELLS " << 6 * this->faces.size() << " " << 3 * 6 * this->faces.size() << std::endl;
+
+ //for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ // for( uint i = 0; i < 6; i++ )
+ // if( this->faces[ faceIdx ].faceToCell[ i ] != INVALID_INDEX )
+ // file << "2 " << this->cells.size() + faceIdx << " " << this->faces[ faceIdx ].faceToCell[ i ] << " " << std::endl;
+ // else
+ // file << "2 " << this->cells.size() + faceIdx << " " << this->cells.size() + faceIdx << " " << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "CELL_TYPES " << 6 * this->faces.size() << std::endl;
+
+ //for ( uint i = 0; i < 6 * this->faces.size(); i++ ){
+ // file << "3" << std::endl;
+ //}
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "\nCELL_DATA " << 6 * this->faces.size() << std::endl;
+
+ //file << "FIELD Label " << 2 << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "FaceIdx 1 " << 6 * this->faces.size() << " int" << std::endl;
+
+ //for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ // for( uint i = 0; i < 6; i++ )
+ // file << faceIdx << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file << "FaceToCell 1 " << 6 * this->faces.size() << " int" << std::endl;
+
+ //for ( uint faceIdx = 0; faceIdx < this->faces.size(); faceIdx++ ){
+
+ // for( uint i = 0; i < 6; i++ )
+ // file << i << std::endl;
+ //}
+
+ ////////////////////////////////////////////////////////////////////////////
+
+ //file.close();
+}
+
+double GksMeshAdapter::getDx(uint level)
+{
+ return dxCoarse / pow( 2.0, level );
+}
diff --git a/src/GksMeshAdapter/GksMeshAdapter.h b/src/GksMeshAdapter/GksMeshAdapter.h
new file mode 100644
index 000000000..ce53d4946
--- /dev/null
+++ b/src/GksMeshAdapter/GksMeshAdapter.h
@@ -0,0 +1,151 @@
+#ifndef GKS_MESH_ADAPTER_H
+#define GKS_MESH_ADAPTER_H
+
+#include <memory>
+#include <array>
+#include <vector>
+
+#include "Core/DataTypes.h"
+#include "Core/PointerDefinitions.h"
+
+#include "MeshCell.h"
+#include "MeshFace.h"
+
+
+class MultipleGridBuilder;
+
+class VF_PUBLIC GksMeshAdapter{
+
+public:
+
+ uint numberOfLevels;
+
+ double dxCoarse;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ std::vector<Vec3> nodes;
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // C e l l s o r t i n g :
+ //
+ // | Level 0 | Level 1 | Level 2 |
+ // | FluidCells | GhostCells | FluidCells | GhostCells | FluidCells | GhostCells |
+ //
+ // GhostCells: not included in Cell update, i.e. BoundaryCells and FCC-Cells
+ // FluidCells: all other, including CFF-Cells
+ //
+ std::vector<MeshCell> cells;
+
+ std::vector<uint> numberOfCellsPerLevel;
+ std::vector<uint> numberOfBulkCellsPerLevel;
+ std::vector<uint> startOfCellsPerLevel;
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // F a c e s o r t i n g :
+ //
+ // | Level 0 | Level 1 | Level 2 |
+ // | x-normal | y-normal | z-normal | x-normal | y-normal | z-normal | x-normal | y-normal | z-normal |
+ //
+ std::vector<MeshFace> faces;
+
+ std::vector<uint> numberOfFacesPerLevelXYZ;
+ std::vector<uint> startOfFacesPerLevelXYZ;
+
+ //std::vector<uint> solidFaces;
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // F i n e T o C o a r s e s o r t i n g :
+ //
+ // | Coarse Cell Idx | Child Idcs | ...
+ //
+ std::vector<uint_9> fineToCoarse;
+
+ std::vector<uint> numberOfFineToCoarsePerLevel;
+ std::vector<uint> startOfFineToCoarsePerLevel;
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // F i n e T o C o a r s e s o r t i n g :
+ //
+ // | Coarse Cell Idx | Coarse Neighbor Idcs | Child Idcs | ...
+ //
+ std::vector<uint_15> coarseToFine;
+
+ std::vector<uint> numberOfCoarseToFinePerLevel;
+ std::vector<uint> startOfCoarseToFinePerLevel;
+
+ //////////////////////////////////////////////////////////////////////////
+ //
+ // Connectivity from LBM grid to GKS Mesh
+ //
+ // cellIdx = gridToMesh[ level ][ gridIdx ]
+ //
+ std::vector< std::vector<uint> > gridToMesh;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ uint_8 cornerCells;
+
+ //////////////////////////////////////////////////////////////////////////
+
+public:
+
+ void inputGrid( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void findQuadtreeConnectivity( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void findCellToCellConnectivity( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void countCells();
+
+ void partitionCells();
+
+ void refreshCellConnectivity(const std::vector<uint>& idxMap);
+
+ void findCornerCells( SPtr<MultipleGridBuilder> gridBuilder, real z0 );
+
+ void generateNodes( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void inputQs( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void morphNodes();
+
+ void repairTriangularCells();
+
+ void computeCellGeometry();
+
+ void generateSolidGhostCells();
+
+ void computeGhostCellGeometry();
+
+ void generateFaces( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void generateSolidFaces( SPtr<MultipleGridBuilder> gridBuilder );
+
+ void sortFaces();
+
+ void countFaces();
+
+ void findSolidFaces();
+
+ void generateInterfaceConnectivity();
+
+ void writeMeshVTK( std::string filename );
+
+ void writeMeshFaceVTK( std::string filename );
+
+ void writeMeshCellToCellVTK( std::string filename );
+
+ void writeMeshFaceToCellVTK( std::string filename );
+
+ //////////////////////////////////////////////////////////////////////////
+
+ double getDx(uint level);
+};
+
+
+#endif
\ No newline at end of file
diff --git a/src/GksMeshAdapter/MeshCell.cpp b/src/GksMeshAdapter/MeshCell.cpp
new file mode 100644
index 000000000..b63c87053
--- /dev/null
+++ b/src/GksMeshAdapter/MeshCell.cpp
@@ -0,0 +1,32 @@
+#include "MeshCell.h"
+
+#include "GridGenerator/grid/NodeValues.h"
+
+MeshCell::MeshCell(){
+
+ level = INVALID_INDEX;
+ gridIdx = INVALID_INDEX;
+
+
+ for( uint& index : this->cellToNode ) index = INVALID_INDEX;
+ for( uint& index : this->cellToEdgeNode ) index = INVALID_INDEX;
+ for( uint& index : this->cellToFaceNode ) index = INVALID_INDEX;
+ for( uint& index : this->cellToCell ) index = INVALID_INDEX;
+ for( uint& index : this->children ) index = INVALID_INDEX;
+
+ parent = INVALID_INDEX;
+
+ for( bool& flag : this->faceExists ) flag = false;
+
+ isGhostCell = false;
+}
+
+bool MeshCell::isCoarseGhostCell()
+{
+ return this->type == FLUID_FCC;
+}
+
+bool MeshCell::isFineGhostCell()
+{
+ return this->type == FLUID_CFF;
+}
diff --git a/src/GksMeshAdapter/MeshCell.h b/src/GksMeshAdapter/MeshCell.h
new file mode 100644
index 000000000..a0e6b0ec9
--- /dev/null
+++ b/src/GksMeshAdapter/MeshCell.h
@@ -0,0 +1,71 @@
+#ifndef MESH_CELL_H
+#define MESH_CELL_H
+
+#include <array>
+
+#include "Core/DataTypes.h"
+#include "Core/VectorTypes.h"
+#include "Core/ArrayTypes.h"
+
+struct VF_PUBLIC MeshCell{
+
+ uint level;
+ uint gridIdx;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ uint_8 cellToNode; uint_12 cellToEdgeNode; uint_6 cellToFaceNode;
+
+ // for sorting see LBM f numbering
+ // 8
+ // 4 o--------o 0 o--------o o--------o
+ // /| /| 7/| 4/| /| 4 /|
+ // / | / | /3| 11 / | / | 2 / |
+ // 6 o--------o 2| o--------o |0 o--------o |
+ // | | | | | | 10 | | |1 | | 0|
+ // |5 o-----|--o 1 1| o-----|--o | o-----|--o
+ // | / | / |5/ 2| /6 | / 3 | /
+ // |/ |/ |/ 9 |/ |/ 5 |/
+ // 7 o--------o 3 o--------o o--------o
+ //
+ // z | / y
+ // |/
+ // +---- x
+ //
+
+ //////////////////////////////////////////////////////////////////////////
+
+ uint_27 cellToCell;
+
+ // for sorting see LBM f numbering
+
+ //////////////////////////////////////////////////////////////////////////
+
+ uint_8 children;
+
+ // for sorting see cellToNode
+
+ uint parent;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ Vec3 cellCenter;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ // order: +x, -x, +y, -y, +z, -z (see cellToCell)
+ bool_6 faceExists;
+
+ bool isGhostCell;
+
+ char type;
+
+ MeshCell();
+
+ bool isCoarseGhostCell();
+
+ bool isFineGhostCell();
+};
+
+
+#endif
\ No newline at end of file
diff --git a/src/GksMeshAdapter/MeshFace.cpp b/src/GksMeshAdapter/MeshFace.cpp
new file mode 100644
index 000000000..6ae2768a8
--- /dev/null
+++ b/src/GksMeshAdapter/MeshFace.cpp
@@ -0,0 +1,12 @@
+#include "MeshFace.h"
+
+MeshFace::MeshFace()
+{
+ for( uint& node : this->faceToNode ) node = INVALID_INDEX;
+
+ isWall = false;
+ posCell = INVALID_INDEX;
+ negCell = INVALID_INDEX;
+ posCellCoarse = INVALID_INDEX;
+ negCellCoarse = INVALID_INDEX;
+}
\ No newline at end of file
diff --git a/src/GksMeshAdapter/MeshFace.h b/src/GksMeshAdapter/MeshFace.h
new file mode 100644
index 000000000..a07f7f979
--- /dev/null
+++ b/src/GksMeshAdapter/MeshFace.h
@@ -0,0 +1,50 @@
+#ifndef MESH_FACE_H
+#define MESH_FACE_H
+
+#include "Core/DataTypes.h"
+#include "Core/VectorTypes.h"
+#include "Core/ArrayTypes.h"
+
+struct VF_PUBLIC MeshFace
+{
+ //////////////////////////////////////////////////////////////////////////
+
+ // o 2
+ // /|
+ // / |
+ // o 3| n
+ // | -+--------->
+ // | o 1
+ // | /
+ // |/
+ // o 0
+ //
+ //
+
+ uint_4 faceToNode;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ uint posCell;
+ uint negCell;
+
+ uint posCellCoarse;
+ uint negCellCoarse;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ Vec3 faceCenter;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ char orientation;
+
+ bool isWall;
+
+ uint level;
+
+ MeshFace();
+};
+
+
+#endif
\ No newline at end of file
diff --git a/src/GksMeshAdapter/package.include b/src/GksMeshAdapter/package.include
new file mode 100644
index 000000000..e69de29bb
diff --git a/targets/apps/GKS/gksTest/CMakeLists.txt b/targets/apps/GKS/gksTest/CMakeLists.txt
index 836b9069e..bf204134e 100644
--- a/targets/apps/GKS/gksTest/CMakeLists.txt
+++ b/targets/apps/GKS/gksTest/CMakeLists.txt
@@ -1,8 +1,10 @@
setTargetNameToFolderName(${CMAKE_CURRENT_LIST_DIR})
set(linkDirectories "")
-set(libsToLink VirtualFluids_GPU GridGenerator)
-set(includeDirectories "${CMAKE_SOURCE_DIR}/src" "${CMAKE_SOURCE_DIR}/src/VirtualFluids_GPU" "${CMAKE_SOURCE_DIR}/src/GridGenerator" "${CMAKE_SOURCE_DIR}/src/VirtualFluidsBasics")
+set(libsToLink GridGenerator GksMeshAdapter)
+set(includeDirectories "${CMAKE_SOURCE_DIR}/src"
+ "${CMAKE_SOURCE_DIR}/src/GridGenerator"
+ "${CMAKE_SOURCE_DIR}/src/GksMeshAdapter")
#glob files and save in MY_SRCS
include(CMakePackage.cmake)
diff --git a/targets/apps/GKS/gksTest/main.cpp b/targets/apps/GKS/gksTest/main.cpp
index c36fe4fa3..7907c851f 100644
--- a/targets/apps/GKS/gksTest/main.cpp
+++ b/targets/apps/GKS/gksTest/main.cpp
@@ -1,1159 +1,104 @@
//#define MPI_LOGGING
-// Stephas Branch
-
-#include <mpi.h>
-#if defined( MPI_LOGGING )
- #include <mpe.h>
-#endif
-
#include <string>
#include <iostream>
-#include <stdexcept>
+#include <exception>
#include <fstream>
#define _USE_MATH_DEFINES
#include <math.h>
-#include "metis.h"
-
-#include "Core/LbmOrGks.h"
-#include "Core/Input/Input.h"
-#include "Core/StringUtilities/StringUtil.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 "core/Logger/Logger.h"
-#include "global.h"
+#include "GridGenerator/geometries/Cuboid/Cuboid.h"
-#include "geometries/Sphere/Sphere.h"
-#include "geometries/VerticalCylinder/VerticalCylinder.h"
-#include "geometries/Cuboid/Cuboid.h"
-#include "geometries/TriangularMesh/TriangularMesh.h"
-#include "geometries/Conglomerate/Conglomerate.h"
-#include "geometries/TriangularMesh/TriangularMeshStrategy.h"
+#include "GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
+#include "GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
+#include "GridGenerator/grid/BoundaryConditions/Side.h"
+#include "GridGenerator/grid/BoundaryConditions/BoundaryCondition.h"
+#include "GridGenerator/grid/GridFactory.h"
-#include "grid/GridBuilder/LevelGridBuilder.h"
-#include "grid/GridBuilder/MultipleGridBuilder.h"
-#include "grid/BoundaryConditions/Side.h"
-#include "grid/BoundaryConditions/BoundaryCondition.h"
-#include "grid/GridFactory.h"
-
-#include "io/SimulationFileWriter/SimulationFileWriter.h"
-#include "io/GridVTKWriter/GridVTKWriter.h"
-#include "io/STLReaderWriter/STLReader.h"
-#include "io/STLReaderWriter/STLWriter.h"
-
-#include "utilities/math/Math.h"
-#include "utilities/communication.h"
-#include "utilities/transformator/TransformatorImp.h"
-
-std::string getGridPath(std::shared_ptr<Parameter> para, std::string Gridpath)
-{
- if (para->getNumprocs() == 1)
- return Gridpath + "/";
-
- return Gridpath + "/" + StringUtil::toString(para->getMyID()) + "/";
-}
+#include "GksMeshAdapter/GksMeshAdapter.h"
-void setParameters(std::shared_ptr<Parameter> para, std::unique_ptr<input::Input> &input)
+void gksTest()
{
- Communicator* comm = Communicator::getInstanz();
-
- para->setMaxDev(StringUtil::toInt(input->getValue("NumberOfDevices")));
- para->setNumprocs(comm->getNummberOfProcess());
- para->setDevices(StringUtil::toVector(input->getValue("Devices")));
- para->setMyID(comm->getPID());
-
- std::string _path = input->getValue("Path");
- std::string _prefix = input->getValue("Prefix");
- std::string _gridpath = input->getValue("GridPath");
- std::string gridPath = getGridPath(para, _gridpath);
- para->setOutputPath(_path);
- para->setOutputPrefix(_prefix);
- para->setFName(_path + "/" + _prefix);
- para->setPrintFiles(false);
- para->setPrintFiles(StringUtil::toBool(input->getValue("WriteGrid")));
- para->setGeometryValues(StringUtil::toBool(input->getValue("GeometryValues")));
- para->setCalc2ndOrderMoments(StringUtil::toBool(input->getValue("calc2ndOrderMoments")));
- para->setCalc3rdOrderMoments(StringUtil::toBool(input->getValue("calc3rdOrderMoments")));
- para->setCalcHighOrderMoments(StringUtil::toBool(input->getValue("calcHigherOrderMoments")));
- para->setReadGeo(StringUtil::toBool(input->getValue("ReadGeometry")));
- para->setCalcMedian(StringUtil::toBool(input->getValue("calcMedian")));
- para->setConcFile(StringUtil::toBool(input->getValue("UseConcFile")));
- para->setUseMeasurePoints(StringUtil::toBool(input->getValue("UseMeasurePoints")));
- para->setUseWale(StringUtil::toBool(input->getValue("UseWale")));
- para->setSimulatePorousMedia(StringUtil::toBool(input->getValue("SimulatePorousMedia")));
- para->setD3Qxx(StringUtil::toInt(input->getValue("D3Qxx")));
- para->setTEnd(StringUtil::toInt(input->getValue("TimeEnd")));
- para->setTOut(StringUtil::toInt(input->getValue("TimeOut")));
- para->setTStartOut(StringUtil::toInt(input->getValue("TimeStartOut")));
- para->setTimeCalcMedStart(StringUtil::toInt(input->getValue("TimeStartCalcMedian")));
- para->setTimeCalcMedEnd(StringUtil::toInt(input->getValue("TimeEndCalcMedian")));
- para->setPressInID(StringUtil::toInt(input->getValue("PressInID")));
- para->setPressOutID(StringUtil::toInt(input->getValue("PressOutID")));
- para->setPressInZ(StringUtil::toInt(input->getValue("PressInZ")));
- para->setPressOutZ(StringUtil::toInt(input->getValue("PressOutZ")));
- //////////////////////////////////////////////////////////////////////////
- para->setDiffOn(StringUtil::toBool(input->getValue("DiffOn")));
- para->setDiffMod(StringUtil::toInt(input->getValue("DiffMod")));
- para->setDiffusivity(StringUtil::toFloat(input->getValue("Diffusivity")));
- para->setTemperatureInit(StringUtil::toFloat(input->getValue("Temp")));
- para->setTemperatureBC(StringUtil::toFloat(input->getValue("TempBC")));
- //////////////////////////////////////////////////////////////////////////
- para->setViscosity(StringUtil::toFloat(input->getValue("Viscosity_LB")));
- para->setVelocity(StringUtil::toFloat(input->getValue("Velocity_LB")));
- para->setViscosityRatio(StringUtil::toFloat(input->getValue("Viscosity_Ratio_World_to_LB")));
- para->setVelocityRatio(StringUtil::toFloat(input->getValue("Velocity_Ratio_World_to_LB")));
- para->setDensityRatio(StringUtil::toFloat(input->getValue("Density_Ratio_World_to_LB")));
- para->setPressRatio(StringUtil::toFloat(input->getValue("Delta_Press")));
- para->setRealX(StringUtil::toFloat(input->getValue("SliceRealX")));
- para->setRealY(StringUtil::toFloat(input->getValue("SliceRealY")));
- para->setFactorPressBC(StringUtil::toFloat(input->getValue("dfpbc")));
- para->setGeometryFileC(input->getValue("GeometryC"));
- para->setGeometryFileM(input->getValue("GeometryM"));
- para->setGeometryFileF(input->getValue("GeometryF"));
- //////////////////////////////////////////////////////////////////////////
- para->setgeoVec(gridPath + input->getValue("geoVec"));
- para->setcoordX(gridPath + input->getValue("coordX"));
- para->setcoordY(gridPath + input->getValue("coordY"));
- para->setcoordZ(gridPath + input->getValue("coordZ"));
- para->setneighborX(gridPath + input->getValue("neighborX"));
- para->setneighborY(gridPath + input->getValue("neighborY"));
- para->setneighborZ(gridPath + input->getValue("neighborZ"));
- para->setscaleCFC(gridPath + input->getValue("scaleCFC"));
- para->setscaleCFF(gridPath + input->getValue("scaleCFF"));
- para->setscaleFCC(gridPath + input->getValue("scaleFCC"));
- para->setscaleFCF(gridPath + input->getValue("scaleFCF"));
- para->setscaleOffsetCF(gridPath + input->getValue("scaleOffsetCF"));
- para->setscaleOffsetFC(gridPath + input->getValue("scaleOffsetFC"));
- para->setgeomBoundaryBcQs(gridPath + input->getValue("geomBoundaryBcQs"));
- para->setgeomBoundaryBcValues(gridPath + input->getValue("geomBoundaryBcValues"));
- para->setinletBcQs(gridPath + input->getValue("inletBcQs"));
- para->setinletBcValues(gridPath + input->getValue("inletBcValues"));
- para->setoutletBcQs(gridPath + input->getValue("outletBcQs"));
- para->setoutletBcValues(gridPath + input->getValue("outletBcValues"));
- para->settopBcQs(gridPath + input->getValue("topBcQs"));
- para->settopBcValues(gridPath + input->getValue("topBcValues"));
- para->setbottomBcQs(gridPath + input->getValue("bottomBcQs"));
- para->setbottomBcValues(gridPath + input->getValue("bottomBcValues"));
- para->setfrontBcQs(gridPath + input->getValue("frontBcQs"));
- para->setfrontBcValues(gridPath + input->getValue("frontBcValues"));
- para->setbackBcQs(gridPath + input->getValue("backBcQs"));
- para->setbackBcValues(gridPath + input->getValue("backBcValues"));
- para->setnumberNodes(gridPath + input->getValue("numberNodes"));
- para->setLBMvsSI(gridPath + input->getValue("LBMvsSI"));
- //////////////////////////////gridPath + ////////////////////////////////////////////
- para->setmeasurePoints(gridPath + input->getValue("measurePoints"));
- para->setpropellerValues(gridPath + input->getValue("propellerValues"));
- para->setclockCycleForMP(StringUtil::toFloat(input->getValue("measureClockCycle")));
- para->settimestepForMP(StringUtil::toInt(input->getValue("measureTimestep")));
- para->setcpTop(gridPath + input->getValue("cpTop"));
- para->setcpBottom(gridPath + input->getValue("cpBottom"));
- para->setcpBottom2(gridPath + input->getValue("cpBottom2"));
- para->setConcentration(gridPath + input->getValue("Concentration"));
- //////////////////////////////////////////////////////////////////////////
- //Normals - Geometry
- para->setgeomBoundaryNormalX(gridPath + input->getValue("geomBoundaryNormalX"));
- para->setgeomBoundaryNormalY(gridPath + input->getValue("geomBoundaryNormalY"));
- para->setgeomBoundaryNormalZ(gridPath + input->getValue("geomBoundaryNormalZ"));
- //Normals - Inlet
- para->setInflowBoundaryNormalX(gridPath + input->getValue("inletBoundaryNormalX"));
- para->setInflowBoundaryNormalY(gridPath + input->getValue("inletBoundaryNormalY"));
- para->setInflowBoundaryNormalZ(gridPath + input->getValue("inletBoundaryNormalZ"));
- //Normals - Outlet
- para->setOutflowBoundaryNormalX(gridPath + input->getValue("outletBoundaryNormalX"));
- para->setOutflowBoundaryNormalY(gridPath + input->getValue("outletBoundaryNormalY"));
- para->setOutflowBoundaryNormalZ(gridPath + input->getValue("outletBoundaryNormalZ"));
- //////////////////////////////////////////////////////////////////////////
- //Forcing
- para->setForcing(StringUtil::toFloat(input->getValue("ForcingX")), StringUtil::toFloat(input->getValue("ForcingY")), StringUtil::toFloat(input->getValue("ForcingZ")));
- //////////////////////////////////////////////////////////////////////////
- //Particles
- para->setCalcParticles(StringUtil::toBool(input->getValue("calcParticles")));
- para->setParticleBasicLevel(StringUtil::toInt(input->getValue("baseLevel")));
- para->setParticleInitLevel(StringUtil::toInt(input->getValue("initLevel")));
- para->setNumberOfParticles(StringUtil::toInt(input->getValue("numberOfParticles")));
- para->setneighborWSB(gridPath + input->getValue("neighborWSB"));
- para->setStartXHotWall(StringUtil::toDouble(input->getValue("startXHotWall")));
- para->setEndXHotWall(StringUtil::toDouble(input->getValue("endXHotWall")));
- //////////////////////////////////////////////////////////////////////////
- //for Multi GPU
- if (para->getNumprocs() > 1)
- {
- ////////////////////////////////////////////////////////////////////////////
- ////1D domain decomposition
- //std::vector<std::string> sendProcNeighbors;
- //std::vector<std::string> recvProcNeighbors;
- //for (int i = 0; i<para->getNumprocs();i++)
- //{
- // sendProcNeighbors.push_back(gridPath + StringUtil::toString(i) + "s.dat");
- // recvProcNeighbors.push_back(gridPath + StringUtil::toString(i) + "r.dat");
- //}
- //para->setPossNeighborFiles(sendProcNeighbors, "send");
- //para->setPossNeighborFiles(recvProcNeighbors, "recv");
- //////////////////////////////////////////////////////////////////////////
- //3D domain decomposition
- std::vector<std::string> sendProcNeighborsX, sendProcNeighborsY, sendProcNeighborsZ;
- std::vector<std::string> recvProcNeighborsX, recvProcNeighborsY, recvProcNeighborsZ;
- for (int i = 0; i < para->getNumprocs(); i++)
- {
- sendProcNeighborsX.push_back(gridPath + StringUtil::toString(i) + "Xs.dat");
- sendProcNeighborsY.push_back(gridPath + StringUtil::toString(i) + "Ys.dat");
- sendProcNeighborsZ.push_back(gridPath + StringUtil::toString(i) + "Zs.dat");
- recvProcNeighborsX.push_back(gridPath + StringUtil::toString(i) + "Xr.dat");
- recvProcNeighborsY.push_back(gridPath + StringUtil::toString(i) + "Yr.dat");
- recvProcNeighborsZ.push_back(gridPath + StringUtil::toString(i) + "Zr.dat");
- }
- para->setPossNeighborFilesX(sendProcNeighborsX, "send");
- para->setPossNeighborFilesY(sendProcNeighborsY, "send");
- para->setPossNeighborFilesZ(sendProcNeighborsZ, "send");
- para->setPossNeighborFilesX(recvProcNeighborsX, "recv");
- para->setPossNeighborFilesY(recvProcNeighborsY, "recv");
- para->setPossNeighborFilesZ(recvProcNeighborsZ, "recv");
- }
- //////////////////////////////////////////////////////////////////////////
- //para->setkFull( input->getValue( "kFull" ));
- //para->setgeoFull( input->getValue( "geoFull" ));
- //para->setnoSlipBcPos( input->getValue( "noSlipBcPos" ));
- //para->setnoSlipBcQs( input->getValue( "noSlipBcQs" ));
- //para->setnoSlipBcValues( input->getValue( "noSlipBcValues" ));
- //para->setnoSlipBcValue( input->getValue( "noSlipBcValue" ));
- //para->setslipBcPos( input->getValue( "slipBcPos" ));
- //para->setslipBcQs( input->getValue( "slipBcQs" ));
- //para->setslipBcValue( input->getValue( "slipBcValue" ));
- //para->setpressBcPos( input->getValue( "pressBcPos" ));
- //para->setpressBcQs( input->getValue( "pressBcQs" ));
- //para->setpressBcValues( input->getValue( "pressBcValues" ));
- //para->setpressBcValue( input->getValue( "pressBcValue" ));
- //para->setvelBcQs( input->getValue( "velBcQs" ));
- //para->setvelBcValues( input->getValue( "velBcValues" ));
- //para->setpropellerCylinder( input->getValue( "propellerCylinder" ));
- //para->setpropellerQs( input->getValue( "propellerQs" ));
- //para->setwallBcQs( input->getValue( "wallBcQs" ));
- //para->setwallBcValues( input->getValue( "wallBcValues" ));
- //para->setperiodicBcQs( input->getValue( "periodicBcQs" ));
- //para->setperiodicBcValues( input->getValue( "periodicBcValues" ));
- //cout << "Try this: " << para->getgeomBoundaryBcValues() << endl;
-
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- //Restart
- para->setTimeDoCheckPoint(StringUtil::toInt(input->getValue("TimeDoCheckPoint")));
- para->setTimeDoRestart(StringUtil::toInt(input->getValue("TimeDoRestart")));
- para->setDoCheckPoint(StringUtil::toBool(input->getValue("DoCheckPoint")));
- para->setDoRestart(StringUtil::toBool(input->getValue("DoRestart")));
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- para->setMaxLevel(StringUtil::toInt(input->getValue("NOGL")));
- para->setGridX(StringUtil::toVector(input->getValue("GridX")));
- para->setGridY(StringUtil::toVector(input->getValue("GridY")));
- para->setGridZ(StringUtil::toVector(input->getValue("GridZ")));
- para->setDistX(StringUtil::toVector(input->getValue("DistX")));
- para->setDistY(StringUtil::toVector(input->getValue("DistY")));
- para->setDistZ(StringUtil::toVector(input->getValue("DistZ")));
-
- para->setNeedInterface(std::vector<bool>{true, true, true, true, true, true});
-}
-
-
-
-void multipleLevel(const std::string& configPath)
-{
- //std::ofstream logFile( "F:/Work/Computations/gridGenerator/grid/gridGeneratorLog.txt" );
- std::ofstream logFile( "grid/gridGeneratorLog.txt" );
- logging::Logger::addStream(&logFile);
-
- 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);
-
- //UbLog::reportingLevel() = UbLog::logLevelFromString("DEBUG5");
-
auto gridFactory = GridFactory::make();
gridFactory->setGridStrategy(Device::CPU);
- //gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::RAYCASTING);
gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_IN_OBJECT);
- //gridFactory->setTriangularMeshDiscretizationMethod(TriangularMeshDiscretizationMethod::POINT_UNDER_TRIANGLE);
auto gridBuilder = MultipleGridBuilder::makeShared(gridFactory);
-
- SPtr<Parameter> para = Parameter::make();
- SPtr<GridProvider> gridGenerator;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- bool useGridGenerator = true;
-
- if(useGridGenerator){
-
- enum testCase{
- DrivAer,
- DLC,
- MultiGPU,
- MetisTest
- };
-
- int testcase = MetisTest;
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- if( testcase == DrivAer )
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- {
- real dx = 0.2;
- real vx = 0.05;
-
- TriangularMesh* DrivAerSTL = TriangularMesh::make("F:/Work/Computations/gridGenerator/stl/DrivAer_Fastback_Coarse.stl");
- //TriangularMesh* triangularMesh = TriangularMesh::make("M:/TestGridGeneration/STL/DrivAer_NoSTLGroups.stl");
- //TriangularMesh* triangularMesh = TriangularMesh::make("M:/TestGridGeneration/STL/DrivAer_Coarse.stl");
- //TriangularMesh* DrivAerSTL = TriangularMesh::make("stl/DrivAer_Fastback_Coarse.stl");
-
- TriangularMesh* DrivAerRefBoxSTL = TriangularMesh::make("F:/Work/Computations/gridGenerator/stl/DrivAer_REF_BOX_Adrea.stl");
- //TriangularMesh* DrivAerRefBoxSTL = TriangularMesh::make("stl/DrivAer_REF_BOX_Adrea.stl");
-
- real z0 = 0.318+0.5*dx;
-
- gridBuilder->addCoarseGrid(- 5.0, -5.0, 0.0 - z0,
- 15.0, 5.0, 5.0 - z0, dx); // DrivAer
-
- //Object* floorBox = new Cuboid( -0.3, -1, -1, 4.0, 1, 0.2 );
- //Object* wakeBox = new Cuboid( 3.5, -1, -1, 5.5, 1, 0.8 );
-
- //Conglomerate* refRegion = new Conglomerate();
-
- //refRegion->add(floorBox);
- //refRegion->add(wakeBox);
- //refRegion->add(DrivAerRefBoxSTL);
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid(DrivAerRefBoxSTL, 2);
-
- //gridBuilder->setNumberOfLayers(10,8);
- //gridBuilder->addGrid(DrivAerSTL, 5);
-
- gridBuilder->addGeometry(DrivAerSTL);
-
- gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
- gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
-
- //////////////////////////////////////////////////////////////////////////
-
- 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);
-
- gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
- gridBuilder->setVelocityBoundaryCondition(SideType::MX, vx, 0.0, 0.0);
-
- gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0);
-
- //////////////////////////////////////////////////////////////////////////
-
- SPtr<Grid> grid = gridBuilder->getGrid(gridBuilder->getNumberOfLevels() - 1);
-
- gridBuilder->getGeometryBoundaryCondition(gridBuilder->getNumberOfLevels() - 1)->setTangentialVelocityForPatch( grid, 4, 0.0075, -2.0, 0.0,
- 0.0075, 2.0, 0.0, -vx, 0.318);
- gridBuilder->getGeometryBoundaryCondition(gridBuilder->getNumberOfLevels() - 1)->setTangentialVelocityForPatch( grid, 3, 2.793 , -2.0, 0.0,
- 2.793 , 2.0, 0.0, -vx, 0.318);
-
- //////////////////////////////////////////////////////////////////////////
-
- gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/DrivAer_Grid");
- gridBuilder->writeArrows ("F:/Work/Computations/gridGenerator/grid/DrivAer_Grid_arrow");
-
- //SimulationFileWriter::write("D:/GRIDGENERATION/files/", gridBuilder, FILEFORMAT::ASCII);
- //SimulationFileWriter::write("C:/Users/lenz/Desktop/Work/gridGenerator/grid/", gridBuilder, FILEFORMAT::ASCII);
- SimulationFileWriter::write("F:/Work/Computations/gridGenerator/grid/", gridBuilder, FILEFORMAT::BINARY);
- //SimulationFileWriter::write("grid/", gridBuilder, FILEFORMAT::ASCII);
-
- return;
-
- gridGenerator = GridGenerator::make(gridBuilder, para);
- }
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- if( testcase == DLC )
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- {
- real velocityRatio = 594.093427;
-
- real dx = 0.2;
- real vx = 0.065272188;
-
- real z0 = 0.24395 + 0.5*dx;
-
- std::vector<uint> ignorePatches = { 152, 153, 154 };
-
- //TriangularMesh* VW370_SERIE_STL = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/VW370_SERIE.stl", ignorePatches);
- TriangularMesh* VW370_SERIE_STL = TriangularMesh::make("stl/VW370_SERIE.stl", ignorePatches);
-
- //TriangularMesh* DLC_RefBox = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC_RefBox.stl");
-
- //TriangularMesh* DLC_RefBox_1 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC_RefBox_withWake/DLC_RefBox_withWake_4m.stl");
- //TriangularMesh* DLC_RefBox_2 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC_RefBox_withWake/DLC_RefBox_withWake_3m.stl");
- //TriangularMesh* DLC_RefBox_3 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC_RefBox_withWake/DLC_RefBox_withWake_2m.stl");
- //TriangularMesh* DLC_RefBox_4 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC_RefBox_withWake/DLC_RefBox_withWake_1m.stl");
-
- //TriangularMesh* DLC_RefBox_Level_3 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC/DLC_RefBox_Level_3.stl");
- //TriangularMesh* DLC_RefBox_Level_4 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC/DLC_RefBox_Level_4.stl");
- //TriangularMesh* DLC_RefBox_Level_5 = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/DLC/DLC_RefBox_Level_5.stl");
-
- TriangularMesh* DLC_RefBox_Level_3 = TriangularMesh::make("stl/DLC/DLC_RefBox_Level_3.stl");
- TriangularMesh* DLC_RefBox_Level_4 = TriangularMesh::make("stl/DLC/DLC_RefBox_Level_4.stl");
- TriangularMesh* DLC_RefBox_Level_5 = TriangularMesh::make("stl/DLC/DLC_RefBox_Level_5.stl");
-
- //TriangularMesh* VW370_SERIE_STL = TriangularMesh::make("stl/VW370_SERIE.stl", ignorePatches);
- //TriangularMesh* DLC_RefBox = TriangularMesh::make("stl/DLC_RefBox.lnx.stl");
- //TriangularMesh* DLC_RefBox_4 = TriangularMesh::make("stl/DLC_RefBox_withWake/DLC_RefBox_withWake_1m.lnx.stl");
-
- gridBuilder->addCoarseGrid(-30.0, -20.0, 0.0 - z0,
- 50.0, 20.0, 25.0 - z0, dx);
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid( new Cuboid( - 6.6, -6, -0.7, 20.6 , 6, 5.3 ), 1 );
- gridBuilder->addGrid( new Cuboid( -3.75, -3, -0.7, 11.75, 3, 2.65 ), 2 );
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid(DLC_RefBox_Level_3, 3);
- gridBuilder->addGrid(DLC_RefBox_Level_4, 4);
-
- Conglomerate* refinement = new Conglomerate();
- refinement->add(DLC_RefBox_Level_5);
- refinement->add(VW370_SERIE_STL);
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid(refinement, 5);
-
- gridBuilder->addGeometry(VW370_SERIE_STL);
-
- gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
- gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
-
- //////////////////////////////////////////////////////////////////////////
-
- 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);
-
- gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
- gridBuilder->setVelocityBoundaryCondition(SideType::MX, vx, 0.0, 0.0);
-
- gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0);
-
- //////////////////////////////////////////////////////////////////////////
-
- SPtr<Grid> grid = gridBuilder->getGrid(gridBuilder->getNumberOfLevels() - 1);
-
- real wheelsFrontX = -0.081;
- real wheelsRearX = 2.5486;
-
- real wheelsFrontZ = 0.0504;
- real wheelsRearZ = 0.057;
-
- real wheelsRadius = 0.318;
-
- real wheelRotationFrequency = 1170.74376 / 60.0;
-
- real wheelTangentialVelocity = -2.0 * M_PI * wheelsRadius * wheelRotationFrequency / velocityRatio;
-
- std::vector<uint> frontWheelPatches = { 71, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 159 };
- std::vector<uint> rearWheelPatches = { 82, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 160 };
-
- for( uint patch : frontWheelPatches ){
- gridBuilder->getGeometryBoundaryCondition(gridBuilder->getNumberOfLevels() - 1)->setTangentialVelocityForPatch( grid, patch, wheelsFrontX, -2.0, wheelsFrontZ,
- wheelsFrontX, 2.0, wheelsFrontZ,
- wheelTangentialVelocity, wheelsRadius);
- }
-
- for( uint patch : rearWheelPatches ){
- gridBuilder->getGeometryBoundaryCondition(gridBuilder->getNumberOfLevels() - 1)->setTangentialVelocityForPatch( grid, patch, wheelsRearX , -2.0, wheelsRearZ ,
- wheelsRearX , 2.0, wheelsRearZ ,
- wheelTangentialVelocity, wheelsRadius);
- }
-
- //////////////////////////////////////////////////////////////////////////
-
- //gridBuilder->writeGridsToVtk("C:/Users/lenz/Desktop/Work/gridGenerator/grid/DLC_Grid");
- //gridBuilder->writeArrows ("C:/Users/lenz/Desktop/Work/gridGenerator/grid/DLC_Grid_arrow");
-
- gridBuilder->writeGridsToVtk("grid/DLC_Grid");
- gridBuilder->writeArrows ("grid/DLC_Grid_arrow");
-
- //SimulationFileWriter::write("D:/GRIDGENERATION/files/", gridBuilder, FILEFORMAT::ASCII);
- //SimulationFileWriter::write("C:/Users/lenz/Desktop/Work/gridGenerator/grid/", gridBuilder, FILEFORMAT::ASCII);
- SimulationFileWriter::write("grid/", gridBuilder, FILEFORMAT::ASCII);
-
- gridGenerator = GridGenerator::make(gridBuilder, para);
- }
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- if( testcase == MultiGPU )
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- {
- //const uint generatePart = 1;
- const uint generatePart = Communicator::getInstanz()->getPID();
-
- std::ofstream logFile2;
-
- if( generatePart == 0 )
- //logFile2.open( "F:/Work/Computations/gridGenerator/grid/0/gridGeneratorLog.txt" );
- logFile2.open( "grid/0/gridGeneratorLog.txt" );
-
- if( generatePart == 1 )
- //logFile2.open( "F:/Work/Computations/gridGenerator/grid/1/gridGeneratorLog.txt" );
- logFile2.open( "grid/1/gridGeneratorLog.txt" );
-
- logging::Logger::addStream(&logFile2);
-
- real dx = 1.0 / 40.0;
- real vx = 0.05;
-
- //TriangularMesh* triangularMesh = TriangularMesh::make("F:/Work/Computations/gridGenerator/stl/ShpereNotOptimal.stl");
- TriangularMesh* triangularMesh = TriangularMesh::make("stl/ShpereNotOptimal.lnx.stl");
-
- // all
- //gridBuilder->addCoarseGrid(-2, -2, -2,
- // 4, 2, 2, dx);
-
- real overlap = 10.0 * dx;
-
- if( generatePart == 0 )
- gridBuilder->addCoarseGrid(-2.0 , -2.0, -2.0,
- 0.5 + overlap, 2.0, 2.0, dx);
-
- if( generatePart == 1 )
- gridBuilder->addCoarseGrid( 0.5 - overlap, -2.0, -2.0,
- 4.0 , 2.0, 2.0, dx);
-
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid(triangularMesh, 1);
-
- gridBuilder->addGeometry(triangularMesh);
-
- if( generatePart == 0 )
- gridBuilder->setSubDomainBox( std::make_shared<BoundingBox>( -2.0, 0.5,
- -2.0, 2.0,
- -2.0, 2.0 ) );
-
- if( generatePart == 1 )
- gridBuilder->setSubDomainBox( std::make_shared<BoundingBox>( 0.5, 4.0,
- -2.0, 2.0,
- -2.0, 2.0 ) );
-
- gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
- gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
-
- if( generatePart == 0 ){
- gridBuilder->findCommunicationIndices(CommunicationDirections::PX);
- gridBuilder->setCommunicationProcess(CommunicationDirections::PX, 1);
- }
-
- if( generatePart == 1 ){
- gridBuilder->findCommunicationIndices(CommunicationDirections::MX);
- gridBuilder->setCommunicationProcess(CommunicationDirections::MX, 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);
-
- if (generatePart == 0) {
- gridBuilder->setVelocityBoundaryCondition(SideType::MX, vx, 0.0, 0.0);
- }
- if (generatePart == 1) {
- gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
- }
-
- gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0);
-
- //////////////////////////////////////////////////////////////////////////
-
- if (generatePart == 0) {
- //gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/0/Test_");
- //gridBuilder->writeArrows ("F:/Work/Computations/gridGenerator/grid/0/Test_Arrow");
- }
- if (generatePart == 1) {
- //gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/1/Test_");
- //gridBuilder->writeArrows ("F:/Work/Computations/gridGenerator/grid/1/Test_Arrow");
- }
-
- if (generatePart == 0)
- //SimulationFileWriter::write("F:/Work/Computations/gridGenerator/grid/0/", gridBuilder, FILEFORMAT::ASCII);
- SimulationFileWriter::write("grid/0/", gridBuilder, FILEFORMAT::ASCII);
- if (generatePart == 1)
- //SimulationFileWriter::write("F:/Work/Computations/gridGenerator/grid/1/", gridBuilder, FILEFORMAT::ASCII);
- SimulationFileWriter::write("grid/1/", gridBuilder, FILEFORMAT::ASCII);
-
- //return;
-
- gridGenerator = GridGenerator::make(gridBuilder, para);
- }
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- if( testcase == MetisTest )
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- {
- //const uint generatePart = 1;
- const uint generatePart = Communicator::getInstanz()->getPID();
-
- real dx = 1.0 / 20.0;
- real vx = 0.05;
-
- TriangularMesh* triangularMesh = TriangularMesh::make("F:/Work/Computations/gridGenerator/stl/ShpereNotOptimal.stl");
- //TriangularMesh* triangularMesh = TriangularMesh::make("stl/ShpereNotOptimal.lnx.stl");
-
- // all
- //gridBuilder->addCoarseGrid(-2, -2, -2,
- // 4, 2, 2, dx);
-
- real overlap = 10.0 * dx;
-
- gridBuilder->addCoarseGrid(-2.0, -2.0, -2.0,
- 4.0, 2.0, 2.0, dx);
-
-
- gridBuilder->setNumberOfLayers(10,8);
- gridBuilder->addGrid(triangularMesh, 1);
-
- gridBuilder->addGeometry(triangularMesh);
-
- gridBuilder->setPeriodicBoundaryCondition(false, false, false);
-
- gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
-
- //////////////////////////////////////////////////////////////////////////
- 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);
+ real dx = 1.0 / 32.0;
- gridBuilder->setVelocityBoundaryCondition(SideType::MX, vx, 0.0, 0.0);
- gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
+ gridBuilder->addCoarseGrid(-0.5, -0.5, -0.5,
+ 0.5, 0.5, 0.5, dx);
- gridBuilder->setVelocityBoundaryCondition(SideType::GEOMETRY, 0.0, 0.0, 0.0);
-
- //////////////////////////////////////////////////////////////////////////
+ Cuboid cube(-1.0, -1.0, 0.45, 1.0, 1.0, 0.55);
- gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/Test_");
- //gridBuilder->writeArrows ("F:/Work/Computations/gridGenerator/grid/Test_Arrow");
+ gridBuilder->setNumberOfLayers(10,8);
+ gridBuilder->addGrid( &cube, 1);
- //SimulationFileWriter::write("F:/Work/Computations/gridGenerator/grid/", gridBuilder, FILEFORMAT::ASCII);
+ gridBuilder->setPeriodicBoundaryCondition(false, false, false);
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- if(false)
- {
+ gridBuilder->buildGrids(GKS, false);
- auto getParentIndex = [&] (uint index, uint level) -> uint
- {
- SPtr<Grid> grid = gridBuilder->getGrid( level );
+ gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/Grid_lev_");
- if( level != 0 )
- {
- real x, y, z;
- grid->transIndexToCoords(index, x, y, z);
-
- SPtr<Grid> coarseGrid = gridBuilder->getGrid(level - 1);
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- if (std::abs(dir[0]) < 0.5 || std::abs(dir[1]) < 0.5 || std::abs(dir[2]) < 0.5) continue;
-
- real coarseX = x + dir[0] * 0.5 * grid->getDelta();
- real coarseY = y + dir[1] * 0.5 * grid->getDelta();
- real coarseZ = z + dir[2] * 0.5 * grid->getDelta();
-
- // check if close enough to coarse grid coordinates
- if( 0.01 * grid->getDelta() < std::abs( (coarseGrid->getStartX() - coarseX) / grid->getDelta()
- - lround( (coarseGrid->getStartX() - coarseX) / grid->getDelta() ) ) ) continue;
- if( 0.01 * grid->getDelta() < std::abs( (coarseGrid->getStartY() - coarseY) / grid->getDelta()
- - lround( (coarseGrid->getStartY() - coarseY) / grid->getDelta() ) ) ) continue;
- if( 0.01 * grid->getDelta() < std::abs( (coarseGrid->getStartZ() - coarseZ) / grid->getDelta()
- - lround( (coarseGrid->getStartZ() - coarseZ) / grid->getDelta() ) ) ) continue;
-
- uint parentIndex = coarseGrid->transCoordToIndex( coarseX, coarseY, coarseZ);
-
- return parentIndex;
- }
- }
-
- return INVALID_INDEX;
- };
-
-
- std::vector<idx_t> xadj;
- std::vector<idx_t> adjncy;
-
- std::vector<idx_t> vwgt;
- std::vector<idx_t> adjwgt;
-
- idx_t vertexCounter = 0;
- uint edgeCounter = 0;
-
- std::cout << "Checkpoint 1:" << std::endl;
-
- std::vector< std::vector<idx_t> > vertexIndex( gridBuilder->getNumberOfLevels() );
-
- std::vector< uint > startVerticesPerLevel;;
-
- for( uint level = 0; level < gridBuilder->getNumberOfLevels(); level++ )
- {
- SPtr<Grid> grid = gridBuilder->getGrid( level );
-
- vertexIndex[level].resize( grid->getSize() );
-
- startVerticesPerLevel.push_back(vertexCounter);
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- if (grid->getSparseIndex(index) == INVALID_INDEX)
- {
- vertexIndex[level][index] = INVALID_INDEX;
- continue;
- }
-
- uint parentIndex = getParentIndex(index, level);
-
- if( parentIndex != INVALID_INDEX )
- {
- SPtr<Grid> coarseGrid = gridBuilder->getGrid(level - 1);
-
- if( coarseGrid->getFieldEntry(parentIndex) == FLUID_CFC ||
- coarseGrid->getFieldEntry(parentIndex) == FLUID_FCC ||
- coarseGrid->getFieldEntry(parentIndex) == STOPPER_COARSE_UNDER_FINE )
- {
- //vertexIndex[level][index] = INVALID_INDEX;
- vertexIndex[level][index] = vertexIndex[level - 1][parentIndex];
- continue;
- }
- }
-
- vertexIndex[level][index] = vertexCounter;
-
- vwgt.push_back( std::pow(2, level) );
- //vwgt.push_back( std::pow(2, 2*level) );
- vertexCounter++;
- }
-
- }
-
- //////////////////////////////////////////////////////////////////////////
- //for( uint level = 0; level < gridBuilder->getNumberOfLevels(); level++ )
- //{
- // SPtr<Grid> grid = gridBuilder->getGrid( level );
-
- // for (uint index = 0; index < grid->getSize(); index++)
- // {
- // grid->setFieldEntry(index, vertexIndex[level][index] >= startVerticesPerLevel[level] && vertexIndex[level][index] != INVALID_INDEX);
- // }
- //}
-
- //gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/VertexIndex_");
-
- //return;
- //////////////////////////////////////////////////////////////////////////
-
-
- std::cout << "Checkpoint 2:" << std::endl;
-
- for( uint level = 0; level < gridBuilder->getNumberOfLevels(); level++ )
- {
- SPtr<Grid> grid = gridBuilder->getGrid( level );
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- //if (grid->getSparseIndex(index) == INVALID_INDEX) continue;
-
- if( vertexIndex[level][index] == INVALID_INDEX ) continue;
-
- if( vertexIndex[level][index] < startVerticesPerLevel[level] ) continue;
-
- xadj.push_back(edgeCounter);
-
- real x, y, z;
- grid->transIndexToCoords(index, x, y, z);
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- const uint neighborIndex = grid->transCoordToIndex(x + dir[0] * grid->getDelta(),
- y + dir[1] * grid->getDelta(),
- z + dir[2] * grid->getDelta());
-
- if (neighborIndex == INVALID_INDEX) continue;
-
- if (neighborIndex == index) continue;
-
- if( vertexIndex[level][neighborIndex] == INVALID_INDEX ) continue;
-
- adjncy.push_back( vertexIndex[level][neighborIndex] );
- adjwgt.push_back( std::pow(2, level) );
-
- edgeCounter++;
- }
-
- if( grid->getFieldEntry(index) == FLUID_CFC ||
- grid->getFieldEntry(index) == FLUID_FCC ||
- grid->getFieldEntry(index) == STOPPER_COARSE_UNDER_FINE )
-
- {
- SPtr<Grid> fineGrid = gridBuilder->getGrid(level + 1);
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- if (std::abs(dir[0]) < 0.5 || std::abs(dir[1]) < 0.5 || std::abs(dir[2]) < 0.5) continue;
-
- real fineX = x + dir[0] * 0.25 * grid->getDelta();
- real fineY = y + dir[1] * 0.25 * grid->getDelta();
- real fineZ = z + dir[2] * 0.25 * grid->getDelta();
-
- uint childIndex = fineGrid->transCoordToIndex(fineX, fineY, fineZ);
-
- if( fineGrid->getFieldEntry(childIndex) == INVALID_INDEX ) continue;
- if( vertexIndex[level + 1][childIndex] == INVALID_INDEX ) continue;
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- const uint neighborIndex = fineGrid->transCoordToIndex( fineX + dir[0] * fineGrid->getDelta(),
- fineY + dir[1] * fineGrid->getDelta(),
- fineZ + dir[2] * fineGrid->getDelta() );
-
- if(neighborIndex == INVALID_INDEX) continue;
-
- if (neighborIndex == childIndex) continue;
-
- if( vertexIndex[level + 1][neighborIndex] == INVALID_INDEX ) continue;
-
- adjncy.push_back( vertexIndex[level + 1][neighborIndex] );
- adjwgt.push_back( std::pow(2, level) );
-
- edgeCounter++;
- }
- }
- }
- }
- }
-
- xadj.push_back( edgeCounter );
-
- std::cout << "Checkpoint 3:" << std::endl;
-
- idx_t nWeights = 1;
- idx_t nParts = 4;
- idx_t objval = 0;
-
- std::vector<idx_t> part( vertexCounter );
-
- std::cout << vertexCounter << std::endl;
- std::cout << edgeCounter << std::endl;
- std::cout << xadj.size() << std::endl;
- std::cout << adjncy.size() << std::endl;
-
- //int ret = METIS_PartGraphRecursive(&vertexCounter, &nWeights, xadj.data(), adjncy.data(),
- // vwgt.data(), NULL, adjwgt.data(), &nParts,
- // NULL, NULL, NULL, &objval, part.data());
-
- int ret = METIS_PartGraphKway(&vertexCounter, &nWeights, xadj.data(), adjncy.data(),
- vwgt.data(), NULL, NULL/*adjwgt.data()*/, &nParts,
- NULL, NULL, NULL, &objval, part.data());
-
- std::cout << "objval:" << objval << std::endl;
-
- std::cout << "Checkpoint 4:" << std::endl;
-
- //uint partCounter = 0;
-
- for( uint level = 0; level < gridBuilder->getNumberOfLevels(); level++ )
- {
- SPtr<Grid> grid = gridBuilder->getGrid( level );
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- if (grid->getSparseIndex(index) == INVALID_INDEX) continue;
-
- grid->setFieldEntry(index, part[vertexIndex[level][index]]);
-
- //partCounter++;
- }
- }
-
- std::cout << "Checkpoint 5:" << std::endl;
-
- gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/Partition_");
-
- }
-
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- {
-
- for( int level = gridBuilder->getNumberOfLevels()-1; level >= 0 ; level-- )
- {
- std::vector< std::vector<idx_t> > vertexIndex( gridBuilder->getNumberOfLevels() );
-
- std::vector<idx_t> xadj;
- std::vector<idx_t> adjncy;
-
- std::vector<idx_t> vwgt;
- std::vector<idx_t> adjwgt;
-
- idx_t vertexCounter = 0;
- uint edgeCounter = 0;
-
- SPtr<Grid> grid = gridBuilder->getGrid( level );
-
- vertexIndex[level].resize( grid->getSize() );
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- if (grid->getSparseIndex(index) == INVALID_INDEX)
- {
- vertexIndex[level][index] = INVALID_INDEX;
- continue;
- }
-
- vertexIndex[level][index] = vertexCounter;
-
- vwgt.push_back( std::pow(2, level) );
- //vwgt.push_back( std::pow(2, 2*level) );
- vertexCounter++;
- }
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- //if (grid->getSparseIndex(index) == INVALID_INDEX) continue;
-
- if( vertexIndex[level][index] == INVALID_INDEX ) continue;
-
- xadj.push_back(edgeCounter);
-
- real x, y, z;
- grid->transIndexToCoords(index, x, y, z);
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- const uint neighborIndex = grid->transCoordToIndex(x + dir[0] * grid->getDelta(),
- y + dir[1] * grid->getDelta(),
- z + dir[2] * grid->getDelta());
-
- if (neighborIndex == INVALID_INDEX) continue;
-
- if (neighborIndex == index) continue;
-
- if( vertexIndex[level][neighborIndex] == INVALID_INDEX ) continue;
-
- adjncy.push_back( vertexIndex[level][neighborIndex] );
- adjwgt.push_back( std::pow(2, level) );
-
- edgeCounter++;
- }
- }
-
- xadj.push_back( edgeCounter );
-
- std::cout << "Checkpoint 3:" << std::endl;
-
- idx_t nWeights = 1;
- idx_t nParts = 4;
- idx_t objval = 0;
-
- std::vector<idx_t> part( vertexCounter );
-
- std::cout << vertexCounter << std::endl;
- std::cout << edgeCounter << std::endl;
- std::cout << xadj.size() << std::endl;
- std::cout << adjncy.size() << std::endl;
-
- int ret = METIS_PartGraphRecursive(&vertexCounter, &nWeights, xadj.data(), adjncy.data(),
- NULL/*vwgt.data()*/, NULL, NULL/*adjwgt.data()*/, &nParts,
- NULL, NULL, NULL, &objval, part.data());
-
- //int ret = METIS_PartGraphKway(&vertexCounter, &nWeights, xadj.data(), adjncy.data(),
- // NULL/*vwgt.data()*/, NULL, NULL/*adjwgt.data()*/, &nParts,
- // NULL, NULL, NULL, &objval, part.data());
-
- std::cout << "objval:" << objval << std::endl;
-
- std::cout << "Checkpoint 4:" << std::endl;
-
- for (uint index = 0; index < grid->getSize(); index++)
- {
- if (vertexIndex[level][index] == INVALID_INDEX) continue;
-
- if( grid->getFieldEntry(index) == FLUID_CFC ||
- grid->getFieldEntry(index) == FLUID_FCC ||
- grid->getFieldEntry(index) == STOPPER_COARSE_UNDER_FINE )
- {
- SPtr<Grid> fineGrid = gridBuilder->getGrid(level+1);
-
- real x, y, z;
- grid->transIndexToCoords(index, x, y, z);
-
- for (const auto dir : DistributionHelper::getDistribution27())
- {
- if (std::abs(dir[0]) < 0.5 || std::abs(dir[1]) < 0.5 || std::abs(dir[2]) < 0.5) continue;
-
- real fineX = x + dir[0] * 0.25 * grid->getDelta();
- real fineY = y + dir[1] * 0.25 * grid->getDelta();
- real fineZ = z + dir[2] * 0.25 * grid->getDelta();
-
- uint childIndex = fineGrid->transCoordToIndex(fineX, fineY, fineZ);
-
- if( childIndex == INVALID_INDEX ) continue;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- fineGrid->setFieldEntry(childIndex, part[vertexIndex[level][index]]);
- //fineGrid->setFieldEntry(childIndex, grid->getFieldEntry(index));
- }
- }
+ GksMeshAdapter meshAdapter;
- grid->setFieldEntry(index, part[vertexIndex[level][index]]);
- }
- }
+ meshAdapter.inputGrid( gridBuilder );
- gridBuilder->writeGridsToVtk("F:/Work/Computations/gridGenerator/grid/Partition_");
+ meshAdapter.findQuadtreeConnectivity( gridBuilder );
- }
+ meshAdapter.findCellToCellConnectivity( gridBuilder );
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ meshAdapter.countCells();
- return;
+ meshAdapter.partitionCells();
- gridGenerator = GridGenerator::make(gridBuilder, para);
- }
- }
- else
- {
- gridGenerator = GridReader::make(FileFormat::BINARY, para);
- //gridGenerator = GridReader::make(FileFormat::ASCII, para);
- }
+ meshAdapter.generateNodes( gridBuilder );
- logFile.close();
+ meshAdapter.computeCellGeometry();
- //return;
-
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ meshAdapter.generateFaces( gridBuilder );
+ meshAdapter.sortFaces();
- std::ifstream stream;
- stream.open(configPath.c_str(), std::ios::in);
- if (stream.fail())
- throw std::runtime_error("can not open config file!");
+ meshAdapter.countFaces();
- UPtr<input::Input> input = input::Input::makeInput(stream, "config");
+ meshAdapter.generateInterfaceConnectivity();
- setParameters(para, input);
+ meshAdapter.writeMeshVTK( "F:/Work/Computations/gridGenerator/grid/Mesh.vtk" );
- Simulation sim;
- SPtr<FileWriter> fileWriter = SPtr<FileWriter>(new FileWriter());
- sim.init(para, gridGenerator, fileWriter);
- sim.run();
+ meshAdapter.writeMeshFaceVTK( "F:/Work/Computations/gridGenerator/grid/MeshFaces.vtk" );
}
-
int main( int argc, char* argv[])
{
- MPI_Init(&argc, &argv);
- std::string str, str2;
- if ( argv != NULL )
- {
- str = static_cast<std::string>(argv[0]);
- if (argc > 1)
- {
- str2 = static_cast<std::string>(argv[1]);
- try
- {
- multipleLevel(str2);
- }
- catch (const std::exception& e)
- {
- *logging::out << logging::Logger::ERROR << e.what() << "\n";
- //MPI_Abort(MPI_COMM_WORLD, -1);
- }
- catch (...)
- {
- std::cout << "unknown exeption" << std::endl;
- }
- }
- else
- {
- try
- {
- multipleLevel("F:/Work/Computations/gridGenerator/inp/configTest.txt");
- }
- catch (const std::exception& e)
- {
-
- *logging::out << logging::Logger::ERROR << e.what() << "\n";
- //std::cout << e.what() << std::flush;
- //MPI_Abort(MPI_COMM_WORLD, -1);
- }
- catch (const std::bad_alloc e)
- {
-
- *logging::out << logging::Logger::ERROR << "Bad Alloc:" << e.what() << "\n";
- //std::cout << e.what() << std::flush;
- //MPI_Abort(MPI_COMM_WORLD, -1);
- }
- catch (...)
- {
- *logging::out << logging::Logger::ERROR << "Unknown exception!\n";
- //std::cout << "unknown exeption" << std::endl;
- }
- std::cout << "\nConfiguration file must be set!: lbmgm <config file>" << std::endl << std::flush;
- //MPI_Abort(MPI_COMM_WORLD, -1);
- }
+ logging::Logger::addStream(&std::cout);
+ logging::Logger::setDebugLevel(logging::Logger::Level::INFO_LOW);
+ logging::Logger::timeStamp(logging::Logger::ENABLE);
+
+ try
+ {
+ gksTest();
+ }
+ catch (const std::exception& e)
+ {
+ *logging::out << logging::Logger::ERROR << e.what() << "\n";
+ }
+ catch (const std::bad_alloc& e)
+ {
+ *logging::out << logging::Logger::ERROR << "Bad Alloc:" << e.what() << "\n";
+ }
+ catch (...)
+ {
+ *logging::out << logging::Logger::ERROR << "Unknown exception!\n";
}
-
- /*
- MPE_Init_log() & MPE_Finish_log() are NOT needed when
- liblmpe.a is linked with this program. In that case,
- MPI_Init() would have called MPE_Init_log() already.
- */
-#if defined( MPI_LOGGING )
- MPE_Init_log();
-#endif
-
-#if defined( MPI_LOGGING )
- if ( argv != NULL )
- MPE_Finish_log( argv[0] );
- if ( str != "" )
- MPE_Finish_log( str.c_str() );
- else
- MPE_Finish_log( "TestLog" );
-#endif
-
- MPI_Finalize();
return 0;
}
diff --git a/targets/libs/GksMeshAdapter/3rdPartyLinking.cmake b/targets/libs/GksMeshAdapter/3rdPartyLinking.cmake
new file mode 100644
index 000000000..6feb9029e
--- /dev/null
+++ b/targets/libs/GksMeshAdapter/3rdPartyLinking.cmake
@@ -0,0 +1,6 @@
+include (${CMAKE_SOURCE_DIR}/${cmakeMacroPath}/Cuda/Link.cmake)
+linkCuda(${targetName})
+include (${CMAKE_SOURCE_DIR}/${cmakeMacroPath}/MPI/Link.cmake)
+linkMPI(${targetName})
+include (${CMAKE_SOURCE_DIR}/${cmakeMacroPath}/Boost/Link.cmake)
+linkBoost(${targetName} "serialization")
\ No newline at end of file
diff --git a/targets/libs/GksMeshAdapter/CMakeLists.txt b/targets/libs/GksMeshAdapter/CMakeLists.txt
new file mode 100644
index 000000000..fadfa8657
--- /dev/null
+++ b/targets/libs/GksMeshAdapter/CMakeLists.txt
@@ -0,0 +1,18 @@
+setTargetNameToFolderName(${CMAKE_CURRENT_LIST_DIR})
+
+set(linkDirectories "")
+set(libsToLink GridGenerator Core)
+
+set(includeDirectories ${CMAKE_SOURCE_DIR}/src/${targetName}
+ ${CMAKE_SOURCE_DIR}/src
+ ${CMAKE_SOURCE_DIR}/src/GridGenerator
+ ${CMAKE_SOURCE_DIR}/src/Core )
+
+#glob files and save in MY_SRCS
+include(CMakePackage.cmake)
+
+buildLib(${targetName} "${MY_SRCS}" "${linkDirectories}" "${libsToLink}" "${includeDirectories}")
+groupTarget(${targetName} ${gksLibraryFolder})
+
+#Specify the linking to 3rdParty libs
+include(3rdPartyLinking.cmake)
diff --git a/targets/libs/GksMeshAdapter/CMakePackage.cmake b/targets/libs/GksMeshAdapter/CMakePackage.cmake
new file mode 100644
index 000000000..7ec316fe0
--- /dev/null
+++ b/targets/libs/GksMeshAdapter/CMakePackage.cmake
@@ -0,0 +1,8 @@
+#FILE ENDINGS
+resetFileEndingsToCollect()
+addCAndCPPFileTypes()
+
+#GLOB SOURCE FILES IN MY_SRCS
+unset(MY_SRCS)
+includeRecursiveAllFilesFrom(${targetName} ${CMAKE_CURRENT_LIST_DIR})
+includeRecursiveProductionFilesFrom(${targetName} ${CMAKE_SOURCE_DIR}/src/${targetName})
\ No newline at end of file
--
GitLab