diff --git a/apps/2020-scaling-workshop/Scaling_Workshop_01_Cube.cpp b/apps/2020-scaling-workshop/Scaling_Workshop_01_Cube.cpp index 4b1452417b6bae4aa4024a8ad93823ab83ae4054..41a9b297dcee60add74d597e30dd80922e6b08ba 100644 --- a/apps/2020-scaling-workshop/Scaling_Workshop_01_Cube.cpp +++ b/apps/2020-scaling-workshop/Scaling_Workshop_01_Cube.cpp @@ -93,11 +93,13 @@ void benchmark( int argc, char** argv ) auto meshInfo = MeshInfo::meshSymmetricCuboid( leftBottom3D, Point3D( { 1, 1, 1 } ), numEdgesPerSide, numEdgesPerSide, numEdgesPerSide ); - auto onBoundary = []( const Point3D& ) { return true; }; - meshInfo.setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); - SetupPrimitiveStorage setupStorage( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); - setupStorage.setMeshBoundaryFlagsInner( 0, true ); + // new code ... + setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true ); + // ... replaces + // auto onBoundary = []( const Point3D& ) { return true; }; + // meshInfo.setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); + // setupStorage.setMeshBoundaryFlagsInner( 0, true ); storage = std::make_shared< PrimitiveStorage >( setupStorage ); } @@ -206,4 +208,4 @@ void benchmark( int argc, char** argv ) int main( int argc, char** argv ) { hyteg::scaling_workshop::benchmark( argc, argv ); -} \ No newline at end of file +} diff --git a/apps/2020-scaling-workshop/Scaling_Workshop_02_Spherical_Shell.cpp b/apps/2020-scaling-workshop/Scaling_Workshop_02_Spherical_Shell.cpp index 31ed962c564751d5f1b103b5baf2ad531bc05702..42b923a74fe1947c2dd16bcd12136d783ae5aab0 100644 --- a/apps/2020-scaling-workshop/Scaling_Workshop_02_Spherical_Shell.cpp +++ b/apps/2020-scaling-workshop/Scaling_Workshop_02_Spherical_Shell.cpp @@ -92,11 +92,13 @@ void benchmark( int argc, char** argv ) { auto meshInfo = MeshInfo::meshSphericalShell( ntan, nrad, rmin, rmax, MeshInfo::SHELLMESH_ON_THE_FLY ); - auto onBoundary = []( const Point3D& ) { return true; }; - meshInfo.setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); - SetupPrimitiveStorage setupStorage( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); - setupStorage.setMeshBoundaryFlagsInner( 0, true ); + // new code ... + setupStorage.setMeshBoundaryFlagsOnBoundary( 1, 0, true ); + // ... replaces + // auto onBoundary = []( const Point3D& ) { return true; }; + // meshInfo.setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); + // setupStorage.setMeshBoundaryFlagsInner( 0, true ); storage = std::make_shared< PrimitiveStorage >( setupStorage ); } @@ -270,4 +272,4 @@ void benchmark( int argc, char** argv ) int main( int argc, char** argv ) { hyteg::scaling_workshop::benchmark( argc, argv ); -} \ No newline at end of file +} diff --git a/apps/2020-tme/Benchmark_01_Cube.cpp b/apps/2020-tme/Benchmark_01_Cube.cpp index 93b88d66d6e1d8ec48520926e4eace5323bbc63f..96699db0bf63d0082b4b9ea63d71c70ee78e0bbe 100644 --- a/apps/2020-tme/Benchmark_01_Cube.cpp +++ b/apps/2020-tme/Benchmark_01_Cube.cpp @@ -105,11 +105,11 @@ void benchmark( int argc, char** argv ) auto meshInfo = MeshInfo::meshSymmetricCuboid( leftBottom3D, Point3D( { 1, 1, 1 } ), numEdgesPerSide, numEdgesPerSide, numEdgesPerSide ); - auto onBoundary = []( const Point3D& ) { return true; }; - meshInfo.setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); - auto setupStorage = std::make_shared< SetupPrimitiveStorage >( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); + + auto onBoundary = []( const Point3D& ) { return true; }; + setupStorage->setMeshBoundaryFlagsByVertexLocation( 1, onBoundary ); setupStorage->setMeshBoundaryFlagsInner( 0, true ); auto storage = std::make_shared< PrimitiveStorage >( *setupStorage ); diff --git a/apps/2020-tme/Benchmark_02_Y-Pipe.cpp b/apps/2020-tme/Benchmark_02_Y-Pipe.cpp index e778049df375ee3b61526036e4ea42b54791afdb..743dfb4dcf1e22c1c4cd42c002e6d6b51728ec07 100644 --- a/apps/2020-tme/Benchmark_02_Y-Pipe.cpp +++ b/apps/2020-tme/Benchmark_02_Y-Pipe.cpp @@ -162,9 +162,6 @@ void benchmark( int argc, char** argv ) return false; }; - meshInfo.setAllMeshBoundaryFlags( 1 ); - meshInfo.setMeshBoundaryFlagsByVertexLocation( 2, outflowBoundary, false ); - auto tiltZCoordinateMap = [&]( const hyteg::Point3D& p ) -> Point3D { if ( p[0] > 1 + eps && p[1] > -eps ) @@ -257,7 +254,15 @@ void benchmark( int argc, char** argv ) auto setupStorage = std::make_shared< SetupPrimitiveStorage >( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); - setupStorage->setMeshBoundaryFlagsInner( 0, true ); + + // new code ... + setupStorage->setMeshBoundaryFlagsOnBoundary( 1, 0, true ); + setupStorage->setMeshBoundaryFlagsByVertexLocation( 2, outflowBoundary, false ); + + // ... replaces old + // meshInfo.setAllMeshBoundaryFlags( 1 ); + // meshInfo.setMeshBoundaryFlagsByVertexLocation( 2, outflowBoundary, false ); + // setupStorage->setMeshBoundaryFlagsInner( 0, true ); auto storage = std::make_shared< PrimitiveStorage >( *setupStorage ); diff --git a/apps/PolarLaplacian.cpp b/apps/PolarLaplacian.cpp index c6908e2e6c15fc069a9a64ff7a65e94ea1b8a2c5..1cd487df9182a418cfdf41f1233ce13d289db3d5 100644 --- a/apps/PolarLaplacian.cpp +++ b/apps/PolarLaplacian.cpp @@ -127,18 +127,8 @@ void solve_using_geometry_map( MeshInfo& meshInfo, walberla::Config::BlockHandle // Prepare storage and set geometry mapping SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); + PolarCoordsMap::setMap( setupStorage ); - for( auto it : setupStorage.getFaces() ) { - setupStorage.setGeometryMap( it.second->getID(), std::make_shared< PolarCoordsMap >() ); - } - - for( auto it : setupStorage.getEdges() ) { - setupStorage.setGeometryMap( it.second->getID(), std::make_shared< PolarCoordsMap >() ); - } - - for( auto it : setupStorage.getVertices() ) { - setupStorage.setGeometryMap( it.second->getID(), std::make_shared< PolarCoordsMap >() ); - } loadbalancing::roundRobin( setupStorage ); std::shared_ptr< PrimitiveStorage > storage = std::make_shared< PrimitiveStorage >( setupStorage ); diff --git a/src/hyteg/geometry/PolarCoordsMap.hpp b/src/hyteg/geometry/PolarCoordsMap.hpp index 98bb1abfd4b504f9eabdf20465836e7d8eb329b7..c56cc3c1e575ab708b538fbe1fa0f0058e9b4a1e 100644 --- a/src/hyteg/geometry/PolarCoordsMap.hpp +++ b/src/hyteg/geometry/PolarCoordsMap.hpp @@ -20,42 +20,64 @@ #pragma once #include -#include "GeometryMap.hpp" + +#include "hyteg/geometry/GeometryMap.hpp" +#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp" namespace hyteg { - /// Class providing geometry mapping based on polar coordinates; convention is x[0] = \f$\rho\f$, x[1] = \f$\varphi\f$ - class PolarCoordsMap : public GeometryMap - { - public: +/// Class providing geometry mapping based on polar coordinates; convention is x[0] = \f$\rho\f$, x[1] = \f$\varphi\f$ +class PolarCoordsMap : public GeometryMap +{ + public: + PolarCoordsMap() {} + + static void setMap( SetupPrimitiveStorage& setupStorage ) + { + auto blend = std::make_shared< PolarCoordsMap >(); + + for ( auto it : setupStorage.getFaces() ) + { + Face& face = *it.second; + setupStorage.setGeometryMap( face.getID(), blend ); + } + + for ( auto it : setupStorage.getEdges() ) + { + Edge& edge = *it.second; + setupStorage.setGeometryMap( edge.getID(), blend ); + } - PolarCoordsMap(){} + for ( auto it : setupStorage.getVertices() ) + { + Vertex& vertex = *it.second; + setupStorage.setGeometryMap( vertex.getID(), blend ); + } + } - void evalF( const Point3D& x, Point3D& Fx ) const { + void evalF( const Point3D& x, Point3D& Fx ) const + { Fx[0] = x[0] * std::cos( x[1] ); Fx[1] = x[0] * std::sin( x[1] ); - } - - void evalDF( const Point3D& x, Matrix2r& DFx ) const - { - DFx( 0, 0 ) = std::cos( x[1] ); - DFx( 0, 1 ) = - x[0] * std::sin( x[1] ); - DFx( 1, 0 ) = std::sin( x[1] ); - DFx( 1, 1 ) = x[0] * std::cos( x[1] ); - } - - void evalDFinv( const Point3D& x, Matrix2r& DFinvx ) const - { - DFinvx( 0, 0 ) = std::cos( x[1] ); - DFinvx( 0, 1 ) = std::sin( x[1] ); - DFinvx( 1, 0 ) = - std::sin( x[1] ) / x[0]; - DFinvx( 1, 1 ) = std::cos( x[1] ) / x[0]; - } - - void serializeSubClass( walberla::mpi::SendBuffer& sendBuffer ) const { - sendBuffer << Type::POLAR_COORDS; - } - - }; + } + + void evalDF( const Point3D& x, Matrix2r& DFx ) const + { + DFx( 0, 0 ) = std::cos( x[1] ); + DFx( 0, 1 ) = -x[0] * std::sin( x[1] ); + DFx( 1, 0 ) = std::sin( x[1] ); + DFx( 1, 1 ) = x[0] * std::cos( x[1] ); + } + + void evalDFinv( const Point3D& x, Matrix2r& DFinvx ) const + { + DFinvx( 0, 0 ) = std::cos( x[1] ); + DFinvx( 0, 1 ) = std::sin( x[1] ); + DFinvx( 1, 0 ) = -std::sin( x[1] ) / x[0]; + DFinvx( 1, 1 ) = std::cos( x[1] ) / x[0]; + } + + void serializeSubClass( walberla::mpi::SendBuffer& sendBuffer ) const { sendBuffer << Type::POLAR_COORDS; } +}; } // namespace hyteg diff --git a/src/hyteg/mesh/MeshGenAnnulus.cpp b/src/hyteg/mesh/MeshGenAnnulus.cpp index 413759a318026f5d3a8f863bd2111bd4507b509e..0d4b293934341cdf8b609ace85fd70a1b976b37f 100644 --- a/src/hyteg/mesh/MeshGenAnnulus.cpp +++ b/src/hyteg/mesh/MeshGenAnnulus.cpp @@ -33,6 +33,11 @@ using walberla::math::pi; namespace hyteg { +// boundary flags for the full annulus +uint_t flagInterior = 0; +uint_t flagInnerBoundary = 1; +uint_t flagOuterBoundary = 2; + MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, const real_t phiLeft, @@ -77,17 +82,47 @@ MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, const WALBERLA_ASSERT_GREATER( nTan, 0 ); WALBERLA_ASSERT_GREATER( nRad, 0 ); - // mesh partial annulus in polar coordinates + // mesh a rectangle representing the annulus in polar coordinates MeshInfo meshInfo = MeshInfo::meshRectangle( Point2D( {rhoMin, real_c( 0.0 )} ), Point2D( {rhoMax, real_c( 2.0 ) * pi} ), flavour, nRad, nTan ); - // close domain at phi = 0 = 2*pi + // determine some tolerances for further operations const real_t tolFactor = 0.1; const real_t xTolerance = ( ( rhoMax - rhoMin ) / real_c( nRad ) ) * tolFactor; const real_t yTolerance = ( 2.0 * pi / real_c( nTan ) ) * tolFactor; - // Gather all vertices on the top boundary and matching - // vertices on the bottom boundary. + // set boundary flags for vertices; + // note that the interior of the left and right edge of the rectangle will be glued together + for ( auto& it : meshInfo.vertices_ ) + { + real_t rho = it.second.getCoordinates()[0]; + if ( std::abs( rhoMax - rho ) < xTolerance ) + { + it.second.setBoundaryFlag( flagOuterBoundary ); + } + else if ( std::abs( rhoMin - rho ) < xTolerance ) + { + it.second.setBoundaryFlag( flagInnerBoundary ); + } + else + { + it.second.setBoundaryFlag( flagInterior ); + } + } + + // now do the same for edges + meshInfo.deduceEdgeFlagsFromVertices( flagInterior ); + + // and (re)set flags for all faces + for ( auto& it : meshInfo.faces_ ) + { + it.second.setBoundaryFlag( flagInterior ); + } + + // close domain at phi = 0 = 2*pi: + + // Gather all vertices on the top boundary (phi = 2*pi) and matching + // vertices on the bottom boundary (phi=0). std::map< IDType, Vertex > bottomBoundaryVertices; // map from top ID to bottom vertex instance for ( const auto& it : meshInfo.vertices_ ) { @@ -136,61 +171,60 @@ MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, const } // Replace vertex IDs in all primitives. - std::vector< std::array< IDType, 2 > > edgesToRemove; + std::vector< std::array< IDType, 2 > > edgesToRemove; std::map< std::array< IDType, 2 >, Edge > newEdges; - for ( const auto & it : meshInfo.edges_ ) + for ( const auto& it : meshInfo.edges_ ) { auto oldEdgeID = it.first; - std::array< IDType, 2 > newKey = oldEdgeID; - bool replaceEdge = false; + std::array< IDType, 2 > newKey = oldEdgeID; + bool replaceEdge = false; if ( bottomBoundaryVertices.count( oldEdgeID[0] ) > 0 ) { - newKey[0] = bottomBoundaryVertices[oldEdgeID[0]].getID(); + newKey[0] = bottomBoundaryVertices[oldEdgeID[0]].getID(); replaceEdge = true; - } if ( bottomBoundaryVertices.count( oldEdgeID[1] ) > 0 ) { - newKey[1] = bottomBoundaryVertices[oldEdgeID[1]].getID(); + newKey[1] = bottomBoundaryVertices[oldEdgeID[1]].getID(); replaceEdge = true; } if ( replaceEdge ) { edgesToRemove.push_back( it.first ); - Edge newEdge( newKey, 0 ); + Edge newEdge( newKey, flagInterior ); newEdges[newKey] = newEdge; } } - for ( const auto & it : edgesToRemove ) + for ( const auto& it : edgesToRemove ) { meshInfo.edges_.erase( it ); } - for ( const auto & it : newEdges ) + for ( const auto& it : newEdges ) { meshInfo.edges_[it.first] = it.second; } - std::vector< std::vector< IDType > > facesToRemove; + std::vector< std::vector< IDType > > facesToRemove; std::map< std::vector< IDType >, Face > newFaces; - for ( const auto & it : meshInfo.faces_ ) + for ( const auto& it : meshInfo.faces_ ) { auto oldFaceID = it.first; - std::vector< IDType > newKey = oldFaceID; - bool replaceFace = false; + std::vector< IDType > newKey = oldFaceID; + bool replaceFace = false; for ( uint_t i = 0; i < 3; i++ ) { if ( bottomBoundaryVertices.count( oldFaceID[i] ) > 0 ) { - newKey[i] = bottomBoundaryVertices[oldFaceID[i]].getID(); + newKey[i] = bottomBoundaryVertices[oldFaceID[i]].getID(); replaceFace = true; } } @@ -198,22 +232,21 @@ MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, const if ( replaceFace ) { facesToRemove.push_back( it.first ); - Face newFace( newKey, 0 ); + Face newFace( newKey, flagInterior ); newFaces[newKey] = newFace; } } - for ( const auto & it : facesToRemove ) + for ( const auto& it : facesToRemove ) { meshInfo.faces_.erase( it ); } - for ( const auto & it : newFaces ) + for ( const auto& it : newFaces ) { meshInfo.faces_[it.first] = it.second; } - // map vertex coordinates to cartesian domain Point3D node; node[2] = 0.0; @@ -257,26 +290,26 @@ MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, uint_t for ( uint_t i = 0; i < nTan; ++i ) { // inner boundary - node[0] = rhoMin * std::cos( (real_t) i * deltaPhi ); - node[1] = rhoMin * std::sin( (real_t) i * deltaPhi ); - meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), 1 ); + node[0] = rhoMin * std::cos( real_c( i ) * deltaPhi ); + node[1] = rhoMin * std::sin( real_c( i ) * deltaPhi ); + meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), flagInnerBoundary ); indices2id.insert( {{i, 0}, id} ); ++id; // interior nodes for ( uint_t j = 1; j < nRad; ++j ) { - node[0] = ( rhoMin + (real_t) j * deltaRho ) * std::cos( (real_t) i * deltaPhi ); - node[1] = ( rhoMin + (real_t) j * deltaRho ) * std::sin( (real_t) i * deltaPhi ); - meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), 0 ); + node[0] = ( rhoMin + real_c( j ) * deltaRho ) * std::cos( real_c( i ) * deltaPhi ); + node[1] = ( rhoMin + real_c( j ) * deltaRho ) * std::sin( real_c( i ) * deltaPhi ); + meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), flagInterior ); indices2id.insert( {{i, j}, id} ); ++id; } // outer boundary - node[0] = rhoMax * std::cos( (real_t) i * deltaPhi ); - node[1] = rhoMax * std::sin( (real_t) i * deltaPhi ); - meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), 1 ); + node[0] = rhoMax * std::cos( real_c( i ) * deltaPhi ); + node[1] = rhoMax * std::sin( real_c( i ) * deltaPhi ); + meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), flagOuterBoundary ); indices2id.insert( {{i, nRad}, id} ); ++id; } @@ -291,93 +324,40 @@ MeshInfo MeshInfo::meshAnnulus( const real_t rhoMin, const real_t rhoMax, uint_t // -------------------- // generate triangles // -------------------- - real_t midPhi = 0.0; - real_t midRho = 0.0; + real_t midPhi = real_c( 0 ); + real_t midRho = real_c( 0 ); for ( uint_t i = 0; i < nTan; ++i ) { for ( uint_t j = 0; j < nRad; ++j ) { // add new central vertex - midPhi = ( (real_t) i + 0.5 ) * deltaPhi; - midRho = rhoMin + ( (real_t) j + 0.5 ) * deltaRho; + midPhi = ( real_c( i ) + real_c( 0.5 ) ) * deltaPhi; + midRho = rhoMin + ( real_c( j ) + real_c( 0.5 ) ) * deltaRho; node[0] = midRho * std::cos( midPhi ); node[1] = midRho * std::sin( midPhi ); - meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), 0 ); + meshInfo.vertices_[id] = MeshInfo::Vertex( id, Point3D( node ), flagInterior ); // add four sub-triangles of cell - meshInfo.addFace( Face( {getIDX( ( i ), j ), id, getIDX( ( i + 1 ) % nTan, j )}, 0 ) ); - meshInfo.addFace( Face( {getIDX( ( i + 1 ) % nTan, j ), id, getIDX( ( i + 1 ) % nTan, j + 1 )}, 0 ) ); - meshInfo.addFace( Face( {getIDX( ( i + 1 ) % nTan, j + 1 ), id, getIDX( ( i ), j + 1 )}, 0 ) ); - meshInfo.addFace( Face( {getIDX( ( i ), j + 1 ), id, getIDX( ( i ), j )}, 0 ) ); + meshInfo.addFace( Face( {getIDX( ( i ), j ), id, getIDX( ( i + 1 ) % nTan, j )}, flagInterior ) ); + meshInfo.addFace( Face( {getIDX( ( i + 1 ) % nTan, j ), id, getIDX( ( i + 1 ) % nTan, j + 1 )}, flagInterior ) ); + meshInfo.addFace( Face( {getIDX( ( i + 1 ) % nTan, j + 1 ), id, getIDX( ( i ), j + 1 )}, flagInterior ) ); + meshInfo.addFace( Face( {getIDX( ( i ), j + 1 ), id, getIDX( ( i ), j )}, flagInterior ) ); ++id; } } // generate edges from faces - meshInfo.deriveEdgesForFullAnnulus( rhoMin + 0.1 * deltaRho, rhoMax - 0.1 * deltaRho ); - - // done - return meshInfo; -} - -void MeshInfo::deriveEdgesForFullAnnulus( real_t minTol, real_t maxTol ) -{ - MeshInfo::FaceContainer faces = this->getFaces(); - MeshInfo::VertexContainer verts = this->getVertices(); - - uint_t edgeBoundaryFlag = 0; - Point3D node1, node2; - real_t radius1, radius2; - - for ( const auto& it : faces ) + for ( const auto& it : meshInfo.faces_ ) { - // extract the three nodes of the face std::vector< IDType > fNode = it.second.getVertices(); - - // determine their position w.r.t. the boundary - std::vector< uint_t > meshBoundaryFlags( 3 ); - meshBoundaryFlags[0] = verts.find( fNode[0] )->second.getBoundaryFlag(); - meshBoundaryFlags[1] = verts.find( fNode[1] )->second.getBoundaryFlag(); - meshBoundaryFlags[2] = verts.find( fNode[2] )->second.getBoundaryFlag(); - - // set the three edges of triangle, edge is on boundary, if both - // its vertices are - edgeBoundaryFlag = 0; - if ( meshBoundaryFlags[0] == 1 && meshBoundaryFlags[1] == 1 ) - { - radius1 = verts.find( fNode[0] )->second.getCoordinates().norm(); - radius2 = verts.find( fNode[1] )->second.getCoordinates().norm(); - if ( ( radius1 < minTol && radius2 < minTol ) || ( radius1 > maxTol && radius2 < maxTol ) ) - { - edgeBoundaryFlag = 1; - } - } - this->addEdge( Edge( {fNode[0], fNode[1]}, edgeBoundaryFlag ) ); - - edgeBoundaryFlag = 0; - if ( meshBoundaryFlags[0] == 1 && meshBoundaryFlags[2] == 1 ) - { - radius1 = verts.find( fNode[0] )->second.getCoordinates().norm(); - radius2 = verts.find( fNode[2] )->second.getCoordinates().norm(); - if ( ( radius1 < minTol && radius2 < minTol ) || ( radius1 > maxTol && radius2 < maxTol ) ) - { - edgeBoundaryFlag = 1; - } - } - this->addEdge( Edge( {fNode[0], fNode[2]}, edgeBoundaryFlag ) ); - - edgeBoundaryFlag = 0; - if ( meshBoundaryFlags[1] == 1 && meshBoundaryFlags[2] == 1 ) - { - radius1 = verts.find( fNode[1] )->second.getCoordinates().norm(); - radius2 = verts.find( fNode[2] )->second.getCoordinates().norm(); - if ( ( radius1 < minTol && radius2 < minTol ) || ( radius1 > maxTol && radius2 < maxTol ) ) - { - edgeBoundaryFlag = 1; - } - } - this->addEdge( Edge( {fNode[1], fNode[2]}, edgeBoundaryFlag ) ); + meshInfo.addEdge( Edge( {fNode[0], fNode[1]}, flagInterior ) ); + meshInfo.addEdge( Edge( {fNode[0], fNode[2]}, flagInterior ) ); + meshInfo.addEdge( Edge( {fNode[1], fNode[2]}, flagInterior ) ); } + meshInfo.deduceEdgeFlagsFromVertices( flagInterior ); + + // done + return meshInfo; } } // namespace hyteg diff --git a/src/hyteg/mesh/MeshInfo.cpp b/src/hyteg/mesh/MeshInfo.cpp index 87aaa4a6b78ead07a1ed5b011c8210657e709e21..000c445a98cdca15b71bc067f4a7cfb32b0713e5 100644 --- a/src/hyteg/mesh/MeshInfo.cpp +++ b/src/hyteg/mesh/MeshInfo.cpp @@ -20,151 +20,100 @@ #include "hyteg/mesh/MeshInfo.hpp" -#include "core/logging/Logging.h" -#include "core/debug/CheckFunctions.h" -#include "core/debug/Debug.h" - #include #include +#include "core/debug/CheckFunctions.h" +#include "core/debug/Debug.h" +#include "core/logging/Logging.h" + namespace hyteg { -void MeshInfo::setAllMeshBoundaryFlags( const uint_t& meshBoundaryFlag ) +void MeshInfo::addVertex( const Vertex& vertex ) { - for ( auto& p : vertices_ ) - { - p.second.setBoundaryFlag( meshBoundaryFlag ); - } + vertices_[vertex.getID()] = vertex; +} + +void MeshInfo::addEdge( const Edge& edge ) +{ + WALBERLA_CHECK_UNEQUAL( edge.getVertices()[0], edge.getVertices()[1], "[Mesh] Mesh contains edge with zero length." ); - for ( auto& p : edges_ ) + std::array< IDType, 2 > sortedVertexIDs; + + if ( edge.getVertices()[0] < edge.getVertices()[1] ) { - p.second.setBoundaryFlag( meshBoundaryFlag ); + sortedVertexIDs = edge.getVertices(); } - - for ( auto& p : faces_ ) + else if ( edge.getVertices()[1] < edge.getVertices()[0] ) { - p.second.setBoundaryFlag( meshBoundaryFlag ); + sortedVertexIDs[0] = edge.getVertices()[1]; + sortedVertexIDs[1] = edge.getVertices()[0]; } - for ( auto& p : cells_ ) + if ( edges_.count( sortedVertexIDs ) == 0 ) { - p.second.setBoundaryFlag( meshBoundaryFlag ); + edges_[sortedVertexIDs] = Edge( sortedVertexIDs, edge.getBoundaryFlag() ); } } -void MeshInfo::setMeshBoundaryFlagsByVertexLocation( const uint_t& meshBoundaryFlag, - const std::function< bool( const Point3D& x ) >& onBoundary, - const bool& allVertices ) +void MeshInfo::addFace( const Face& face ) { - auto cond = [allVertices, onBoundary]( const std::vector< Point3D >& coordinates ) { - if ( allVertices ) - return std::all_of( coordinates.begin(), coordinates.end(), onBoundary ); - else - return std::any_of( coordinates.begin(), coordinates.end(), onBoundary ); - }; - - for ( auto& p : vertices_ ) - { - if ( cond( {p.second.getCoordinates()} ) ) - p.second.setBoundaryFlag( meshBoundaryFlag ); - } + std::vector< IDType > sortedVertexIDs = face.getVertices(); + std::sort( sortedVertexIDs.begin(), sortedVertexIDs.end() ); - for ( auto& p : edges_ ) - { - std::vector< Point3D > coordinates; - for ( auto v : p.second.getVertices() ) - { - coordinates.push_back( vertices_[v].getCoordinates() ); - } - if ( cond( coordinates ) ) - p.second.setBoundaryFlag( meshBoundaryFlag ); - } - - for ( auto& p : faces_ ) - { - std::vector< Point3D > coordinates; - for ( auto v : p.second.getVertices() ) - { - coordinates.push_back( vertices_[v].getCoordinates() ); - } - if ( cond( coordinates ) ) - p.second.setBoundaryFlag( meshBoundaryFlag ); - } + WALBERLA_CHECK_EQUAL( std::set< IDType >( sortedVertexIDs.begin(), sortedVertexIDs.end() ).size(), + sortedVertexIDs.size(), + "[Mesh] Mesh contains face with duplicate vertices." ); - for ( auto& p : cells_ ) + if ( faces_.count( sortedVertexIDs ) == 0 ) { - std::vector< Point3D > coordinates; - for ( auto v : p.second.getVertices() ) - { - coordinates.push_back( vertices_[v].getCoordinates() ); - } - if ( cond( coordinates ) ) - p.second.setBoundaryFlag( meshBoundaryFlag ); + faces_[sortedVertexIDs] = Face( sortedVertexIDs, face.getBoundaryFlag() ); } } -void MeshInfo::addVertex( const Vertex & vertex ) -{ - vertices_[vertex.getID()] = vertex; -} - -void MeshInfo::addEdge( const Edge & edge ) +void MeshInfo::addCellAndAllEdgesAndFaces( const Cell& cell ) { - WALBERLA_CHECK_UNEQUAL( edge.getVertices()[0], edge.getVertices()[1], "[Mesh] Mesh contains edge with zero length." ); - - std::array< IDType, 2 > sortedVertexIDs; - - if ( edge.getVertices()[0] < edge.getVertices()[1] ) - { - sortedVertexIDs = edge.getVertices(); - } - else if ( edge.getVertices()[1] < edge.getVertices()[0] ) - { - sortedVertexIDs[0] = edge.getVertices()[1]; - sortedVertexIDs[1] = edge.getVertices()[0]; - } - - if (edges_.count( sortedVertexIDs ) == 0) - { - edges_[ sortedVertexIDs ] = Edge( sortedVertexIDs, edge.getBoundaryFlag() ); - } + const auto cellCoordinates = cell.getVertices(); + + WALBERLA_ASSERT_EQUAL( cells_.count( cellCoordinates ), 0 ); + cells_[cellCoordinates] = cell; + + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[0], cellCoordinates[1]}} ), 0 ) ); + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[0], cellCoordinates[2]}} ), 0 ) ); + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[0], cellCoordinates[3]}} ), 0 ) ); + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[1], cellCoordinates[2]}} ), 0 ) ); + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[1], cellCoordinates[3]}} ), 0 ) ); + addEdge( Edge( std::array< IDType, 2 >( {{cellCoordinates[2], cellCoordinates[3]}} ), 0 ) ); + + addFace( Face( std::vector< IDType >( {{cellCoordinates[0], cellCoordinates[1], cellCoordinates[2]}} ), 0 ) ); + addFace( Face( std::vector< IDType >( {{cellCoordinates[0], cellCoordinates[1], cellCoordinates[3]}} ), 0 ) ); + addFace( Face( std::vector< IDType >( {{cellCoordinates[0], cellCoordinates[2], cellCoordinates[3]}} ), 0 ) ); + addFace( Face( std::vector< IDType >( {{cellCoordinates[1], cellCoordinates[2], cellCoordinates[3]}} ), 0 ) ); } - -void MeshInfo::addFace( const Face & face ) +void MeshInfo::deduceEdgeFlagsFromVertices( uint_t flagInconsistent ) { - std::vector< IDType > sortedVertexIDs = face.getVertices(); - std::sort( sortedVertexIDs.begin(), sortedVertexIDs.end() ); - - WALBERLA_CHECK_EQUAL( std::set< IDType >( sortedVertexIDs.begin(), sortedVertexIDs.end() ).size(), sortedVertexIDs.size(), - "[Mesh] Mesh contains face with duplicate vertices." ); - - if ( faces_.count( sortedVertexIDs ) == 0 ) - { - faces_[ sortedVertexIDs ] = Face( sortedVertexIDs, face.getBoundaryFlag() ); - } - + for ( auto& it : edges_ ) + { + Edge& edge = it.second; + auto vtxID = edge.getVertices(); + auto f1 = vertices_[vtxID[0]].getBoundaryFlag(); + auto f2 = vertices_[vtxID[1]].getBoundaryFlag(); + edge.setBoundaryFlag( f1 == f2 ? f1 : flagInconsistent ); + } } - -void MeshInfo::addCellAndAllEdgesAndFaces( const Cell & cell ) +void MeshInfo::deduceFaceFlagsFromVertices( uint_t flagInconsistent ) { - const auto cellCoordinates = cell.getVertices(); - - WALBERLA_ASSERT_EQUAL( cells_.count( cellCoordinates ), 0 ); - cells_[cellCoordinates] = cell; - - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[0], cellCoordinates[1] }} ), 0 ) ); - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[0], cellCoordinates[2] }} ), 0 ) ); - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[0], cellCoordinates[3] }} ), 0 ) ); - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[1], cellCoordinates[2] }} ), 0 ) ); - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[1], cellCoordinates[3] }} ), 0 ) ); - addEdge( Edge( std::array< IDType, 2 >( {{ cellCoordinates[2], cellCoordinates[3] }} ), 0 ) ); - - addFace( Face( std::vector< IDType >( {{ cellCoordinates[0], cellCoordinates[1], cellCoordinates[2] }} ), 0 ) ); - addFace( Face( std::vector< IDType >( {{ cellCoordinates[0], cellCoordinates[1], cellCoordinates[3] }} ), 0 ) ); - addFace( Face( std::vector< IDType >( {{ cellCoordinates[0], cellCoordinates[2], cellCoordinates[3] }} ), 0 ) ); - addFace( Face( std::vector< IDType >( {{ cellCoordinates[1], cellCoordinates[2], cellCoordinates[3] }} ), 0 ) ); + for ( auto& it : faces_ ) + { + Face& face = it.second; + auto vtxID = face.getVertices(); + auto f1 = vertices_[vtxID[0]].getBoundaryFlag(); + auto f2 = vertices_[vtxID[1]].getBoundaryFlag(); + auto f3 = vertices_[vtxID[2]].getBoundaryFlag(); + face.setBoundaryFlag( f1 == f2 && f1 == f3 ? f1 : flagInconsistent ); + } } } // namespace hyteg diff --git a/src/hyteg/mesh/MeshInfo.hpp b/src/hyteg/mesh/MeshInfo.hpp index 8fee9e1a71729f6c8e3630895a08f039ae46f5fe..ba1fc77373b1cdcb340953a025ede54885a402b0 100644 --- a/src/hyteg/mesh/MeshInfo.hpp +++ b/src/hyteg/mesh/MeshInfo.hpp @@ -383,13 +383,6 @@ class MeshInfo } } - void setAllMeshBoundaryFlags( const uint_t& meshBoundaryFlag ); - - /// Every primitive for which onBoundary() returns true for all / any (if allVertices == true / false) of the primitve's vertices is assigned the passed mesh boundary flag. - void setMeshBoundaryFlagsByVertexLocation( const uint_t& meshBoundaryFlag, - const std::function< bool( const Point3D& x ) >& onBoundary, - const bool& allVertices = true ); - typedef std::map< IDType, Vertex > VertexContainer; typedef std::map< std::array< IDType, 2 >, Edge > EdgeContainer; typedef std::map< std::vector< IDType >, Face > FaceContainer; @@ -595,14 +588,19 @@ class MeshInfo /// \param tol parameter used to determine, whether an edge is part of the boundary, or not void deriveEdgesForRectangles( const Point2D lowerLeft, const Point2D upperRight, real_t tol ); - /// Derive information on edges from vertices and faces (for full annulus) + /// Set the boundaryFlag of all edges automatically from that of their vertices - /// This method is used in the 2D inline mesh generator for the full annulus. The latter - /// provides only information on vertices and faces. The information on the edges is then - /// derived from the faces, as each face has three edges. - /// \param minTol parameter used to determine, whether an edge is part of the inner boundary, or not - /// \param maxTol parameter used to determine, whether an edge is part of the outer boundary, or not - void deriveEdgesForFullAnnulus( real_t minTol, real_t maxTol ); + /// Calling the method results in the boundaryFlag of all edges in the MeshInfo being (re)set. + /// The flag of a single edge is set to that of its two vertices, if the latter agree. Otherwise + /// it will be set to the value of the undecided argument. + void deduceEdgeFlagsFromVertices( uint_t flagInconsistent = 0 ); + + /// Set the boundaryFlag of all faces automatically from that of their vertices + + /// Calling the method results in the boundaryFlag of all faces in the MeshInfo being (re)set. + /// The flag of a single face is set to that of its three vertices, if the latter agree. Otherwise + /// it will be set to the value of the undecided argument. + void deduceFaceFlagsFromVertices( uint_t flagInconsistent = 0 ); VertexContainer vertices_; EdgeContainer edges_; diff --git a/src/hyteg/primitivestorage/SetupPrimitiveStorage.cpp b/src/hyteg/primitivestorage/SetupPrimitiveStorage.cpp index 438be39fa8ee9c577d4c065c590a1a96beb624ec..5292009e62905684bfe046ba45ca4dad4d00ccf0 100644 --- a/src/hyteg/primitivestorage/SetupPrimitiveStorage.cpp +++ b/src/hyteg/primitivestorage/SetupPrimitiveStorage.cpp @@ -18,378 +18,403 @@ * along with this program. If not, see . */ +#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp" + +#include +#include +#include + #include "core/debug/CheckFunctions.h" #include "core/debug/Debug.h" #include "core/logging/Logging.h" -#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp" #include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp" -#include -#include -#include +using walberla::real_c; namespace hyteg { -SetupPrimitiveStorage::SetupPrimitiveStorage( const MeshInfo & meshInfo, const uint_t & numberOfProcesses ) : - numberOfProcesses_( numberOfProcesses ) +SetupPrimitiveStorage::SetupPrimitiveStorage( const MeshInfo& meshInfo, const uint_t& numberOfProcesses ) +: numberOfProcesses_( numberOfProcesses ) { - WALBERLA_ASSERT_GREATER( numberOfProcesses_, 0, "Number of processes must be positive" ); - - // since the MeshInfo IDs of the vertices do not necessarily - // match the primitive IDs of the vertices in the SetupStorage, we need an assignment - std::map< uint_t, PrimitiveID > meshVertexIDToPrimitiveID; - - // We cache the inserted primitives (edges, faces and cells) by filling - // these maps with the surrounding vertexIDs as keys and the inserted - // PrimitiveIDs as values. - // This way we do not need to search for the neighboring lower level - // primitives when building inner primitives. - std::map< std::vector< PrimitiveID >, PrimitiveID > vertexIDsToEdgeIDs; - std::map< std::vector< PrimitiveID >, PrimitiveID > vertexIDsToFaceIDs; - auto findCachedPrimitiveID = []( const std::vector< PrimitiveID > & unsortedPrimitiveIDs, const std::map< std::vector< PrimitiveID >, PrimitiveID > & cache ) -> PrimitiveID - { - std::vector< PrimitiveID > sortedKey( unsortedPrimitiveIDs ); - std::sort( sortedKey.begin(), sortedKey.end() ); - WALBERLA_ASSERT_GREATER( cache.count( sortedKey ), 0, "Could not find primitive in cache during SetupStorage construction." ); - return cache.at( sortedKey ); - }; - - // Adding vertices to storage - const MeshInfo::VertexContainer vertices = meshInfo.getVertices(); - for ( const auto & it : vertices ) - { - const MeshInfo::Vertex meshInfoVertex = it.second; - - PrimitiveID vertexID = generatePrimitiveID(); - - meshVertexIDToPrimitiveID[ meshInfoVertex.getID() ] = vertexID; - - Point3D coordinates( meshInfoVertex.getCoordinates() ); - vertices_[ vertexID.getID() ] = std::make_shared< Vertex >( vertexID, coordinates ); - - setMeshBoundaryFlag( vertexID.getID(), meshInfoVertex.getBoundaryFlag() ); - } - - // Adding edges to storage - const MeshInfo::EdgeContainer edges = meshInfo.getEdges(); - for ( const auto & it : edges ) - { - const MeshInfo::Edge meshInfoEdge = it.second; - - PrimitiveID edgeID = generatePrimitiveID(); - - WALBERLA_ASSERT_EQUAL( meshInfoEdge.getVertices().size(), 2, "Edges are expected to have two vertices." ); - PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[ meshInfoEdge.getVertices().at( 0 ) ]; - PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[ meshInfoEdge.getVertices().at( 1 ) ]; - - std::array coords; - - coords[0] = vertices_[ vertexID0.getID() ]->getCoordinates(); - coords[1] = vertices_[ vertexID1.getID() ]->getCoordinates(); - - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID.getID() ), 0 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); - edges_[ edgeID.getID() ] = std::make_shared< Edge >( edgeID, vertexID0, vertexID1, coords); - - setMeshBoundaryFlag( edgeID.getID(), meshInfoEdge.getBoundaryFlag() ); - - // Adding edge ID as neighbor to SetupVertices - vertices_[ vertexID0.getID() ]->addEdge( edgeID ); - vertices_[ vertexID1.getID() ]->addEdge( edgeID ); - - // Caching neighboring vertices - std::vector< PrimitiveID > vertexIDs = {{ vertexID0, vertexID1 }}; - std::sort( vertexIDs.begin(), vertexIDs.end() ); - vertexIDsToEdgeIDs[ vertexIDs ] = edgeID; - } - - // Adding faces to storage - const MeshInfo::FaceContainer faces = meshInfo.getFaces(); - for ( const auto & it : faces ) - { - const MeshInfo::Face meshInfoFace = it.second; - - PrimitiveID faceID = generatePrimitiveID(); - - WALBERLA_ASSERT_EQUAL( meshInfoFace.getVertices().size(), 3, "Only supporting triangle faces." ); - PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[ meshInfoFace.getVertices().at( 0 ) ]; - PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[ meshInfoFace.getVertices().at( 1 ) ]; - PrimitiveID vertexID2 = meshVertexIDToPrimitiveID[ meshInfoFace.getVertices().at( 2 ) ]; - - WALBERLA_ASSERT_EQUAL( faces_.count( faceID.getID() ), 0 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID2.getID() ), 1 ); - - PrimitiveID edgeID0 = findCachedPrimitiveID( {{ vertexID0, vertexID1 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID1 = findCachedPrimitiveID( {{ vertexID0, vertexID2 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID2 = findCachedPrimitiveID( {{ vertexID1, vertexID2 }}, vertexIDsToEdgeIDs ); - - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID1.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID2.getID() ), 1 ); - - // Edge Orientation - std::array< int, 3 > edgeOrientation; - - PrimitiveID edge0Vertex0 = edges_[ edgeID0.getID() ]->getVertexID0(); - PrimitiveID edge0Vertex1 = edges_[ edgeID0.getID() ]->getVertexID1(); - PrimitiveID edge1Vertex0 = edges_[ edgeID1.getID() ]->getVertexID0(); - PrimitiveID edge1Vertex1 = edges_[ edgeID1.getID() ]->getVertexID1(); - PrimitiveID edge2Vertex0 = edges_[ edgeID2.getID() ]->getVertexID0(); - PrimitiveID edge2Vertex1 = edges_[ edgeID2.getID() ]->getVertexID1(); - - WALBERLA_UNUSED( edge2Vertex1 ); - - if ( edge0Vertex0 == vertexID0 ) - { - WALBERLA_ASSERT( edge0Vertex1 == vertexID1 ); - edgeOrientation[0] = 1; - } - else - { - WALBERLA_ASSERT( edge0Vertex0 == vertexID1 ); - WALBERLA_ASSERT( edge0Vertex1 == vertexID0 ); - edgeOrientation[0] = -1; - } - - if ( edge1Vertex0 == vertexID0 ) - { - WALBERLA_ASSERT( edge1Vertex1 == vertexID2 ); - edgeOrientation[1] = 1; - } - else - { - WALBERLA_ASSERT( edge1Vertex0 == vertexID2 ); - WALBERLA_ASSERT( edge1Vertex1 == vertexID0 ); - edgeOrientation[1] = -1; - } - - if ( edge2Vertex0 == vertexID1 ) - { - WALBERLA_ASSERT( edge2Vertex1 == vertexID2 ); - edgeOrientation[2] = 1; - } - else - { - WALBERLA_ASSERT( edge2Vertex0 == vertexID2 ); - WALBERLA_ASSERT( edge2Vertex1 == vertexID1 ); - edgeOrientation[2] = -1; - } - - // Corner coordinates - std::array< Point3D, 3 > coordinates; - std::array< PrimitiveID, 3 > vertexIDs; - std::vector< PrimitiveID > verticesOnBoundary; - std::vector< PrimitiveID > edgesOnBoundary; - - if (edgeOrientation[0] == 1) - { - coordinates[0] = vertices_[ edge0Vertex0.getID() ]->getCoordinates(); - coordinates[1] = vertices_[ edge0Vertex1.getID() ]->getCoordinates(); - - vertexIDs[0] = edge0Vertex0.getID(); - vertexIDs[1] = edge0Vertex1.getID(); - } - else - { - coordinates[0] = vertices_[ edge0Vertex1.getID() ]->getCoordinates(); - coordinates[1] = vertices_[ edge0Vertex0.getID() ]->getCoordinates(); - - vertexIDs[0] = edge0Vertex1.getID(); - vertexIDs[1] = edge0Vertex0.getID(); - } - - if (edgeOrientation[1] == 1) - { - coordinates[2] = vertices_[ edge1Vertex1.getID() ]->getCoordinates(); - - vertexIDs[2] = edge1Vertex1.getID(); - } - else - { - coordinates[2] = vertices_[ edge1Vertex0.getID() ]->getCoordinates(); - - vertexIDs[2] = edge1Vertex0.getID(); - } - - faces_[ faceID.getID() ] = std::shared_ptr< Face >( new Face( faceID, vertexIDs, {{edgeID0, edgeID1, edgeID2}}, edgeOrientation, coordinates ) ); - - setMeshBoundaryFlag( faceID.getID(), meshInfoFace.getBoundaryFlag() ); - - // Adding face ID to vertices as neighbors - vertices_[vertexIDs[0].getID()]->addFace(faceID); - vertices_[vertexIDs[1].getID()]->addFace(faceID); - vertices_[vertexIDs[2].getID()]->addFace(faceID); - - // Adding face ID to edges as neighbors - edges_[ edgeID0.getID() ]->addFace( faceID ); - edges_[ edgeID1.getID() ]->addFace( faceID ); - edges_[ edgeID2.getID() ]->addFace( faceID ); - - // Caching neighboring vertices - std::vector< PrimitiveID > neighboringVertexIDs = {{ vertexID0, vertexID1, vertexID2 }}; - std::sort( neighboringVertexIDs.begin(), neighboringVertexIDs.end() ); - vertexIDsToFaceIDs[ neighboringVertexIDs ] = faceID; - } - - // Adding cells to storage - const MeshInfo::CellContainer cells = meshInfo.getCells(); - for ( const auto & it : cells ) - { - const MeshInfo::Cell meshInfoCell = it.second; - - PrimitiveID cellID = generatePrimitiveID(); - - WALBERLA_ASSERT_EQUAL( meshInfoCell.getVertices().size(), 4, "Only supporting tetrahedron cells." ); - - PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[ meshInfoCell.getVertices().at( 0 ) ]; - PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[ meshInfoCell.getVertices().at( 1 ) ]; - PrimitiveID vertexID2 = meshVertexIDToPrimitiveID[ meshInfoCell.getVertices().at( 2 ) ]; - PrimitiveID vertexID3 = meshVertexIDToPrimitiveID[ meshInfoCell.getVertices().at( 3 ) ]; - - PrimitiveID edgeID0 = findCachedPrimitiveID( {{ vertexID0, vertexID1 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID1 = findCachedPrimitiveID( {{ vertexID0, vertexID2 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID2 = findCachedPrimitiveID( {{ vertexID1, vertexID2 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID3 = findCachedPrimitiveID( {{ vertexID0, vertexID3 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID4 = findCachedPrimitiveID( {{ vertexID1, vertexID3 }}, vertexIDsToEdgeIDs ); - PrimitiveID edgeID5 = findCachedPrimitiveID( {{ vertexID2, vertexID3 }}, vertexIDsToEdgeIDs ); - - PrimitiveID faceID0 = findCachedPrimitiveID( {{ vertexID0, vertexID1, vertexID2 }}, vertexIDsToFaceIDs ); - PrimitiveID faceID1 = findCachedPrimitiveID( {{ vertexID0, vertexID1, vertexID3 }}, vertexIDsToFaceIDs ); - PrimitiveID faceID2 = findCachedPrimitiveID( {{ vertexID0, vertexID2, vertexID3 }}, vertexIDsToFaceIDs ); - PrimitiveID faceID3 = findCachedPrimitiveID( {{ vertexID1, vertexID2, vertexID3 }}, vertexIDsToFaceIDs ); - - WALBERLA_ASSERT_EQUAL( cells_.count( cellID.getID() ), 0 ); - - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID2.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID3.getID() ), 1 ); - - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID1.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID2.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID3.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID4.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( edges_.count( edgeID5.getID() ), 1 ); - - WALBERLA_ASSERT_EQUAL( faces_.count( faceID0.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( faces_.count( faceID1.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( faces_.count( faceID2.getID() ), 1 ); - WALBERLA_ASSERT_EQUAL( faces_.count( faceID3.getID() ), 1 ); - - std::vector< PrimitiveID > cellVertices = {{ vertexID0, vertexID1, vertexID2, vertexID3 }}; - std::vector< PrimitiveID > cellEdges = {{ edgeID0, edgeID1, edgeID2, edgeID3, edgeID4, edgeID5 }}; - std::vector< PrimitiveID > cellFaces = {{ faceID0, faceID1, faceID2, faceID3 }}; - - for ( const auto & id : cellVertices ) { WALBERLA_ASSERT( vertexExists( id ) ); vertices_[ id.getID() ]->addCell( cellID ); } - for ( const auto & id : cellEdges ) { WALBERLA_ASSERT( edgeExists( id ) ); edges_[ id.getID() ]->addCell( cellID ); } - for ( const auto & id : cellFaces ) { WALBERLA_ASSERT( faceExists( id ) ); faces_[ id.getID() ]->addCell( cellID ); } - - std::array< Point3D, 4 > cellCoordinates = {{ vertices_.at( vertexID0.getID() )->getCoordinates(), - vertices_.at( vertexID1.getID() )->getCoordinates(), - vertices_.at( vertexID2.getID() )->getCoordinates(), - vertices_.at( vertexID3.getID() )->getCoordinates() }}; - - std::array< std::map< uint_t, uint_t >, 6 > edgeLocalVertexToCellLocalVertexMaps; - - // edgeLocalVertexToCellLocalVertexMaps[ cellLocalEdgeID ][ edgeLocalVertexID ] = cellLocalVertexID; - - edgeLocalVertexToCellLocalVertexMaps[0][ edges_.at( edgeID0.getID() )->vertex_index( vertexID0 ) ] = 0; - edgeLocalVertexToCellLocalVertexMaps[0][ edges_.at( edgeID0.getID() )->vertex_index( vertexID1 ) ] = 1; - - edgeLocalVertexToCellLocalVertexMaps[1][ edges_.at( edgeID1.getID() )->vertex_index( vertexID0 ) ] = 0; - edgeLocalVertexToCellLocalVertexMaps[1][ edges_.at( edgeID1.getID() )->vertex_index( vertexID2 ) ] = 2; - - edgeLocalVertexToCellLocalVertexMaps[2][ edges_.at( edgeID2.getID() )->vertex_index( vertexID1 ) ] = 1; - edgeLocalVertexToCellLocalVertexMaps[2][ edges_.at( edgeID2.getID() )->vertex_index( vertexID2 ) ] = 2; - - edgeLocalVertexToCellLocalVertexMaps[3][ edges_.at( edgeID3.getID() )->vertex_index( vertexID0 ) ] = 0; - edgeLocalVertexToCellLocalVertexMaps[3][ edges_.at( edgeID3.getID() )->vertex_index( vertexID3 ) ] = 3; - - edgeLocalVertexToCellLocalVertexMaps[4][ edges_.at( edgeID4.getID() )->vertex_index( vertexID1 ) ] = 1; - edgeLocalVertexToCellLocalVertexMaps[4][ edges_.at( edgeID4.getID() )->vertex_index( vertexID3 ) ] = 3; - - edgeLocalVertexToCellLocalVertexMaps[5][ edges_.at( edgeID5.getID() )->vertex_index( vertexID2 ) ] = 2; - edgeLocalVertexToCellLocalVertexMaps[5][ edges_.at( edgeID5.getID() )->vertex_index( vertexID3 ) ] = 3; - - std::array< std::map< uint_t, uint_t >, 4 > faceLocalVertexToCellLocalVertexMaps; - - // faceLocalVertexToCellLocalVertexMaps[ cellLocalFaceID ][ faceLocalVertexID ] = cellLocalVertexID; - - faceLocalVertexToCellLocalVertexMaps[0][ faces_.at( faceID0.getID() )->vertex_index( vertexID0 ) ] = 0; - faceLocalVertexToCellLocalVertexMaps[0][ faces_.at( faceID0.getID() )->vertex_index( vertexID1 ) ] = 1; - faceLocalVertexToCellLocalVertexMaps[0][ faces_.at( faceID0.getID() )->vertex_index( vertexID2 ) ] = 2; - - faceLocalVertexToCellLocalVertexMaps[1][ faces_.at( faceID1.getID() )->vertex_index( vertexID0 ) ] = 0; - faceLocalVertexToCellLocalVertexMaps[1][ faces_.at( faceID1.getID() )->vertex_index( vertexID1 ) ] = 1; - faceLocalVertexToCellLocalVertexMaps[1][ faces_.at( faceID1.getID() )->vertex_index( vertexID3 ) ] = 3; - - faceLocalVertexToCellLocalVertexMaps[2][ faces_.at( faceID2.getID() )->vertex_index( vertexID0 ) ] = 0; - faceLocalVertexToCellLocalVertexMaps[2][ faces_.at( faceID2.getID() )->vertex_index( vertexID2 ) ] = 2; - faceLocalVertexToCellLocalVertexMaps[2][ faces_.at( faceID2.getID() )->vertex_index( vertexID3 ) ] = 3; - - faceLocalVertexToCellLocalVertexMaps[3][ faces_.at( faceID3.getID() )->vertex_index( vertexID1 ) ] = 1; - faceLocalVertexToCellLocalVertexMaps[3][ faces_.at( faceID3.getID() )->vertex_index( vertexID2 ) ] = 2; - faceLocalVertexToCellLocalVertexMaps[3][ faces_.at( faceID3.getID() )->vertex_index( vertexID3 ) ] = 3; - - cells_[ cellID.getID() ] = std::make_shared< Cell >( cellID, cellVertices, cellEdges, cellFaces, cellCoordinates, edgeLocalVertexToCellLocalVertexMaps, faceLocalVertexToCellLocalVertexMaps ); - - setMeshBoundaryFlag( cellID.getID(), meshInfoCell.getBoundaryFlag() ); - } + WALBERLA_ASSERT_GREATER( numberOfProcesses_, 0, "Number of processes must be positive" ); + + // since the MeshInfo IDs of the vertices do not necessarily + // match the primitive IDs of the vertices in the SetupStorage, we need an assignment + std::map< uint_t, PrimitiveID > meshVertexIDToPrimitiveID; + + // We cache the inserted primitives (edges, faces and cells) by filling + // these maps with the surrounding vertexIDs as keys and the inserted + // PrimitiveIDs as values. + // This way we do not need to search for the neighboring lower level + // primitives when building inner primitives. + std::map< std::vector< PrimitiveID >, PrimitiveID > vertexIDsToEdgeIDs; + std::map< std::vector< PrimitiveID >, PrimitiveID > vertexIDsToFaceIDs; + auto findCachedPrimitiveID = []( const std::vector< PrimitiveID >& unsortedPrimitiveIDs, + const std::map< std::vector< PrimitiveID >, PrimitiveID >& cache ) -> PrimitiveID { + std::vector< PrimitiveID > sortedKey( unsortedPrimitiveIDs ); + std::sort( sortedKey.begin(), sortedKey.end() ); + WALBERLA_ASSERT_GREATER( + cache.count( sortedKey ), 0, "Could not find primitive in cache during SetupStorage construction." ); + return cache.at( sortedKey ); + }; + + // Adding vertices to storage + const MeshInfo::VertexContainer vertices = meshInfo.getVertices(); + for ( const auto& it : vertices ) + { + const MeshInfo::Vertex meshInfoVertex = it.second; + + PrimitiveID vertexID = generatePrimitiveID(); + + meshVertexIDToPrimitiveID[meshInfoVertex.getID()] = vertexID; + + Point3D coordinates( meshInfoVertex.getCoordinates() ); + vertices_[vertexID.getID()] = std::make_shared< Vertex >( vertexID, coordinates ); + + setMeshBoundaryFlag( vertexID.getID(), meshInfoVertex.getBoundaryFlag() ); + } + + // Adding edges to storage + const MeshInfo::EdgeContainer edges = meshInfo.getEdges(); + for ( const auto& it : edges ) + { + const MeshInfo::Edge meshInfoEdge = it.second; + + PrimitiveID edgeID = generatePrimitiveID(); + + WALBERLA_ASSERT_EQUAL( meshInfoEdge.getVertices().size(), 2, "Edges are expected to have two vertices." ); + PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[meshInfoEdge.getVertices().at( 0 )]; + PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[meshInfoEdge.getVertices().at( 1 )]; + + std::array< Point3D, 2 > coords; + + coords[0] = vertices_[vertexID0.getID()]->getCoordinates(); + coords[1] = vertices_[vertexID1.getID()]->getCoordinates(); + + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID.getID() ), 0 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); + edges_[edgeID.getID()] = std::make_shared< Edge >( edgeID, vertexID0, vertexID1, coords ); + + setMeshBoundaryFlag( edgeID.getID(), meshInfoEdge.getBoundaryFlag() ); + + // Adding edge ID as neighbor to SetupVertices + vertices_[vertexID0.getID()]->addEdge( edgeID ); + vertices_[vertexID1.getID()]->addEdge( edgeID ); + + // Caching neighboring vertices + std::vector< PrimitiveID > vertexIDs = {{vertexID0, vertexID1}}; + std::sort( vertexIDs.begin(), vertexIDs.end() ); + vertexIDsToEdgeIDs[vertexIDs] = edgeID; + } + + // Adding faces to storage + const MeshInfo::FaceContainer faces = meshInfo.getFaces(); + for ( const auto& it : faces ) + { + const MeshInfo::Face meshInfoFace = it.second; + + PrimitiveID faceID = generatePrimitiveID(); + + WALBERLA_ASSERT_EQUAL( meshInfoFace.getVertices().size(), 3, "Only supporting triangle faces." ); + PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[meshInfoFace.getVertices().at( 0 )]; + PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[meshInfoFace.getVertices().at( 1 )]; + PrimitiveID vertexID2 = meshVertexIDToPrimitiveID[meshInfoFace.getVertices().at( 2 )]; + + WALBERLA_ASSERT_EQUAL( faces_.count( faceID.getID() ), 0 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID2.getID() ), 1 ); + + PrimitiveID edgeID0 = findCachedPrimitiveID( {{vertexID0, vertexID1}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID1 = findCachedPrimitiveID( {{vertexID0, vertexID2}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID2 = findCachedPrimitiveID( {{vertexID1, vertexID2}}, vertexIDsToEdgeIDs ); + + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID1.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID2.getID() ), 1 ); + + // Edge Orientation + std::array< int, 3 > edgeOrientation; + + PrimitiveID edge0Vertex0 = edges_[edgeID0.getID()]->getVertexID0(); + PrimitiveID edge0Vertex1 = edges_[edgeID0.getID()]->getVertexID1(); + PrimitiveID edge1Vertex0 = edges_[edgeID1.getID()]->getVertexID0(); + PrimitiveID edge1Vertex1 = edges_[edgeID1.getID()]->getVertexID1(); + PrimitiveID edge2Vertex0 = edges_[edgeID2.getID()]->getVertexID0(); + PrimitiveID edge2Vertex1 = edges_[edgeID2.getID()]->getVertexID1(); + + WALBERLA_UNUSED( edge2Vertex1 ); + + if ( edge0Vertex0 == vertexID0 ) + { + WALBERLA_ASSERT( edge0Vertex1 == vertexID1 ); + edgeOrientation[0] = 1; + } + else + { + WALBERLA_ASSERT( edge0Vertex0 == vertexID1 ); + WALBERLA_ASSERT( edge0Vertex1 == vertexID0 ); + edgeOrientation[0] = -1; + } + + if ( edge1Vertex0 == vertexID0 ) + { + WALBERLA_ASSERT( edge1Vertex1 == vertexID2 ); + edgeOrientation[1] = 1; + } + else + { + WALBERLA_ASSERT( edge1Vertex0 == vertexID2 ); + WALBERLA_ASSERT( edge1Vertex1 == vertexID0 ); + edgeOrientation[1] = -1; + } + + if ( edge2Vertex0 == vertexID1 ) + { + WALBERLA_ASSERT( edge2Vertex1 == vertexID2 ); + edgeOrientation[2] = 1; + } + else + { + WALBERLA_ASSERT( edge2Vertex0 == vertexID2 ); + WALBERLA_ASSERT( edge2Vertex1 == vertexID1 ); + edgeOrientation[2] = -1; + } + + // Corner coordinates + std::array< Point3D, 3 > coordinates; + std::array< PrimitiveID, 3 > vertexIDs; + std::vector< PrimitiveID > verticesOnBoundary; + std::vector< PrimitiveID > edgesOnBoundary; - // add indirect neighbor faces - for ( auto& it : faces_ ) - { - auto faceID = it.first; - auto face = it.second; - - std::set< PrimitiveID > indirectNeighborsSet; - - for ( const auto& vertexID : face->neighborVertices() ) - { - auto vertex = getVertex( vertexID ); - for ( const auto& neighborFaceID : vertex->neighborFaces() ) - { - if ( neighborFaceID != faceID ) - { - indirectNeighborsSet.insert( neighborFaceID ); - } - } - } + if ( edgeOrientation[0] == 1 ) + { + coordinates[0] = vertices_[edge0Vertex0.getID()]->getCoordinates(); + coordinates[1] = vertices_[edge0Vertex1.getID()]->getCoordinates(); + + vertexIDs[0] = edge0Vertex0.getID(); + vertexIDs[1] = edge0Vertex1.getID(); + } + else + { + coordinates[0] = vertices_[edge0Vertex1.getID()]->getCoordinates(); + coordinates[1] = vertices_[edge0Vertex0.getID()]->getCoordinates(); + + vertexIDs[0] = edge0Vertex1.getID(); + vertexIDs[1] = edge0Vertex0.getID(); + } + + if ( edgeOrientation[1] == 1 ) + { + coordinates[2] = vertices_[edge1Vertex1.getID()]->getCoordinates(); + + vertexIDs[2] = edge1Vertex1.getID(); + } + else + { + coordinates[2] = vertices_[edge1Vertex0.getID()]->getCoordinates(); + + vertexIDs[2] = edge1Vertex0.getID(); + } + + faces_[faceID.getID()] = + std::shared_ptr< Face >( new Face( faceID, vertexIDs, {{edgeID0, edgeID1, edgeID2}}, edgeOrientation, coordinates ) ); + + setMeshBoundaryFlag( faceID.getID(), meshInfoFace.getBoundaryFlag() ); + + // Adding face ID to vertices as neighbors + vertices_[vertexIDs[0].getID()]->addFace( faceID ); + vertices_[vertexIDs[1].getID()]->addFace( faceID ); + vertices_[vertexIDs[2].getID()]->addFace( faceID ); + + // Adding face ID to edges as neighbors + edges_[edgeID0.getID()]->addFace( faceID ); + edges_[edgeID1.getID()]->addFace( faceID ); + edges_[edgeID2.getID()]->addFace( faceID ); + + // Caching neighboring vertices + std::vector< PrimitiveID > neighboringVertexIDs = {{vertexID0, vertexID1, vertexID2}}; + std::sort( neighboringVertexIDs.begin(), neighboringVertexIDs.end() ); + vertexIDsToFaceIDs[neighboringVertexIDs] = faceID; + } + + // Adding cells to storage + const MeshInfo::CellContainer cells = meshInfo.getCells(); + for ( const auto& it : cells ) + { + const MeshInfo::Cell meshInfoCell = it.second; + + PrimitiveID cellID = generatePrimitiveID(); + + WALBERLA_ASSERT_EQUAL( meshInfoCell.getVertices().size(), 4, "Only supporting tetrahedron cells." ); + + PrimitiveID vertexID0 = meshVertexIDToPrimitiveID[meshInfoCell.getVertices().at( 0 )]; + PrimitiveID vertexID1 = meshVertexIDToPrimitiveID[meshInfoCell.getVertices().at( 1 )]; + PrimitiveID vertexID2 = meshVertexIDToPrimitiveID[meshInfoCell.getVertices().at( 2 )]; + PrimitiveID vertexID3 = meshVertexIDToPrimitiveID[meshInfoCell.getVertices().at( 3 )]; + + PrimitiveID edgeID0 = findCachedPrimitiveID( {{vertexID0, vertexID1}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID1 = findCachedPrimitiveID( {{vertexID0, vertexID2}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID2 = findCachedPrimitiveID( {{vertexID1, vertexID2}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID3 = findCachedPrimitiveID( {{vertexID0, vertexID3}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID4 = findCachedPrimitiveID( {{vertexID1, vertexID3}}, vertexIDsToEdgeIDs ); + PrimitiveID edgeID5 = findCachedPrimitiveID( {{vertexID2, vertexID3}}, vertexIDsToEdgeIDs ); + + PrimitiveID faceID0 = findCachedPrimitiveID( {{vertexID0, vertexID1, vertexID2}}, vertexIDsToFaceIDs ); + PrimitiveID faceID1 = findCachedPrimitiveID( {{vertexID0, vertexID1, vertexID3}}, vertexIDsToFaceIDs ); + PrimitiveID faceID2 = findCachedPrimitiveID( {{vertexID0, vertexID2, vertexID3}}, vertexIDsToFaceIDs ); + PrimitiveID faceID3 = findCachedPrimitiveID( {{vertexID1, vertexID2, vertexID3}}, vertexIDsToFaceIDs ); + + WALBERLA_ASSERT_EQUAL( cells_.count( cellID.getID() ), 0 ); - face->indirectNeighborFaceIDs_.clear(); - face->indirectNeighborFaceIDs_.insert( face->indirectNeighborFaceIDs_.begin(), indirectNeighborsSet.begin(), indirectNeighborsSet.end() ); - } - - // add indirect neighbor cells - for ( auto& it : cells_ ) - { - auto cellID = it.first; - auto cell = it.second; - - std::set< PrimitiveID > indirectNeighborsSet; + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID1.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID2.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( vertices_.count( vertexID3.getID() ), 1 ); - for ( const auto& vertexID : cell->neighborVertices() ) - { - auto vertex = getVertex( vertexID ); - for ( const auto& neighborCellID : vertex->neighborCells() ) - { - if ( neighborCellID != cellID ) - { - indirectNeighborsSet.insert( neighborCellID ); - } - } - } + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID1.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID2.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID3.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID4.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( edges_.count( edgeID5.getID() ), 1 ); - cell->indirectNeighborCellIDs_.clear(); - cell->indirectNeighborCellIDs_.insert( cell->indirectNeighborCellIDs_.begin(), indirectNeighborsSet.begin(), indirectNeighborsSet.end() ); - } + WALBERLA_ASSERT_EQUAL( faces_.count( faceID0.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( faces_.count( faceID1.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( faces_.count( faceID2.getID() ), 1 ); + WALBERLA_ASSERT_EQUAL( faces_.count( faceID3.getID() ), 1 ); - loadbalancing::roundRobin( *this ); + std::vector< PrimitiveID > cellVertices = {{vertexID0, vertexID1, vertexID2, vertexID3}}; + std::vector< PrimitiveID > cellEdges = {{edgeID0, edgeID1, edgeID2, edgeID3, edgeID4, edgeID5}}; + std::vector< PrimitiveID > cellFaces = {{faceID0, faceID1, faceID2, faceID3}}; + + for ( const auto& id : cellVertices ) + { + WALBERLA_ASSERT( vertexExists( id ) ); + vertices_[id.getID()]->addCell( cellID ); + } + for ( const auto& id : cellEdges ) + { + WALBERLA_ASSERT( edgeExists( id ) ); + edges_[id.getID()]->addCell( cellID ); + } + for ( const auto& id : cellFaces ) + { + WALBERLA_ASSERT( faceExists( id ) ); + faces_[id.getID()]->addCell( cellID ); + } + + std::array< Point3D, 4 > cellCoordinates = {{vertices_.at( vertexID0.getID() )->getCoordinates(), + vertices_.at( vertexID1.getID() )->getCoordinates(), + vertices_.at( vertexID2.getID() )->getCoordinates(), + vertices_.at( vertexID3.getID() )->getCoordinates()}}; + + std::array< std::map< uint_t, uint_t >, 6 > edgeLocalVertexToCellLocalVertexMaps; + + // edgeLocalVertexToCellLocalVertexMaps[ cellLocalEdgeID ][ edgeLocalVertexID ] = cellLocalVertexID; + + edgeLocalVertexToCellLocalVertexMaps[0][edges_.at( edgeID0.getID() )->vertex_index( vertexID0 )] = 0; + edgeLocalVertexToCellLocalVertexMaps[0][edges_.at( edgeID0.getID() )->vertex_index( vertexID1 )] = 1; + + edgeLocalVertexToCellLocalVertexMaps[1][edges_.at( edgeID1.getID() )->vertex_index( vertexID0 )] = 0; + edgeLocalVertexToCellLocalVertexMaps[1][edges_.at( edgeID1.getID() )->vertex_index( vertexID2 )] = 2; + + edgeLocalVertexToCellLocalVertexMaps[2][edges_.at( edgeID2.getID() )->vertex_index( vertexID1 )] = 1; + edgeLocalVertexToCellLocalVertexMaps[2][edges_.at( edgeID2.getID() )->vertex_index( vertexID2 )] = 2; + + edgeLocalVertexToCellLocalVertexMaps[3][edges_.at( edgeID3.getID() )->vertex_index( vertexID0 )] = 0; + edgeLocalVertexToCellLocalVertexMaps[3][edges_.at( edgeID3.getID() )->vertex_index( vertexID3 )] = 3; + + edgeLocalVertexToCellLocalVertexMaps[4][edges_.at( edgeID4.getID() )->vertex_index( vertexID1 )] = 1; + edgeLocalVertexToCellLocalVertexMaps[4][edges_.at( edgeID4.getID() )->vertex_index( vertexID3 )] = 3; + + edgeLocalVertexToCellLocalVertexMaps[5][edges_.at( edgeID5.getID() )->vertex_index( vertexID2 )] = 2; + edgeLocalVertexToCellLocalVertexMaps[5][edges_.at( edgeID5.getID() )->vertex_index( vertexID3 )] = 3; + + std::array< std::map< uint_t, uint_t >, 4 > faceLocalVertexToCellLocalVertexMaps; + + // faceLocalVertexToCellLocalVertexMaps[ cellLocalFaceID ][ faceLocalVertexID ] = cellLocalVertexID; + + faceLocalVertexToCellLocalVertexMaps[0][faces_.at( faceID0.getID() )->vertex_index( vertexID0 )] = 0; + faceLocalVertexToCellLocalVertexMaps[0][faces_.at( faceID0.getID() )->vertex_index( vertexID1 )] = 1; + faceLocalVertexToCellLocalVertexMaps[0][faces_.at( faceID0.getID() )->vertex_index( vertexID2 )] = 2; + + faceLocalVertexToCellLocalVertexMaps[1][faces_.at( faceID1.getID() )->vertex_index( vertexID0 )] = 0; + faceLocalVertexToCellLocalVertexMaps[1][faces_.at( faceID1.getID() )->vertex_index( vertexID1 )] = 1; + faceLocalVertexToCellLocalVertexMaps[1][faces_.at( faceID1.getID() )->vertex_index( vertexID3 )] = 3; + + faceLocalVertexToCellLocalVertexMaps[2][faces_.at( faceID2.getID() )->vertex_index( vertexID0 )] = 0; + faceLocalVertexToCellLocalVertexMaps[2][faces_.at( faceID2.getID() )->vertex_index( vertexID2 )] = 2; + faceLocalVertexToCellLocalVertexMaps[2][faces_.at( faceID2.getID() )->vertex_index( vertexID3 )] = 3; + + faceLocalVertexToCellLocalVertexMaps[3][faces_.at( faceID3.getID() )->vertex_index( vertexID1 )] = 1; + faceLocalVertexToCellLocalVertexMaps[3][faces_.at( faceID3.getID() )->vertex_index( vertexID2 )] = 2; + faceLocalVertexToCellLocalVertexMaps[3][faces_.at( faceID3.getID() )->vertex_index( vertexID3 )] = 3; + + cells_[cellID.getID()] = std::make_shared< Cell >( cellID, + cellVertices, + cellEdges, + cellFaces, + cellCoordinates, + edgeLocalVertexToCellLocalVertexMaps, + faceLocalVertexToCellLocalVertexMaps ); + + setMeshBoundaryFlag( cellID.getID(), meshInfoCell.getBoundaryFlag() ); + } + + // add indirect neighbor faces + for ( auto& it : faces_ ) + { + auto faceID = it.first; + auto face = it.second; + + std::set< PrimitiveID > indirectNeighborsSet; + + for ( const auto& vertexID : face->neighborVertices() ) + { + auto vertex = getVertex( vertexID ); + for ( const auto& neighborFaceID : vertex->neighborFaces() ) + { + if ( neighborFaceID != faceID ) + { + indirectNeighborsSet.insert( neighborFaceID ); + } + } + } + + face->indirectNeighborFaceIDs_.clear(); + face->indirectNeighborFaceIDs_.insert( + face->indirectNeighborFaceIDs_.begin(), indirectNeighborsSet.begin(), indirectNeighborsSet.end() ); + } + + // add indirect neighbor cells + for ( auto& it : cells_ ) + { + auto cellID = it.first; + auto cell = it.second; + + std::set< PrimitiveID > indirectNeighborsSet; + + for ( const auto& vertexID : cell->neighborVertices() ) + { + auto vertex = getVertex( vertexID ); + for ( const auto& neighborCellID : vertex->neighborCells() ) + { + if ( neighborCellID != cellID ) + { + indirectNeighborsSet.insert( neighborCellID ); + } + } + } + + cell->indirectNeighborCellIDs_.clear(); + cell->indirectNeighborCellIDs_.insert( + cell->indirectNeighborCellIDs_.begin(), indirectNeighborsSet.begin(), indirectNeighborsSet.end() ); + } + + loadbalancing::roundRobin( *this ); } SetupPrimitiveStorage::SetupPrimitiveStorage( const VertexMap& vertices, @@ -406,28 +431,52 @@ SetupPrimitiveStorage::SetupPrimitiveStorage( const VertexMap& vertices, loadbalancing::roundRobin( *this ); } -Primitive * SetupPrimitiveStorage::getPrimitive( const PrimitiveID & id ) +Primitive* SetupPrimitiveStorage::getPrimitive( const PrimitiveID& id ) { - if ( vertexExists( id ) ) { return getVertex( id ); } - if ( edgeExists( id ) ) { return getEdge( id ); } - if ( faceExists( id ) ) { return getFace( id ); } - if ( cellExists( id ) ) { return getCell( id ); } - return nullptr; + if ( vertexExists( id ) ) + { + return getVertex( id ); + } + if ( edgeExists( id ) ) + { + return getEdge( id ); + } + if ( faceExists( id ) ) + { + return getFace( id ); + } + if ( cellExists( id ) ) + { + return getCell( id ); + } + return nullptr; } -const Primitive * SetupPrimitiveStorage::getPrimitive( const PrimitiveID & id ) const +const Primitive* SetupPrimitiveStorage::getPrimitive( const PrimitiveID& id ) const { - if ( vertexExists( id ) ) { return getVertex( id ); } - if ( edgeExists( id ) ) { return getEdge( id ); } - if ( faceExists( id ) ) { return getFace( id ); } - if ( cellExists( id ) ) { return getCell( id ); } - return nullptr; + if ( vertexExists( id ) ) + { + return getVertex( id ); + } + if ( edgeExists( id ) ) + { + return getEdge( id ); + } + if ( faceExists( id ) ) + { + return getFace( id ); + } + if ( cellExists( id ) ) + { + return getCell( id ); + } + return nullptr; } uint_t SetupPrimitiveStorage::getNumCellsOnRank( uint_t rank ) const { uint_t n = 0; - for ( const auto & it : cells_ ) + for ( const auto& it : cells_ ) { if ( getTargetRank( it.first ) == rank ) { @@ -440,7 +489,7 @@ uint_t SetupPrimitiveStorage::getNumCellsOnRank( uint_t rank ) const uint_t SetupPrimitiveStorage::getNumFacesOnRank( uint_t rank ) const { uint_t n = 0; - for ( const auto & it : faces_ ) + for ( const auto& it : faces_ ) { if ( getTargetRank( it.first ) == rank ) { @@ -453,7 +502,7 @@ uint_t SetupPrimitiveStorage::getNumFacesOnRank( uint_t rank ) const uint_t SetupPrimitiveStorage::getNumEdgesOnRank( uint_t rank ) const { uint_t n = 0; - for ( const auto & it : edges_ ) + for ( const auto& it : edges_ ) { if ( getTargetRank( it.first ) == rank ) { @@ -466,7 +515,7 @@ uint_t SetupPrimitiveStorage::getNumEdgesOnRank( uint_t rank ) const uint_t SetupPrimitiveStorage::getNumVerticesOnRank( uint_t rank ) const { uint_t n = 0; - for ( const auto & it : vertices_ ) + for ( const auto& it : vertices_ ) { if ( getTargetRank( it.first ) == rank ) { @@ -476,317 +525,390 @@ uint_t SetupPrimitiveStorage::getNumVerticesOnRank( uint_t rank ) const return n; } -void SetupPrimitiveStorage::assembleRankToSetupPrimitivesMap( RankToSetupPrimitivesMap & rankToSetupPrimitivesMap ) const +void SetupPrimitiveStorage::assembleRankToSetupPrimitivesMap( RankToSetupPrimitivesMap& rankToSetupPrimitivesMap ) const { - rankToSetupPrimitivesMap.clear(); + rankToSetupPrimitivesMap.clear(); - PrimitiveMap setupPrimitives; - getSetupPrimitives( setupPrimitives ); - for ( uint_t rank = 0; rank < numberOfProcesses_; rank++ ) - { - rankToSetupPrimitivesMap[ rank ] = std::vector< PrimitiveID::IDType >(); - for ( auto setupPrimitive : setupPrimitives ) - { - if ( rank == getTargetRank( setupPrimitive.first ) ) + PrimitiveMap setupPrimitives; + getSetupPrimitives( setupPrimitives ); + for ( uint_t rank = 0; rank < numberOfProcesses_; rank++ ) + { + rankToSetupPrimitivesMap[rank] = std::vector< PrimitiveID::IDType >(); + for ( auto setupPrimitive : setupPrimitives ) { - rankToSetupPrimitivesMap[ rank ].push_back( setupPrimitive.first ); + if ( rank == getTargetRank( setupPrimitive.first ) ) + { + rankToSetupPrimitivesMap[rank].push_back( setupPrimitive.first ); + } } - } - } + } - WALBERLA_ASSERT_LESS_EQUAL( rankToSetupPrimitivesMap.size(), numberOfProcesses_ ); + WALBERLA_ASSERT_LESS_EQUAL( rankToSetupPrimitivesMap.size(), numberOfProcesses_ ); } uint_t SetupPrimitiveStorage::getNumberOfEmptyProcesses() const { - uint_t numberOfEmptyProcesses = 0; - RankToSetupPrimitivesMap rankToSetupPrimitivesMap; - assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); - for ( auto const & rankToSetupPrimitives : rankToSetupPrimitivesMap ) - { - if ( rankToSetupPrimitives.second.size() == 0 ) - { - numberOfEmptyProcesses++; - } - } - return numberOfEmptyProcesses; + uint_t numberOfEmptyProcesses = 0; + RankToSetupPrimitivesMap rankToSetupPrimitivesMap; + assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); + for ( auto const& rankToSetupPrimitives : rankToSetupPrimitivesMap ) + { + if ( rankToSetupPrimitives.second.size() == 0 ) + { + numberOfEmptyProcesses++; + } + } + return numberOfEmptyProcesses; } uint_t SetupPrimitiveStorage::getMinPrimitivesPerRank() const { - uint_t minNumberOfPrimitives = std::numeric_limits< uint_t >::max(); - RankToSetupPrimitivesMap rankToSetupPrimitivesMap; - assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); - for ( auto const & rankToSetupPrimitives : rankToSetupPrimitivesMap ) - { - minNumberOfPrimitives = std::min( rankToSetupPrimitives.second.size(), minNumberOfPrimitives ); - } - return minNumberOfPrimitives; + uint_t minNumberOfPrimitives = std::numeric_limits< uint_t >::max(); + RankToSetupPrimitivesMap rankToSetupPrimitivesMap; + assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); + for ( auto const& rankToSetupPrimitives : rankToSetupPrimitivesMap ) + { + minNumberOfPrimitives = std::min( rankToSetupPrimitives.second.size(), minNumberOfPrimitives ); + } + return minNumberOfPrimitives; } uint_t SetupPrimitiveStorage::getMaxPrimitivesPerRank() const { - uint_t maxNumberOfPrimitives = 0; - RankToSetupPrimitivesMap rankToSetupPrimitivesMap; - assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); - for ( auto const & rankToSetupPrimitives : rankToSetupPrimitivesMap ) - { - maxNumberOfPrimitives = std::max( rankToSetupPrimitives.second.size(), maxNumberOfPrimitives ); - } - return maxNumberOfPrimitives; + uint_t maxNumberOfPrimitives = 0; + RankToSetupPrimitivesMap rankToSetupPrimitivesMap; + assembleRankToSetupPrimitivesMap( rankToSetupPrimitivesMap ); + for ( auto const& rankToSetupPrimitives : rankToSetupPrimitivesMap ) + { + maxNumberOfPrimitives = std::max( rankToSetupPrimitives.second.size(), maxNumberOfPrimitives ); + } + return maxNumberOfPrimitives; } real_t SetupPrimitiveStorage::getAvgPrimitivesPerRank() const { - return real_t( getNumberOfPrimitives() ) / real_t( numberOfProcesses_ ); + return real_t( getNumberOfPrimitives() ) / real_t( numberOfProcesses_ ); } uint_t SetupPrimitiveStorage::getNumVerticesOnBoundary() const { - return uint_c( std::count_if( vertices_.begin(), vertices_.end(), - [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Vertex > > & mapEntry ) - { return onBoundary( PrimitiveID( mapEntry.first ), true ); } ) ); + return uint_c( std::count_if( + vertices_.begin(), vertices_.end(), [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Vertex > >& mapEntry ) { + return onBoundary( PrimitiveID( mapEntry.first ), true ); + } ) ); } uint_t SetupPrimitiveStorage::getNumEdgesOnBoundary() const { - return uint_c( std::count_if( edges_.begin(), edges_.end(), - [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Edge > > & mapEntry ) - { return onBoundary( PrimitiveID( mapEntry.first ), true ); } ) ); + return uint_c( std::count_if( + edges_.begin(), edges_.end(), [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Edge > >& mapEntry ) { + return onBoundary( PrimitiveID( mapEntry.first ), true ); + } ) ); } uint_t SetupPrimitiveStorage::getNumFacesOnBoundary() const { - return uint_c( std::count_if( faces_.begin(), faces_.end(), - [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Face > > & mapEntry ) - { return onBoundary( PrimitiveID( mapEntry.first ), true ); } ) ); + return uint_c( std::count_if( + faces_.begin(), faces_.end(), [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Face > >& mapEntry ) { + return onBoundary( PrimitiveID( mapEntry.first ), true ); + } ) ); } uint_t SetupPrimitiveStorage::getNumCellsOnBoundary() const { - return uint_c( std::count_if( cells_.begin(), cells_.end(), - [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Cell > > & mapEntry ) - { return onBoundary( PrimitiveID( mapEntry.first ), true ); } ) ); + return uint_c( std::count_if( + cells_.begin(), cells_.end(), [this]( const std::pair< PrimitiveID::IDType, std::shared_ptr< Cell > >& mapEntry ) { + return onBoundary( PrimitiveID( mapEntry.first ), true ); + } ) ); } -void SetupPrimitiveStorage::toStream( std::ostream & os, bool verbose ) const +void SetupPrimitiveStorage::toStream( std::ostream& os, bool verbose ) const { - os << "SetupPrimitiveStorage:\n"; - - os << " - Processes (overall): " << std::setw(10) << numberOfProcesses_ << "\n"; - os << " - Processes (empty) : " << std::setw(10) << getNumberOfEmptyProcesses() << "\n"; - - os << " - Number of...\n" - << " + Vertices: " << std::setw(10) << vertices_.size() << " | " << getNumVerticesOnBoundary() << "/" << vertices_.size() << " on boundary\n" - << " + Edges: " << std::setw(10) << edges_.size() << " | " << getNumEdgesOnBoundary() << "/" << edges_.size() << " on boundary\n" - << " + Faces: " << std::setw(10) << faces_.size() << " | " << getNumFacesOnBoundary() << "/" << faces_.size() << " on boundary\n" - << " + Cells: " << std::setw(10) << cells_.size() << " | " << getNumCellsOnBoundary() << "/" << cells_.size() << " on boundary\n"; - - os << " - Primitives per process...\n" - << " + min: " << std::setw(10) << getMinPrimitivesPerRank() << "\n" - << " + max: " << std::setw(10) << getMaxPrimitivesPerRank() << "\n" - << " + avg: " << std::setw(10) << getAvgPrimitivesPerRank() << ""; - - if ( verbose ) { - os << "\n"; - os << "Vertices: ID | Target Rank | Position | Neighbor Edges \n" - << "---------------------------------------------------------\n"; - for ( auto it = vertices_.begin(); it != vertices_.end(); it++ ) { - Point3D coordinates = it->second->getCoordinates(); - os << " " << std::setw( 4 ) << it->first << " | " - << std::setw( 11 ) << getTargetRank( it->first ) << " | " - << coordinates << " | "; - for ( const auto & neighborEdgeID : it->second->getHigherDimNeighbors()) { - os << neighborEdgeID.getID() << " "; + os << "SetupPrimitiveStorage:\n"; + + os << " - Processes (overall): " << std::setw( 10 ) << numberOfProcesses_ << "\n"; + os << " - Processes (empty) : " << std::setw( 10 ) << getNumberOfEmptyProcesses() << "\n"; + + os << " - Number of...\n" + << " + Vertices: " << std::setw( 10 ) << vertices_.size() << " | " << getNumVerticesOnBoundary() << "/" + << vertices_.size() << " on boundary\n" + << " + Edges: " << std::setw( 10 ) << edges_.size() << " | " << getNumEdgesOnBoundary() << "/" << edges_.size() + << " on boundary\n" + << " + Faces: " << std::setw( 10 ) << faces_.size() << " | " << getNumFacesOnBoundary() << "/" << faces_.size() + << " on boundary\n" + << " + Cells: " << std::setw( 10 ) << cells_.size() << " | " << getNumCellsOnBoundary() << "/" << cells_.size() + << " on boundary\n"; + + os << " - Primitives per process...\n" + << " + min: " << std::setw( 10 ) << getMinPrimitivesPerRank() << "\n" + << " + max: " << std::setw( 10 ) << getMaxPrimitivesPerRank() << "\n" + << " + avg: " << std::setw( 10 ) << getAvgPrimitivesPerRank() << ""; + + if ( verbose ) + { + os << "\n"; + os << " | | mesh | | \n" + << "Vertices: ID | Target Rank | boundary | Position | Neighbor Edges\n" + << " | | flag | | \n" + << "-----------------------------------------------------------------------------------------------\n"; + for ( auto it = vertices_.begin(); it != vertices_.end(); it++ ) + { + Point3D coordinates = it->second->getCoordinates(); + os << " " << std::setw( 4 ) << it->first << " | " << std::setw( 11 ) << getTargetRank( it->first ) << " | " + << std::setw( 8 ) << it->second->getMeshBoundaryFlag() << " | " << std::scientific << std::showpos + << std::setprecision( 3 ) << coordinates << " | "; + for ( const auto& neighborEdgeID : it->second->getHigherDimNeighbors() ) + { + os << neighborEdgeID.getID() << " "; + } + os << "\n"; } os << "\n"; - } - os << "\n"; - - os << "Edges: ID | Target Rank | VertexID_0 | VertexID_1 | mesh boundary flag | Neighbor Faces \n" - << "----------------------------------------------------------------------------------------------\n"; - for ( auto it = edges_.begin(); it != edges_.end(); it++ ) { - os << " " << std::setw( 4 ) << it->first << " | " - << std::setw( 11 ) << getTargetRank( it->first ) << " | " - << std::setw( 10 ) << it->second->getVertexID0().getID() << " | " - << std::setw( 10 ) << it->second->getVertexID1().getID() << " | " - << std::setw( 20 ) << it->second->getMeshBoundaryFlag() << " | "; - for ( const auto & neighborFaceID : it->second->getHigherDimNeighbors()) { - os << neighborFaceID.getID() << " "; + os << " | | | | mesh |\n" + << "Edges: ID | Target Rank | VertexID_0 | VertexID_1 | boundary | Neighbor Faces\n" + << " | | | | flag |\n" + << "----------------------------------------------------------------------------------\n"; + for ( auto it = edges_.begin(); it != edges_.end(); it++ ) + { + os << " " << std::setw( 4 ) << it->first << " | " << std::setw( 11 ) << getTargetRank( it->first ) << " | " + << std::setw( 10 ) << it->second->getVertexID0().getID() << " | " << std::setw( 10 ) + << it->second->getVertexID1().getID() << " | " << std::setw( 8 ) << it->second->getMeshBoundaryFlag() << " | "; + for ( const auto& neighborFaceID : it->second->getHigherDimNeighbors() ) + { + os << neighborFaceID.getID() << " "; + } + os << "\n"; } os << "\n"; - } - os << "\n"; - - os << "Faces: ID | Target Rank | EdgeID_0 | EdgeID_1 | EdgeID_2\n" - << "-------------------------------------------------------------\n"; - for ( auto it = faces_.begin(); it != faces_.end(); it++ ) { - os << " " << std::setw( 4 ) << it->first << " | " - << std::setw( 11 ) << getTargetRank( it->first ) << " | " - << std::setw( 8 ) << it->second->getEdgeID0().getID() << " | " - << std::setw( 8 ) << it->second->getEdgeID1().getID() << " | " - << std::setw( 8 ) << it->second->getEdgeID2().getID() << "\n"; - } - os << "\n"; - - if ( cells_.size() > 0 ) { - os << "Cells: ID | Target Rank | FaceID_0 | FaceID_1 | FaceID_2 | FaceID_3\n" + + os << " | | mesh | | | \n" + << "Faces: ID | Target Rank | boundary | EdgeID_0 | EdgeID_1 | EdgeID_2\n" + << " | | flag | | | \n" << "------------------------------------------------------------------------\n"; - for ( auto it = cells_.begin(); it != cells_.end(); it++ ) { - os << " " << std::setw( 4 ) << it->first << " | " - << std::setw( 11 ) << getTargetRank( it->first ) << " | " - << std::setw( 8 ) << it->second->neighborFaces()[0].getID() << " | " - << std::setw( 8 ) << it->second->neighborFaces()[1].getID() << " | " - << std::setw( 8 ) << it->second->neighborFaces()[2].getID() << " | " - << std::setw( 8 ) << it->second->neighborFaces()[3].getID() << "\n"; + for ( auto it = faces_.begin(); it != faces_.end(); it++ ) + { + os << " " << std::setw( 4 ) << it->first << " | " << std::setw( 11 ) << getTargetRank( it->first ) << " | " + << std::setw( 8 ) << it->second->getMeshBoundaryFlag() << " | " << std::setw( 8 ) << it->second->getEdgeID0().getID() + << " | " << std::setw( 8 ) << it->second->getEdgeID1().getID() << " | " << std::setw( 8 ) + << it->second->getEdgeID2().getID() << "\n"; } - } - } + os << "\n"; + if ( cells_.size() > 0 ) + { + os << " | | mesh | | | |\n" + << "Cells: ID | Target Rank | boundary | FaceID_0 | FaceID_1 | FaceID_2 | FaceID_3\n" + << " | | flag | | | |\n" + << "--------------------------------------------------------------------------------------\n"; + for ( auto it = cells_.begin(); it != cells_.end(); it++ ) + { + os << " " << std::setw( 4 ) << it->first << " | " << std::setw( 11 ) << getTargetRank( it->first ) << " | " + << std::setw( 8 ) << it->second->getMeshBoundaryFlag() << " | " << std::setw( 8 ) + << it->second->neighborFaces()[0].getID() << " | " << std::setw( 8 ) << it->second->neighborFaces()[1].getID() + << " | " << std::setw( 8 ) << it->second->neighborFaces()[2].getID() << " | " << std::setw( 8 ) + << it->second->neighborFaces()[3].getID() << "\n"; + } + } + } } - -void SetupPrimitiveStorage::getSetupPrimitives( PrimitiveMap & setupPrimitiveMap ) const +void SetupPrimitiveStorage::getSetupPrimitives( PrimitiveMap& setupPrimitiveMap ) const { - setupPrimitiveMap.clear(); + setupPrimitiveMap.clear(); - setupPrimitiveMap.insert( vertices_.begin(), vertices_.end() ); - setupPrimitiveMap.insert( edges_.begin(), edges_.end() ); - setupPrimitiveMap.insert( faces_.begin(), faces_.end() ); - setupPrimitiveMap.insert( cells_.begin(), cells_.end() ); + setupPrimitiveMap.insert( vertices_.begin(), vertices_.end() ); + setupPrimitiveMap.insert( edges_.begin(), edges_.end() ); + setupPrimitiveMap.insert( faces_.begin(), faces_.end() ); + setupPrimitiveMap.insert( cells_.begin(), cells_.end() ); - WALBERLA_ASSERT_EQUAL( setupPrimitiveMap.size(), vertices_.size() + edges_.size() + faces_.size() + cells_.size() ); + WALBERLA_ASSERT_EQUAL( setupPrimitiveMap.size(), vertices_.size() + edges_.size() + faces_.size() + cells_.size() ); } - PrimitiveID SetupPrimitiveStorage::generatePrimitiveID() const { - PrimitiveID newID( getNumberOfPrimitives() ); - WALBERLA_ASSERT( !primitiveExists( newID ) ); - return newID; + PrimitiveID newID( getNumberOfPrimitives() ); + WALBERLA_ASSERT( !primitiveExists( newID ) ); + return newID; } -void SetupPrimitiveStorage::setMeshBoundaryFlagsOnBoundary( const uint_t & meshBoundaryFlagOnBoundary, const uint_t & meshBoundaryFlagInner, const bool & highestDimensionAlwaysInner ) +void SetupPrimitiveStorage::setMeshBoundaryFlagsOnBoundary( const uint_t& meshBoundaryFlagOnBoundary, + const uint_t& meshBoundaryFlagInner, + const bool& highestDimensionAlwaysInner ) { - PrimitiveMap primitives; - getSetupPrimitives( primitives ); - for ( auto & primitive : primitives ) - { - primitive.second->meshBoundaryFlag_ = onBoundary( primitive.first, highestDimensionAlwaysInner ) ? meshBoundaryFlagOnBoundary : meshBoundaryFlagInner; - } + PrimitiveMap primitives; + getSetupPrimitives( primitives ); + for ( auto& primitive : primitives ) + { + primitive.second->meshBoundaryFlag_ = + onBoundary( primitive.first, highestDimensionAlwaysInner ) ? meshBoundaryFlagOnBoundary : meshBoundaryFlagInner; + } } -void SetupPrimitiveStorage::setMeshBoundaryFlagsInner( const uint_t & meshBoundaryFlagInner, const bool & highestDimensionAlwaysInner ) +void SetupPrimitiveStorage::setMeshBoundaryFlagsInner( const uint_t& meshBoundaryFlagInner, + const bool& highestDimensionAlwaysInner ) { - PrimitiveMap primitives; - getSetupPrimitives( primitives ); - for ( auto & primitive : primitives ) - { - primitive.second->meshBoundaryFlag_ = onBoundary( primitive.first, highestDimensionAlwaysInner ) ? primitive.second->getMeshBoundaryFlag() : meshBoundaryFlagInner; - } + PrimitiveMap primitives; + getSetupPrimitives( primitives ); + for ( auto& primitive : primitives ) + { + primitive.second->meshBoundaryFlag_ = onBoundary( primitive.first, highestDimensionAlwaysInner ) ? + primitive.second->getMeshBoundaryFlag() : + meshBoundaryFlagInner; + } } - -void SetupPrimitiveStorage::setMeshBoundaryFlagsByVertexLocation( const uint_t & meshBoundaryFlag, - const std::function< bool( const Point3D & x ) > & onBoundary, - const bool & allVertices ) +void SetupPrimitiveStorage::setMeshBoundaryFlagsByVertexLocation( const uint_t& meshBoundaryFlag, + const std::function< bool( const Point3D& x ) >& onBoundary, + const bool& allVertices ) { - auto cond = [ allVertices, onBoundary ]( const std::vector< Point3D > & coordinates ) { - if ( allVertices ) - return std::all_of( coordinates.begin(), coordinates.end(), onBoundary ); - else - return std::any_of( coordinates.begin(), coordinates.end(), onBoundary ); - }; - - for ( const auto & p : vertices_ ) - { - if ( cond( { p.second->getCoordinates() } ) ) - setMeshBoundaryFlag( p.first, meshBoundaryFlag ); - } - - for ( const auto & p : edges_ ) - { - if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) - setMeshBoundaryFlag( p.first, meshBoundaryFlag ); - } - - for ( const auto & p : faces_ ) - { - if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) - setMeshBoundaryFlag( p.first, meshBoundaryFlag ); - } - - for ( const auto & p : cells_ ) - { - if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) - setMeshBoundaryFlag( p.first, meshBoundaryFlag ); - } -} + auto cond = [allVertices, onBoundary]( const std::vector< Point3D >& coordinates ) { + if ( allVertices ) + return std::all_of( coordinates.begin(), coordinates.end(), onBoundary ); + else + return std::any_of( coordinates.begin(), coordinates.end(), onBoundary ); + }; + + for ( const auto& p : vertices_ ) + { + if ( cond( {p.second->getCoordinates()} ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } + for ( const auto& p : edges_ ) + { + if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } + + for ( const auto& p : faces_ ) + { + if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } -bool SetupPrimitiveStorage::onBoundary( const PrimitiveID & primitiveID, const bool & highestDimensionAlwaysInner ) const + for ( const auto& p : cells_ ) + { + if ( cond( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ) ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } +} + +void SetupPrimitiveStorage::setMeshBoundaryFlagsByCentroidLocation( const uint_t& meshBoundaryFlag, + const std::function< bool( const Point3D& x ) >& onBoundary, + bool useGeometryMap ) { - WALBERLA_ASSERT( primitiveExists( primitiveID ) ); - - if ( getNumberOfCells() == 0 ) - { - // 2D - if ( highestDimensionAlwaysInner && faceExists( primitiveID ) ) - { - return false; - } - if ( edgeExists( primitiveID ) ) - { - const auto edge = getEdge( primitiveID ); - WALBERLA_ASSERT_GREATER( edge->getNumNeighborFaces(), 0 ); - WALBERLA_ASSERT_LESS_EQUAL( edge->getNumNeighborFaces(), 2 ); - return edge->getNumNeighborFaces() == 1; - } - else - { - const auto primitive = getPrimitive( primitiveID ); - std::vector< PrimitiveID > neighborEdges; - primitive->getNeighborEdges( neighborEdges ); - for ( auto it : neighborEdges ) - { - if ( onBoundary( it ) ) - { - return true; - } - } - return false; - } - } - else - { - // 3D - if ( highestDimensionAlwaysInner && cellExists( primitiveID ) ) - { - return false; - } - if ( faceExists( primitiveID ) ) - { - const auto face = getFace( primitiveID ); - WALBERLA_ASSERT_GREATER( face->getNumNeighborCells(), 0 ); - WALBERLA_ASSERT_LESS_EQUAL( face->getNumNeighborCells(), 2 ); - return face->getNumNeighborCells() == 1; - } - else - { - const auto primitive = getPrimitive( primitiveID ); - std::vector< PrimitiveID > neighborFaces; - primitive->getNeighborFaces( neighborFaces ); - for ( auto it : neighborFaces ) - { - if ( onBoundary( it ) ) - { - return true; - } - } - return false; - } - } + auto centroid = [useGeometryMap]( const std::vector< Point3D >& coordinates, + const std::shared_ptr< GeometryMap >& map ) -> Point3D { + Point3D c( {real_c( 0 ), real_c( 0 ), real_c( 0 )} ); + for ( const auto& p : coordinates ) + { + c += p; + } + c *= real_c( 1 ) / real_c( coordinates.size() ); + if ( useGeometryMap ) + { + Point3D cMapped; + map->evalF( c, cMapped ); + return cMapped; + } + return c; + }; + + for ( const auto& p : vertices_ ) + { + const auto c = centroid( {p.second->getCoordinates()}, p.second->getGeometryMap() ); + if ( onBoundary( c ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } + + for ( const auto& p : edges_ ) + { + const auto c = centroid( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ), + p.second->getGeometryMap() ); + if ( onBoundary( c ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } + + for ( const auto& p : faces_ ) + { + const auto c = centroid( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ), + p.second->getGeometryMap() ); + if ( onBoundary( c ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } + + for ( const auto& p : cells_ ) + { + const auto c = centroid( std::vector< Point3D >( p.second->getCoordinates().begin(), p.second->getCoordinates().end() ), + p.second->getGeometryMap() ); + if ( onBoundary( c ) ) + setMeshBoundaryFlag( p.first, meshBoundaryFlag ); + } } +bool SetupPrimitiveStorage::onBoundary( const PrimitiveID& primitiveID, const bool& highestDimensionAlwaysInner ) const +{ + WALBERLA_ASSERT( primitiveExists( primitiveID ) ); + + if ( getNumberOfCells() == 0 ) + { + // 2D + if ( highestDimensionAlwaysInner && faceExists( primitiveID ) ) + { + return false; + } + if ( edgeExists( primitiveID ) ) + { + const auto edge = getEdge( primitiveID ); + WALBERLA_ASSERT_GREATER( edge->getNumNeighborFaces(), 0 ); + WALBERLA_ASSERT_LESS_EQUAL( edge->getNumNeighborFaces(), 2 ); + return edge->getNumNeighborFaces() == 1; + } + else + { + const auto primitive = getPrimitive( primitiveID ); + std::vector< PrimitiveID > neighborEdges; + primitive->getNeighborEdges( neighborEdges ); + for ( auto it : neighborEdges ) + { + if ( onBoundary( it ) ) + { + return true; + } + } + return false; + } + } + else + { + // 3D + if ( highestDimensionAlwaysInner && cellExists( primitiveID ) ) + { + return false; + } + if ( faceExists( primitiveID ) ) + { + const auto face = getFace( primitiveID ); + WALBERLA_ASSERT_GREATER( face->getNumNeighborCells(), 0 ); + WALBERLA_ASSERT_LESS_EQUAL( face->getNumNeighborCells(), 2 ); + return face->getNumNeighborCells() == 1; + } + else + { + const auto primitive = getPrimitive( primitiveID ); + std::vector< PrimitiveID > neighborFaces; + primitive->getNeighborFaces( neighborFaces ); + for ( auto it : neighborFaces ) + { + if ( onBoundary( it ) ) + { + return true; + } + } + return false; + } + } } + +} // namespace hyteg diff --git a/src/hyteg/primitivestorage/SetupPrimitiveStorage.hpp b/src/hyteg/primitivestorage/SetupPrimitiveStorage.hpp index e8e94a91c7e599e0c303753c89a257424cb82a78..9f019608823ad4fa752b5e260464da7c9e500960 100644 --- a/src/hyteg/primitivestorage/SetupPrimitiveStorage.hpp +++ b/src/hyteg/primitivestorage/SetupPrimitiveStorage.hpp @@ -20,155 +20,178 @@ #pragma once -#include "core/debug/Debug.h" -#include "hyteg/mesh/MeshInfo.hpp" -#include "hyteg/PrimitiveID.hpp" -#include "hyteg/primitives/Primitive.hpp" -#include "hyteg/primitives/Vertex.hpp" -#include "hyteg/primitives/Edge.hpp" -#include "hyteg/primitives/Face.hpp" -#include "hyteg/primitives/Cell.hpp" - #include #include #include #include +#include "core/debug/Debug.h" + +#include "hyteg/PrimitiveID.hpp" +#include "hyteg/mesh/MeshInfo.hpp" +#include "hyteg/primitives/Cell.hpp" +#include "hyteg/primitives/Edge.hpp" +#include "hyteg/primitives/Face.hpp" +#include "hyteg/primitives/Primitive.hpp" +#include "hyteg/primitives/Vertex.hpp" + namespace hyteg { -using walberla::real_t; using walberla::memory_t; +using walberla::real_t; class SetupPrimitiveStorage { -public: - - typedef std::map< PrimitiveID::IDType, std::shared_ptr< Primitive > > PrimitiveMap; - typedef std::map< PrimitiveID::IDType, std::shared_ptr< Vertex > > VertexMap; - typedef std::map< PrimitiveID::IDType, std::shared_ptr< Edge > > EdgeMap; - typedef std::map< PrimitiveID::IDType, std::shared_ptr< Face > > FaceMap; - typedef std::map< PrimitiveID::IDType, std::shared_ptr< Cell > > CellMap; - - SetupPrimitiveStorage( const MeshInfo & meshInfo, const uint_t & numberOfProcesses ); - - SetupPrimitiveStorage( const VertexMap& vertices, - const EdgeMap& edges, - const FaceMap& faces, - const CellMap& cells, - const uint_t& numberOfProcesses ); - - void toStream( std::ostream& os, bool verbose = false ) const; - - uint_t getNumberOfProcesses() const { return numberOfProcesses_; } - uint_t getNumberOfEmptyProcesses() const; - - bool primitiveExists( const PrimitiveID & id ) const { return vertexExists( id ) || edgeExists( id ) || faceExists( id ) || cellExists( id ); } - bool vertexExists ( const PrimitiveID & id ) const { return vertices_.count( id.getID() ) > 0; } - bool edgeExists ( const PrimitiveID & id ) const { return edges_.count( id.getID() ) > 0; } - bool faceExists ( const PrimitiveID & id ) const { return faces_.count( id.getID() ) > 0; } - bool cellExists ( const PrimitiveID & id ) const { return cells_.count( id.getID() ) > 0; } - - Primitive * getPrimitive( const PrimitiveID & id ); - const Primitive * getPrimitive( const PrimitiveID & id ) const; - Vertex * getVertex( const PrimitiveID & id ) { return vertexExists( id ) ? vertices_.at( id.getID() ).get() : nullptr; } - const Vertex * getVertex( const PrimitiveID & id ) const { return vertexExists( id ) ? vertices_.at( id.getID() ).get() : nullptr; } - Edge * getEdge ( const PrimitiveID & id ) { return edgeExists( id ) ? edges_.at( id.getID() ).get() : nullptr; } - const Edge * getEdge ( const PrimitiveID & id ) const { return edgeExists( id ) ? edges_.at( id.getID() ).get() : nullptr; } - Face * getFace ( const PrimitiveID & id ) { return faceExists( id ) ? faces_.at( id.getID() ).get() : nullptr; } - const Face * getFace ( const PrimitiveID & id ) const { return faceExists( id ) ? faces_.at( id.getID() ).get() : nullptr; } - Cell * getCell ( const PrimitiveID & id ) { return cellExists( id ) ? cells_.at( id.getID() ).get() : nullptr; } - const Cell * getCell ( const PrimitiveID & id ) const { return cellExists( id ) ? cells_.at( id.getID() ).get() : nullptr; } - - void getSetupPrimitives( PrimitiveMap & setupPrimitiveMap ) const; - - uint_t getNumberOfPrimitives() const { return getNumberOfVertices() + getNumberOfEdges() + getNumberOfFaces() + getNumberOfCells(); } - uint_t getNumberOfVertices () const { return vertices_.size(); } - uint_t getNumberOfEdges () const { return edges_.size(); } - uint_t getNumberOfFaces () const { return faces_.size(); } - uint_t getNumberOfCells () const { return cells_.size(); } - - /// Returns a reference to a map of \ref Vertex instances - const VertexMap & getVertices() const { return vertices_; } - - /// Returns a reference to a map of \ref Edge instances - const EdgeMap & getEdges() const { return edges_; } - - /// Returns a reference to a map of \ref Face instances - const FaceMap & getFaces() const { return faces_; } - - /// Returns a reference to a map of \ref Cell instances - const CellMap & getCells() const { return cells_; } - - void setTargetRank( const PrimitiveID & primitiveID, const uint_t & targetRank ) { primitiveIDToTargetRankMap_[ primitiveID.getID() ] = targetRank; } - uint_t getTargetRank( const PrimitiveID & primitiveID ) const { return primitiveIDToTargetRankMap_.at( primitiveID.getID() ); } - - uint_t getNumCellsOnRank( uint_t rank ) const; - uint_t getNumFacesOnRank( uint_t rank ) const; - uint_t getNumEdgesOnRank( uint_t rank ) const; - uint_t getNumVerticesOnRank( uint_t rank ) const; - - void setGeometryMap( const PrimitiveID & primitiveID, const std::shared_ptr& map) { getPrimitive(primitiveID)->geometryMap_ = map; } - - void setMeshBoundaryFlag( const PrimitiveID & primitiveID, const uint_t & meshBoundaryFlag ) { getPrimitive(primitiveID)->meshBoundaryFlag_ = meshBoundaryFlag; } - - /// Sets the mesh boundary flag of the primitives to a specified value if they are located on the boundary of the domain - /// \param meshBoundaryFlagOnBoundary the flag the primitives are set to if they are located on the boundary - /// \param meshBoundaryFlagInner the flag the primitives are set to if they are not located on the boundary - /// \param highestDimensionAlwaysInner if true, cells in 3D meshes and faces in 2D meshes are treated as inner primitives - void setMeshBoundaryFlagsOnBoundary( const uint_t & meshBoundaryFlagOnBoundary, const uint_t & meshBoundaryFlagInner, const bool & highestDimensionAlwaysInner ); - - /// Sets the mesh boundary flag of the primitives to a specified value if they are NOT located on the boundary of the domain. - /// This does not change mesh boundary flags on primitives that are not located on the boundary - /// \param meshBoundaryFlagInner the flag the primitives are set to if they are not located on the boundary - /// \param highestDimensionAlwaysInner if true, cells in 3D meshes and faces in 2D meshes are treated as inner primitives - void setMeshBoundaryFlagsInner( const uint_t & meshBoundaryFlagInner, const bool & highestDimensionAlwaysInner ); - - /// Every primitive for which onBoundary() returns true for all / any (if allVertices == true / false) of the primitve's vertices is assigned the passed mesh boundary flag. - void setMeshBoundaryFlagsByVertexLocation( const uint_t & meshBoundaryFlag, - const std::function< bool( const Point3D & x ) > & onBoundary, - const bool & allVertices = true ); - - /// Returns true, if the primitive lies on the boundary. - /// \param primitiveID the ID of the primitive to be tested - /// \param highestDimensionAlwaysInner if true, this method always returns false if the targeted primitive is of highest dimension (cells for 3D, faces for 2D meshes) - bool onBoundary( const PrimitiveID & primitiveID, const bool & highestDimensionAlwaysInner = false ) const; - -private: - - typedef std::map< uint_t, std::vector< PrimitiveID::IDType > > RankToSetupPrimitivesMap; - - PrimitiveID generatePrimitiveID() const; - - void assembleRankToSetupPrimitivesMap( RankToSetupPrimitivesMap & rankToSetupPrimitivesMap ) const; - - /// Returns the number of primitives on the target rank with the least number of primitives - uint_t getMinPrimitivesPerRank() const; - /// Returns the number of primitives on the target rank with the largest number of primitives - uint_t getMaxPrimitivesPerRank() const; - /// Returns the average number of primitives per rank - real_t getAvgPrimitivesPerRank() const; - - /// Returns the number of primitives that lie on the boundary - uint_t getNumVerticesOnBoundary() const; - uint_t getNumEdgesOnBoundary() const; - uint_t getNumFacesOnBoundary() const; - uint_t getNumCellsOnBoundary() const; - - uint_t numberOfProcesses_; - - VertexMap vertices_; - EdgeMap edges_; - FaceMap faces_; - CellMap cells_; - - std::map< PrimitiveID::IDType, uint_t > primitiveIDToTargetRankMap_; - + public: + typedef std::map< PrimitiveID::IDType, std::shared_ptr< Primitive > > PrimitiveMap; + typedef std::map< PrimitiveID::IDType, std::shared_ptr< Vertex > > VertexMap; + typedef std::map< PrimitiveID::IDType, std::shared_ptr< Edge > > EdgeMap; + typedef std::map< PrimitiveID::IDType, std::shared_ptr< Face > > FaceMap; + typedef std::map< PrimitiveID::IDType, std::shared_ptr< Cell > > CellMap; + + SetupPrimitiveStorage( const MeshInfo& meshInfo, const uint_t& numberOfProcesses ); + + SetupPrimitiveStorage( const VertexMap& vertices, + const EdgeMap& edges, + const FaceMap& faces, + const CellMap& cells, + const uint_t& numberOfProcesses ); + + void toStream( std::ostream& os, bool verbose = false ) const; + + uint_t getNumberOfProcesses() const { return numberOfProcesses_; } + uint_t getNumberOfEmptyProcesses() const; + + bool primitiveExists( const PrimitiveID& id ) const + { + return vertexExists( id ) || edgeExists( id ) || faceExists( id ) || cellExists( id ); + } + bool vertexExists( const PrimitiveID& id ) const { return vertices_.count( id.getID() ) > 0; } + bool edgeExists( const PrimitiveID& id ) const { return edges_.count( id.getID() ) > 0; } + bool faceExists( const PrimitiveID& id ) const { return faces_.count( id.getID() ) > 0; } + bool cellExists( const PrimitiveID& id ) const { return cells_.count( id.getID() ) > 0; } + + Primitive* getPrimitive( const PrimitiveID& id ); + const Primitive* getPrimitive( const PrimitiveID& id ) const; + Vertex* getVertex( const PrimitiveID& id ) { return vertexExists( id ) ? vertices_.at( id.getID() ).get() : nullptr; } + const Vertex* getVertex( const PrimitiveID& id ) const + { + return vertexExists( id ) ? vertices_.at( id.getID() ).get() : nullptr; + } + Edge* getEdge( const PrimitiveID& id ) { return edgeExists( id ) ? edges_.at( id.getID() ).get() : nullptr; } + const Edge* getEdge( const PrimitiveID& id ) const { return edgeExists( id ) ? edges_.at( id.getID() ).get() : nullptr; } + Face* getFace( const PrimitiveID& id ) { return faceExists( id ) ? faces_.at( id.getID() ).get() : nullptr; } + const Face* getFace( const PrimitiveID& id ) const { return faceExists( id ) ? faces_.at( id.getID() ).get() : nullptr; } + Cell* getCell( const PrimitiveID& id ) { return cellExists( id ) ? cells_.at( id.getID() ).get() : nullptr; } + const Cell* getCell( const PrimitiveID& id ) const { return cellExists( id ) ? cells_.at( id.getID() ).get() : nullptr; } + + void getSetupPrimitives( PrimitiveMap& setupPrimitiveMap ) const; + + uint_t getNumberOfPrimitives() const + { + return getNumberOfVertices() + getNumberOfEdges() + getNumberOfFaces() + getNumberOfCells(); + } + uint_t getNumberOfVertices() const { return vertices_.size(); } + uint_t getNumberOfEdges() const { return edges_.size(); } + uint_t getNumberOfFaces() const { return faces_.size(); } + uint_t getNumberOfCells() const { return cells_.size(); } + + /// Returns a reference to a map of \ref Vertex instances + const VertexMap& getVertices() const { return vertices_; } + + /// Returns a reference to a map of \ref Edge instances + const EdgeMap& getEdges() const { return edges_; } + + /// Returns a reference to a map of \ref Face instances + const FaceMap& getFaces() const { return faces_; } + + /// Returns a reference to a map of \ref Cell instances + const CellMap& getCells() const { return cells_; } + + void setTargetRank( const PrimitiveID& primitiveID, const uint_t& targetRank ) + { + primitiveIDToTargetRankMap_[primitiveID.getID()] = targetRank; + } + uint_t getTargetRank( const PrimitiveID& primitiveID ) const { return primitiveIDToTargetRankMap_.at( primitiveID.getID() ); } + + uint_t getNumCellsOnRank( uint_t rank ) const; + uint_t getNumFacesOnRank( uint_t rank ) const; + uint_t getNumEdgesOnRank( uint_t rank ) const; + uint_t getNumVerticesOnRank( uint_t rank ) const; + + void setGeometryMap( const PrimitiveID& primitiveID, const std::shared_ptr< GeometryMap >& map ) + { + getPrimitive( primitiveID )->geometryMap_ = map; + } + + void setMeshBoundaryFlag( const PrimitiveID& primitiveID, const uint_t& meshBoundaryFlag ) + { + getPrimitive( primitiveID )->meshBoundaryFlag_ = meshBoundaryFlag; + } + + /// Sets the mesh boundary flag of the primitives to a specified value if they are located on the boundary of the domain + /// \param meshBoundaryFlagOnBoundary the flag the primitives are set to if they are located on the boundary + /// \param meshBoundaryFlagInner the flag the primitives are set to if they are not located on the boundary + /// \param highestDimensionAlwaysInner if true, cells in 3D meshes and faces in 2D meshes are treated as inner primitives + void setMeshBoundaryFlagsOnBoundary( const uint_t& meshBoundaryFlagOnBoundary, + const uint_t& meshBoundaryFlagInner, + const bool& highestDimensionAlwaysInner ); + + /// Sets the mesh boundary flag of the primitives to a specified value if they are NOT located on the boundary of the domain. + /// This does not change mesh boundary flags on primitives that are not located on the boundary + /// \param meshBoundaryFlagInner the flag the primitives are set to if they are not located on the boundary + /// \param highestDimensionAlwaysInner if true, cells in 3D meshes and faces in 2D meshes are treated as inner primitives + void setMeshBoundaryFlagsInner( const uint_t& meshBoundaryFlagInner, const bool& highestDimensionAlwaysInner ); + + /// Every primitive for which onBoundary() returns true for all / any (if allVertices == true / false) of the primitve's vertices is assigned the passed mesh boundary flag. + void setMeshBoundaryFlagsByVertexLocation( const uint_t& meshBoundaryFlag, + const std::function< bool( const Point3D& x ) >& onBoundary, + const bool& allVertices = true ); + + /// Every primitive for which onBoundary() returns true for the primitive's centroid is assigned the passed mesh boundary flag. + void setMeshBoundaryFlagsByCentroidLocation( const uint_t& meshBoundaryFlag, + const std::function< bool( const Point3D& x ) >& onBoundary, + bool useGeometryMap = true ); + + /// Returns true, if the primitive lies on the boundary. + /// \param primitiveID the ID of the primitive to be tested + /// \param highestDimensionAlwaysInner if true, this method always returns false if the targeted primitive is of highest dimension (cells for 3D, faces for 2D meshes) + bool onBoundary( const PrimitiveID& primitiveID, const bool& highestDimensionAlwaysInner = false ) const; + + private: + typedef std::map< uint_t, std::vector< PrimitiveID::IDType > > RankToSetupPrimitivesMap; + + PrimitiveID generatePrimitiveID() const; + + void assembleRankToSetupPrimitivesMap( RankToSetupPrimitivesMap& rankToSetupPrimitivesMap ) const; + + /// Returns the number of primitives on the target rank with the least number of primitives + uint_t getMinPrimitivesPerRank() const; + /// Returns the number of primitives on the target rank with the largest number of primitives + uint_t getMaxPrimitivesPerRank() const; + /// Returns the average number of primitives per rank + real_t getAvgPrimitivesPerRank() const; + + /// Returns the number of primitives that lie on the boundary + uint_t getNumVerticesOnBoundary() const; + uint_t getNumEdgesOnBoundary() const; + uint_t getNumFacesOnBoundary() const; + uint_t getNumCellsOnBoundary() const; + + uint_t numberOfProcesses_; + + VertexMap vertices_; + EdgeMap edges_; + FaceMap faces_; + CellMap cells_; + + std::map< PrimitiveID::IDType, uint_t > primitiveIDToTargetRankMap_; }; -inline std::ostream & operator<<( std::ostream & os, const SetupPrimitiveStorage & storage ) +inline std::ostream& operator<<( std::ostream& os, const SetupPrimitiveStorage& storage ) { - storage.toStream( os ); - return os; + storage.toStream( os ); + return os; } } // namespace hyteg diff --git a/tests/hyteg/boundaryUID/BoundaryUIDTest.cpp b/tests/hyteg/boundaryUID/BoundaryUIDTest.cpp index 187587b54e8b7779382c5c8de1068cb8751eb56d..5f84441da52d773f107d01e66cb12182ecfe815a 100644 --- a/tests/hyteg/boundaryUID/BoundaryUIDTest.cpp +++ b/tests/hyteg/boundaryUID/BoundaryUIDTest.cpp @@ -24,6 +24,7 @@ #include "hyteg/boundary/BoundaryConditions.hpp" #include "hyteg/dataexport/VTKOutput.hpp" #include "hyteg/geometry/AnnulusMap.hpp" +#include "hyteg/geometry/PolarCoordsMap.hpp" #include "hyteg/mesh/MeshInfo.hpp" #include "hyteg/primitivestorage/PrimitiveStorage.hpp" #include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp" @@ -40,10 +41,12 @@ using namespace hyteg; // Dirichlet BC values. template < typename func_t > -void runTest() +void runTest( bool useCentroids ) { WALBERLA_LOG_INFO_ON_ROOT( "Running BoundaryUIDTest for " << FunctionTrait< func_t >::getTypeName() ); + bool beVerbose = false; + // ----------------------------------------- // Define markers for geometric boundaries // ----------------------------------------- @@ -57,9 +60,13 @@ void runTest() real_t outerRad = 2.0; uint_t nLayers = 2; real_t boundaryRad = 0.0; - real_t tol = real_c( 0.5 ) * ( outerRad - innerRad ) / real_c( nLayers ); + real_t tol = real_c( 0.1 ) * ( outerRad - innerRad ) / real_c( nLayers ); - MeshInfo meshInfo = MeshInfo::meshAnnulus( innerRad, outerRad, 0.0, 2.0 * pi, MeshInfo::CROSS, 8, nLayers ); + MeshInfo meshInfo = MeshInfo::meshAnnulus( innerRad, outerRad, MeshInfo::CROSS, 8, nLayers ); + + // generate the setupStorage and associate blending map + SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); + AnnulusMap::setMap( setupStorage ); // flag the inner and outer boundary by assigning different values auto onBoundary = [&boundaryRad, tol]( const Point3D& x ) { @@ -67,18 +74,33 @@ void runTest() return std::abs( boundaryRad - radius ) < tol; }; - boundaryRad = outerRad; - meshInfo.setMeshBoundaryFlagsByVertexLocation( markerOuterBoundary, onBoundary, true ); + if ( !useCentroids ) + { + boundaryRad = outerRad; + setupStorage.setMeshBoundaryFlagsByVertexLocation( markerOuterBoundary, onBoundary, true ); - boundaryRad = innerRad; - meshInfo.setMeshBoundaryFlagsByVertexLocation( markerInnerBoundary, onBoundary, true ); + boundaryRad = innerRad; + setupStorage.setMeshBoundaryFlagsByVertexLocation( markerInnerBoundary, onBoundary, true ); + } + else + { + boundaryRad = outerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerOuterBoundary, onBoundary, true ); - // generate the setupStorage and associate blending map - SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); - AnnulusMap::setMap( setupStorage ); + boundaryRad = innerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerInnerBoundary, onBoundary, true ); + } auto storage = std::make_shared< PrimitiveStorage >( setupStorage ); + // report primitives and their flags + if ( beVerbose ) + { + std::stringstream sStr; + setupStorage.toStream( sStr, true ); + WALBERLA_LOG_INFO_ON_ROOT( "" << sStr.str() ); + } + // ----------------------- // Function Manipulation // ----------------------- @@ -172,7 +194,6 @@ void runTest() // ------------ // Output VTK // ------------ - bool beVerbose = false; if ( beVerbose ) { std::string fPath = "../../output"; @@ -187,6 +208,196 @@ void runTest() } } +void captureTheFlags( bool useCentroids ) +{ + WALBERLA_LOG_INFO_ON_ROOT( "Playing 'capture the flags'" ); + + // ----------------------------------------- + // Define markers for geometric boundaries + // ----------------------------------------- + uint_t markerInnerBoundary = 11; + uint_t markerOuterBoundary = 22; + + // -------------- + // SetupStorage + // -------------- + real_t innerRad = 1.0; + real_t outerRad = 2.0; + uint_t nLayers = 2; + real_t boundaryRad = 0.0; + real_t tol = real_c( 0.1 ) * ( outerRad - innerRad ) / real_c( nLayers ); + + MeshInfo meshInfo = MeshInfo::meshAnnulus( innerRad, outerRad, MeshInfo::CROSS, 8, nLayers ); + + // generate the setupStorage and associate blending map + SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); + AnnulusMap::setMap( setupStorage ); + + // flag the inner and outer boundary by assigning different values + auto onBoundary = [&boundaryRad, tol]( const Point3D& x ) { + real_t radius = std::sqrt( x[0] * x[0] + x[1] * x[1] ); + return std::abs( boundaryRad - radius ) < tol; + }; + + if ( !useCentroids ) + { + boundaryRad = outerRad; + setupStorage.setMeshBoundaryFlagsByVertexLocation( markerOuterBoundary, onBoundary, true ); + + boundaryRad = innerRad; + setupStorage.setMeshBoundaryFlagsByVertexLocation( markerInnerBoundary, onBoundary, true ); + } + else + { + boundaryRad = outerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerOuterBoundary, onBoundary, true ); + + boundaryRad = innerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerInnerBoundary, onBoundary, true ); + } + + std::stringstream sStr; + setupStorage.toStream( sStr, true ); + WALBERLA_LOG_INFO_ON_ROOT( "Here we go:\n" << sStr.str() ); +} + +// we mesh a rectangle and use a PolarCoordsMap to map it to a half annulus +template < typename func_t > +void centroidHardBlendingTest() +{ + WALBERLA_LOG_INFO_ON_ROOT( "Running BoundaryUIDTest::centroidHardBlendingTest for " + << FunctionTrait< func_t >::getTypeName() ); + + bool beVerbose = false; + + // ----------------------------------------- + // Define markers for geometric boundaries + // ----------------------------------------- + uint_t markerInnerBoundary = 33; + uint_t markerOuterBoundary = 66; + uint_t markerSideBoundary = 99; + + // -------------- + // SetupStorage + // -------------- + real_t innerRad = real_c( 1 ); + real_t outerRad = real_c( 2 ); + real_t phiMin = real_c( 0 ); + real_t phiMax = real_c( pi ); + real_t boundaryRad = 0.0; + real_t tol = real_c( 1e-5 ); + + MeshInfo meshInfo = + MeshInfo::meshRectangle( Point2D( {innerRad, phiMin} ), Point2D( {outerRad, phiMax} ), MeshInfo::CRISS, 3, 4 ); + // MeshInfo::meshRectangle( Point2D( {innerRad, phiMin} ), Point2D( {outerRad, phiMax} ), MeshInfo::DIAMOND, 3, 2 ); + SetupPrimitiveStorage setupStorage( meshInfo, uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); + PolarCoordsMap::setMap( setupStorage ); + + // --------------- + // BoundaryFlags + // --------------- + + // oracle for curved boundaries + auto onCurvedBoundary = [&boundaryRad, tol]( const Point3D& x ) { + real_t radius = std::sqrt( x[0] * x[0] + x[1] * x[1] ); + return std::abs( boundaryRad - radius ) < tol; + }; + + // oracle for straight boundaries + auto onStraightBoundary = [tol]( const Point3D& x ) { return std::abs( x[1] ) < tol; }; + + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerSideBoundary, onStraightBoundary, true ); + + boundaryRad = outerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerOuterBoundary, onCurvedBoundary, true ); + + boundaryRad = innerRad; + setupStorage.setMeshBoundaryFlagsByCentroidLocation( markerInnerBoundary, onCurvedBoundary, true ); + + auto storage = std::make_shared< PrimitiveStorage >( setupStorage ); + + // report primitives and their flags + if ( beVerbose ) + { + std::stringstream sStr; + setupStorage.toStream( sStr, true ); + WALBERLA_LOG_INFO_ON_ROOT( "" << sStr.str() ); + } + + // ----------------------- + // Function Manipulation + // ----------------------- + uint_t minLevel = 2; + uint_t maxLevel = 2; + func_t test( "Test Func", storage, minLevel, maxLevel ); + func_t ctrl( "Ctrl Func", storage, minLevel, maxLevel ); + + // generate bc object and set different conditions on inner, outer, and straight boundaries + BoundaryCondition bcs; + BoundaryUID outerBC = bcs.createDirichletBC( "Dirichlet on outer radius", markerOuterBoundary ); + BoundaryUID innerBC = bcs.createDirichletBC( "Dirichlet on inner radius", markerInnerBoundary ); + BoundaryUID sideBC = bcs.createDirichletBC( "Dirichlet on side boundary radius", markerSideBoundary ); + + test.setBoundaryCondition( bcs ); + ctrl.setBoundaryCondition( bcs ); + + // assign functions values + real_t iValue = real_c( 30 ); // on inner boundary + real_t mValue = real_c( 20 ); // in the interior + real_t oValue = real_c( 10 ); // on outer boundary + real_t sValue = real_c( -4 ); // on side boundary + + test.interpolate( mValue, maxLevel, All ); + test.interpolate( iValue, maxLevel, innerBC ); + test.interpolate( oValue, maxLevel, outerBC ); + test.interpolate( sValue, maxLevel, sideBC ); + + // ------------------ + // Control Function + // ------------------ + std::function< real_t( const Point3D& ) > controlValues = + [innerRad, outerRad, iValue, mValue, oValue, sValue]( const Point3D& x ) { + real_t radius = std::sqrt( x[0] * x[0] + x[1] * x[1] ); + real_t mytol = 1e-14; + if ( std::abs( innerRad - radius ) < mytol ) + { + return iValue; + } + else if ( std::abs( outerRad - radius ) < mytol ) + { + return oValue; + } + else if ( std::abs( x[1] ) < mytol ) + { + return sValue; + } + return mValue; + }; + ctrl.interpolate( controlValues, maxLevel, All ); + + // ------------ + // Output VTK + // ------------ + if ( beVerbose ) + { + std::string fPath = "../../output"; + std::string fName = "centroidTest"; + WALBERLA_LOG_INFO_ON_ROOT( "Exporting to '" << fPath << "/" << fName << "'" ); + VTKOutput vtkOutput( fPath, fName, storage ); + vtkOutput.add( test ); + vtkOutput.add( ctrl ); + vtkOutput.write( maxLevel ); + } + + // ----------------------- + // check for differences + // ----------------------- + func_t diff( "Diff Func", storage, minLevel, maxLevel ); + diff.assign( {-1.0, 1.0}, {test, ctrl}, maxLevel, All ); + real_t check = diff.getMaxMagnitude( maxLevel, All ); + WALBERLA_CHECK_FLOAT_EQUAL( check, real_c( 0 ) ); +} + int main( int argc, char* argv[] ) { // ------------- @@ -196,13 +407,33 @@ int main( int argc, char* argv[] ) walberla::logging::Logging::instance()->setLogLevel( walberla::logging::Logging::PROGRESS ); walberla::MPIManager::instance()->useWorldComm(); - // ------------------------------------------ - // Run test with different function classes - // ------------------------------------------ - runTest< P1Function< real_t > >(); - runTest< P2Function< real_t > >(); - runTest< P1VectorFunction< real_t > >(); - runTest< P2VectorFunction< real_t > >(); + // ----------------------------------------------------------- + // Run test with different function classes and flag setters + // ----------------------------------------------------------- + + bool useCentroids = false; + WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------------------------" ); + WALBERLA_LOG_INFO_ON_ROOT( "setupStorage::setMeshBoundaryFlagsByVertexLocation" ); + WALBERLA_LOG_INFO_ON_ROOT( "--------------------------------------------------" ); + captureTheFlags( useCentroids ); + runTest< P1Function< real_t > >( useCentroids ); + runTest< P2Function< real_t > >( useCentroids ); + runTest< P1VectorFunction< real_t > >( useCentroids ); + runTest< P2VectorFunction< real_t > >( useCentroids ); + + useCentroids = true; + WALBERLA_LOG_INFO_ON_ROOT( "\n\n----------------------------------------------------" ); + WALBERLA_LOG_INFO_ON_ROOT( "setupStorage::setMeshBoundaryFlagsByCentroidLocation" ); + WALBERLA_LOG_INFO_ON_ROOT( "----------------------------------------------------" ); + captureTheFlags( useCentroids ); + runTest< P1Function< real_t > >( useCentroids ); + runTest< P2Function< real_t > >( useCentroids ); + runTest< P1VectorFunction< real_t > >( useCentroids ); + runTest< P2VectorFunction< real_t > >( useCentroids ); + + // use a strongly blended geometry + centroidHardBlendingTest< P1Function< real_t > >(); + centroidHardBlendingTest< P2Function< real_t > >(); return 0; } diff --git a/tutorials/09_BoundaryConditions/09_BoundaryConditions.cpp b/tutorials/09_BoundaryConditions/09_BoundaryConditions.cpp index b1c6ac6da8bede04f3239c3bf34505f7956563db..25009c504b092a33916144fffc05115eb10a35ef 100644 --- a/tutorials/09_BoundaryConditions/09_BoundaryConditions.cpp +++ b/tutorials/09_BoundaryConditions/09_BoundaryConditions.cpp @@ -89,9 +89,14 @@ * * The flags will be set on the \link hyteg::SetupPrimitiveStorage `SetupPrimitiveStorage`\endlink object we * generate from our mesh using one of the different `setMeshBoundaryFlags` methods. In our case - * \link hyteg::SetupPrimitiveStorage::setMeshBoundaryFlagsByVertexLocation() `setMeshBoundaryFlagsByVertexLocation()`\endlink - * is the appropriate one. We need to supply to it a callback that returns true or false, depending on whether a - * macro-vertex belongs to the boundary part we want to flag, or not. + * \link hyteg::SetupPrimitiveStorage::setMeshBoundaryFlagsByCentroidLocation() `setMeshBoundaryFlagsByCentroidLocation()`\endlink + * is the appropriate one. + * We need to supply to it a callback that returns true or false, depending on whether the centroid of a macro primitive + * belongs to the boundary part we want to flag, or not. + * Note that we might also use + * \link hyteg::SetupPrimitiveStorage::setMeshBoundaryFlagsByVertexLocation() `setMeshBoundaryFlagsByVertexLocation()`\endlink. + * However, this method requires a little more consideration w.r.t. the callbacks as the flags for higher-dimensional + * primitives are derived from those of their associated vertices. * * \snippet tutorials/09_BoundaryConditions/09_BoundaryConditions.cpp Flag_Boundaries_Parts * @@ -196,10 +201,10 @@ */ #include "core/DataTypes.h" +#include "core/Environment.h" #include "core/config/Config.h" #include "core/math/Constants.h" #include "core/mpi/MPIManager.h" -#include "core/Environment.h" #include "hyteg/composites/P2P1TaylorHoodFunction.hpp" #include "hyteg/dataexport/VTKOutput.hpp" @@ -215,7 +220,8 @@ using namespace hyteg; int main( int argc, char** argv ) { - walberla::Environment env( argc, argv ); // TODO: this needs to be added for MPI builds I think (did not build on my machine ;)) + walberla::Environment env( argc, + argv ); // TODO: this needs to be added for MPI builds I think (did not build on my machine ;)) walberla::MPIManager::instance()->useWorldComm(); const uint_t nx = 3; @@ -225,7 +231,7 @@ int main( int argc, char** argv ) /// [Flag_Boundaries_Parts] // generate mesh and setup storage - auto meshInfo = MeshInfo::meshRectangle( Point2D( { 0, 0 } ), Point2D( { 1.5, 1 } ), MeshInfo::CROSS, nx, ny ); + auto meshInfo = MeshInfo::meshRectangle( Point2D( {0, 0} ), Point2D( {1.5, 1} ), MeshInfo::CROSS, nx, ny ); auto setupStorage = std::make_shared< SetupPrimitiveStorage >( meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) ); @@ -237,20 +243,24 @@ int main( int argc, char** argv ) setupStorage->setMeshBoundaryFlagsOnBoundary( 1, 0, true ); // callbacks for marking parts of the boundary - real_t eps = 1e-3; - auto topBoundary = [ eps ]( const Point3D& x ) { return std::abs( x[1] - real_c( 1 ) ) < eps; }; + real_t eps = 1e-3; + auto topBoundary = [eps]( const Point3D& x ) { return std::abs( x[1] - real_c( 1 ) ) < eps; }; - auto bottomBoundary = [ eps ]( const Point3D& x ) { return x[1] < eps; }; + auto bottomBoundary = [eps]( const Point3D& x ) { return x[1] < eps; }; - auto leftBoundary = [ eps ]( const Point3D& x ) { return x[0] < eps; }; + auto leftBoundary = [eps]( const Point3D& x ) { + return x[0] < eps && !( std::abs( x[1] - real_c( 1 ) ) < eps ) && !( x[1] < eps ); + }; - auto rightBoundary = [ eps ]( const Point3D& x ) { return std::abs( x[0] - real_c( 1.5 ) ) < eps; }; + auto rightBoundary = [eps]( const Point3D& x ) { + return std::abs( x[0] - real_c( 1.5 ) ) < eps && !( std::abs( x[1] - real_c( 1 ) ) < eps ) && !( x[1] < eps ); + }; // assigning mesh boundary flags to different parts of the boundary - setupStorage->setMeshBoundaryFlagsByVertexLocation( 1, topBoundary ); - setupStorage->setMeshBoundaryFlagsByVertexLocation( 2, bottomBoundary ); - setupStorage->setMeshBoundaryFlagsByVertexLocation( 3, leftBoundary ); - setupStorage->setMeshBoundaryFlagsByVertexLocation( 4, rightBoundary ); + setupStorage->setMeshBoundaryFlagsByCentroidLocation( 1, topBoundary ); + setupStorage->setMeshBoundaryFlagsByCentroidLocation( 2, bottomBoundary ); + setupStorage->setMeshBoundaryFlagsByCentroidLocation( 3, leftBoundary ); + setupStorage->setMeshBoundaryFlagsByCentroidLocation( 4, rightBoundary ); // now create the actual storage auto storage = std::make_shared< PrimitiveStorage >( *setupStorage ); @@ -260,13 +270,13 @@ int main( int argc, char** argv ) BoundaryCondition bcVelocity; BoundaryCondition bcTemperature; - BoundaryUID idTopWallVel = bcVelocity.createDirichletBC( "topWall", { 1 } ); - BoundaryUID idBotWallVel = bcVelocity.createDirichletBC( "botWall", { 2 } ); - bcVelocity.createFreeslipBC( "sideWalls", { 3, 4 } ); + BoundaryUID idTopWallVel = bcVelocity.createDirichletBC( "topWall", {1} ); + BoundaryUID idBotWallVel = bcVelocity.createDirichletBC( "botWall", {2} ); + bcVelocity.createFreeslipBC( "sideWalls", {3, 4} ); - BoundaryUID idTopWallTemp = bcTemperature.createDirichletBC( "dirichletWallTop", { 1 } ); - BoundaryUID idBotWallTemp = bcTemperature.createDirichletBC( "dirichletWallBot", { 2 } ); - bcTemperature.createNeumannBC( "neumannWalls", { 3, 4 } ); + BoundaryUID idTopWallTemp = bcTemperature.createDirichletBC( "dirichletWallTop", {1} ); + BoundaryUID idBotWallTemp = bcTemperature.createDirichletBC( "dirichletWallBot", {2} ); + bcTemperature.createNeumannBC( "neumannWalls", {3, 4} ); /// [BC_Objects] /// [Function_Creation] @@ -282,8 +292,8 @@ int main( int argc, char** argv ) /// [Setting_BC_Temp] // fixed temperatures on top and bottom edge - real_t topT = real_c( 0.0 ); - real_t botT = real_c( 1.0 ); + real_t topT = real_c( 0.0 ); + real_t botT = real_c( 1.0 ); // as the values are a constant we do not need to provide a callback function to the interpolation method temperature.interpolate( botT, maxLevel, idBotWallTemp ); @@ -296,20 +306,20 @@ int main( int argc, char** argv ) // per component (note that we cannot mix the types, though) // set velocity field on bottom wall to no-slip - velocity.interpolate( { real_c( 0 ), real_c( 0 ) }, maxLevel, idBotWallVel ); + velocity.interpolate( {real_c( 0 ), real_c( 0 )}, maxLevel, idBotWallVel ); // set velocity field on top wall to no-outflow + non-zero tangential component - velocity.interpolate( { real_c( 2 ), real_c( 0 ) }, maxLevel, idTopWallVel ); + velocity.interpolate( {real_c( 2 ), real_c( 0 )}, maxLevel, idTopWallVel ); // to set the boundary values on our TaylorHood funcion, we extract its velocity component - uAndp.uvw.interpolate( { real_c( 0 ), real_c( 0 ) }, maxLevel, idBotWallVel ); - uAndp.uvw.interpolate( { real_c( 2 ), real_c( 0 ) }, maxLevel, idTopWallVel ); + uAndp.uvw.interpolate( {real_c( 0 ), real_c( 0 )}, maxLevel, idBotWallVel ); + uAndp.uvw.interpolate( {real_c( 2 ), real_c( 0 )}, maxLevel, idTopWallVel ); /// [Setting_BC_Vel] /// [Setting_ICs] auto initialTemperature = []( const Point3D& x ) { - return ( 1.0 - x[1] ) + 0.01 * std::cos( pi * x[0] / real_c(1.5) ) * std::sin( pi * x[1] ); + return ( 1.0 - x[1] ) + 0.01 * std::cos( pi * x[0] / real_c( 1.5 ) ) * std::sin( pi * x[1] ); }; temperature.interpolate( initialTemperature, maxLevel, Inner | NeumannBoundary | FreeslipBoundary ); /// [Setting_ICs]