diff --git a/doc/grids/gmsh/junction2d3dbseg.geo b/doc/grids/gmsh/junction2d3dbseg.geo new file mode 100644 index 0000000000000000000000000000000000000000..b5c5cc00546efc23db99e2e2bf24dd1aad181dac --- /dev/null +++ b/doc/grids/gmsh/junction2d3dbseg.geo @@ -0,0 +1,22 @@ +Point(1) = { 0, 0, 0, 10.0}; +Point(2) = { 0, 0, 1, 10.0}; +Point(3) = {-1, 0, 0, 10.0}; +Point(4) = { 1, 0, 1, 10.0}; +Point(5) = { 1, 1, 1, 10.0}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 1}; +Line(4) = {1, 4}; +Line(5) = {4, 2}; +Line(6) = {2, 5}; +Line(7) = {5, 1}; + +Line Loop(1) = {1, 2, 3}; +Plane Surface(1) = {1}; +Line Loop(2) = {-1, 4, 5}; +Plane Surface(2) = {2}; +Line Loop(3) = {1, 6, 7}; +Plane Surface(3) = {3}; +Physical Surface(1) = {1:3}; + +Physical Line(42) = {2,3,4}; diff --git a/doc/grids/gmsh/junction2d3dbseg.msh b/doc/grids/gmsh/junction2d3dbseg.msh new file mode 100644 index 0000000000000000000000000000000000000000..d685615a38c142718e4c0b8b1ab7eb714de84677 --- /dev/null +++ b/doc/grids/gmsh/junction2d3dbseg.msh @@ -0,0 +1,20 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +5 +1 0 0 0 +2 0 0 1 +3 -1 0 0 +4 1 0 1 +5 1 1 1 +$EndNodes +$Elements +6 +1 1 2 42 2 2 3 +2 1 2 42 3 3 1 +3 1 2 42 4 1 4 +4 2 2 1 1 1 2 3 +5 2 2 1 2 1 4 2 +6 2 2 1 3 1 2 5 +$EndElements diff --git a/doc/grids/gmsh/line1d2dbseg.geo b/doc/grids/gmsh/line1d2dbseg.geo new file mode 100644 index 0000000000000000000000000000000000000000..d8fd97d98b0710dc54e9affcc2b4d56d1f091f5e --- /dev/null +++ b/doc/grids/gmsh/line1d2dbseg.geo @@ -0,0 +1,12 @@ +Point(1) = {0, 0, 0, 1.0}; +Point(2) = {1, 1, 0, 1.0}; +Point(3) = {-1, 0, 0, 1.0}; +Point(4) = {-2, 0, 0, 1.0}; +Point(5) = {-3, 1, 0, 1.0}; +Line(1) = {3, 1}; +Line(2) = {1, 2}; +Line(3) = {3, 4}; +Line(4) = {4, 5}; + +Physical Line(1) = {1:4}; +Physical Point(42) = {2}; diff --git a/doc/grids/gmsh/line1d2dbseg.msh b/doc/grids/gmsh/line1d2dbseg.msh new file mode 100644 index 0000000000000000000000000000000000000000..7c2d67a7956a8792906064cad0b6350b2042feda --- /dev/null +++ b/doc/grids/gmsh/line1d2dbseg.msh @@ -0,0 +1,19 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +5 +1 0 0 0 +2 1 1 0 +3 -1 0 0 +4 -2 0 0 +5 -3 1 0 +$EndNodes +$Elements +5 +1 15 2 42 2 2 +2 1 2 1 1 3 1 +3 1 2 1 2 1 2 +4 1 2 1 3 3 4 +5 1 2 1 4 4 5 +$EndElements diff --git a/dune/foamgrid/foamgrid/foamgridfactory.hh b/dune/foamgrid/foamgrid/foamgridfactory.hh index 2787f55a2be7b5bad4e046604ddc9a7d185b2729..1d0e9c700f27a01c62745acc28d7efe29de04655 100644 --- a/dune/foamgrid/foamgrid/foamgridfactory.hh +++ b/dune/foamgrid/foamgrid/foamgridfactory.hh @@ -9,8 +9,8 @@ */ #include -#include #include +#include #include #include @@ -75,25 +75,6 @@ template vertexArray_.push_back(&*std::get<0>(grid_->entityImps_[0]).rbegin()); } - /** \brief Insert a boundary segment. - - This is only needed if you want to control the numbering of the boundary segments - */ - void insertBoundarySegment(const std::vector& vertices) override { - DUNE_THROW(Dune::NotImplemented, "insertBoundarySegment not implemented yet!"); - } - - /** \brief Insert a boundary segment (== a line) and the boundary segment geometry - * - This influences the ordering of the boundary segments. - Currently, the BoundarySegment object does not actually have any effect. - */ - void insertBoundarySegment(const std::vector& vertices, - const std::shared_ptr >& boundarySegment) override - { - insertBoundarySegment(vertices); - } - /** \brief Return the number of the element to insert parameters */ unsigned int @@ -121,6 +102,9 @@ template /** \brief Array containing all vertices */ std::vector*> vertexArray_; + + /** \brief Counter that creates the boundary segment indices */ + unsigned int boundarySegmentCounter_ = 0; }; template @@ -147,6 +131,35 @@ template GridFactoryBase<1,dimworld>(grid) {} + /** \brief Insert a boundary segment. + This is only needed if you want to control the numbering of the boundary segments + */ + void insertBoundarySegment(const std::vector& vertices) override + { + boundarySegmentIndices_[ vertices[0] ] = this->boundarySegmentCounter_++; + } + + /** \brief Insert a boundary segment and the boundary segment geometry + This influences the ordering of the boundary segments. + */ + void insertBoundarySegment(const std::vector& vertices, + const std::shared_ptr >& boundarySegment) override + { + DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented"); + } + + /** \brief Return true if leaf intersection was inserted as boundary segment + */ + bool wasInserted( const typename FoamGrid::LeafIntersection &intersection ) const override + { + if ( !intersection.boundary() || intersection.inside().level() != 0 ) + return false; + + const auto& vertex = intersection.inside().template subEntity<1>(intersection.indexInInside()); + const auto& it = boundarySegmentIndices_.find( this->insertionIndex(vertex) ); + return (it != boundarySegmentIndices_.end()); + } + /** \brief Insert an element into the coarse grid \param type The GeometryType of the new element \param vertices The vertices of the new element, using the DUNE numbering @@ -221,25 +234,31 @@ template // Set the boundary ids // //////////////////////////////////////////////// - unsigned int boundaryFacetCounter = 0; - // Iterate over all facets (=vertices in 1d) FacetIterator fIt = std::get<0>(this->grid_->entityImps_[0]).begin(); const FacetIterator fEndIt = std::get<0>(this->grid_->entityImps_[0]).end(); for (; fIt != fEndIt; ++fIt) if(fIt->elements_.size()==1) // if boundary facet - fIt->boundarySegmentIndex_ = boundaryFacetCounter++; + { + const auto& it = boundarySegmentIndices_.find( fIt->vertex_[0]->leafIndex_ ); + if (it != boundarySegmentIndices_.end()) + fIt->boundarySegmentIndex_ = it->second; + else + fIt->boundarySegmentIndex_ = this->boundarySegmentCounter_++; + } // //////////////////////////////////////////////// // Hand over the new grid // //////////////////////////////////////////////// Dune::FoamGrid* tmp = this->grid_; - tmp->numBoundarySegments_ = boundaryFacetCounter; + tmp->numBoundarySegments_ = this->boundarySegmentCounter_; this->grid_ = nullptr; return tmp; } + private: + std::unordered_map boundarySegmentIndices_; }; /** \brief Specialization of GridFactoryBase for 2D-FoamGrid<2, dimworld> @@ -261,6 +280,57 @@ template GridFactoryBase<2,dimworld>(grid) {} + /** \brief Insert a boundary segment. + This is only needed if you want to control the numbering of the boundary segments + */ + void insertBoundarySegment(const std::vector& vertices) override + { + std::array vertexIndices {{ vertices[0], vertices[1] }}; + + // sort the indices + if ( vertexIndices[0] > vertexIndices[1] ) + std::swap( vertexIndices[0], vertexIndices[1] ); + + boundarySegmentIndices_[ vertexIndices ] = this->boundarySegmentCounter_++; + } + + /** \brief Insert a boundary segment (== a line) and the boundary segment geometry + This influences the ordering of the boundary segments. + */ + void insertBoundarySegment(const std::vector& vertices, + const std::shared_ptr >& boundarySegment) override + { + DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented"); + } + + /** \brief Return true if leaf intersection was inserted as boundary segment + */ + bool wasInserted( const typename FoamGrid::LeafIntersection &intersection ) const override + { + if ( !intersection.boundary() || intersection.inside().level() != 0 ) + return false; + + // Get the vertices of the intersection + const auto refElement = ReferenceElements::general(intersection.inside().type()); + + const int subIdx0 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/0, dimgrid); + const auto vertex0 = intersection.inside().template subEntity<2>( subIdx0 ); + const int subIdx1 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/1, dimgrid); + const auto vertex1 = intersection.inside().template subEntity<2>( subIdx1 ); + + std::array vertexIndices {{ + this->insertionIndex( vertex0 ), + this->insertionIndex( vertex1 ) + }}; + + // sort the indices + if ( vertexIndices[0] > vertexIndices[1] ) + std::swap( vertexIndices[0], vertexIndices[1] ); + + const auto& it = boundarySegmentIndices_.find( vertexIndices ); + return (it != boundarySegmentIndices_.end()); + } + /** \brief Insert an element into the coarse grid \param type The GeometryType of the new element \param vertices The vertices of the new element, using the DUNE numbering @@ -381,25 +451,44 @@ template // Set the boundary ids // //////////////////////////////////////////////// - unsigned int boundaryFacetCounter = 0; - // Iterate over all facets (=edges in 2D) FacetIterator fIt = std::get<1>(this->grid_->entityImps_[0]).begin(); const FacetIterator fEndIt = std::get<1>(this->grid_->entityImps_[0]).end(); for (; fIt!=fEndIt; ++fIt) if(fIt->elements_.size()==1) //if boundary facet - fIt->boundarySegmentIndex_ = boundaryFacetCounter++; - + { + std::array vertexIndices {{ fIt->vertex_[0]->leafIndex_, fIt->vertex_[1]->leafIndex_ }}; + + // sort the indices + if ( vertexIndices[0] > vertexIndices[1] ) + std::swap( vertexIndices[0], vertexIndices[1] ); + + auto it = boundarySegmentIndices_.find( vertexIndices ); + if (it != boundarySegmentIndices_.end()) { + fIt->boundarySegmentIndex_ = it->second; + } else { // edge was not inserted as boundary segment + fIt->boundarySegmentIndex_ = this->boundarySegmentCounter_++; + } + } // //////////////////////////////////////////////// // Hand over the new grid // //////////////////////////////////////////////// Dune::FoamGrid* tmp = this->grid_; - tmp->numBoundarySegments_ = boundaryFacetCounter; + tmp->numBoundarySegments_ = this->boundarySegmentCounter_; this->grid_ = nullptr; return tmp; } + + private: + struct HashUIntArray { + std::size_t operator() (const std::array& a) const { + return std::hash{}(a[0]) ^ (std::hash{}(a[1]) << 1); + } + }; + + std::unordered_map, unsigned int, HashUIntArray> boundarySegmentIndices_; }; } // end namespace Dune diff --git a/dune/foamgrid/test/CMakeLists.txt b/dune/foamgrid/test/CMakeLists.txt index 2bdeaccb5975c364d886a6c21c7b24d593fc5856..71f74ca15435e4a7f2385daa9821a2c933a05d13 100644 --- a/dune/foamgrid/test/CMakeLists.txt +++ b/dune/foamgrid/test/CMakeLists.txt @@ -6,9 +6,14 @@ if(${DUNE_GRID_VERSION} VERSION_LESS 2.7) # in the grid test about the sum over all outer normals of a grid view dune_add_test(SOURCES foamgrid-test.cc EXPECT_FAIL) dune_add_test(SOURCES local-refine-test.cc EXPECT_FAIL) + # before version 2.7 this test fails to run because of a bug + # in the GmshReader that failed to read boundary segments for 1d grids + # (fixed in dune-grid commit eaadf1b1) + dune_add_test(SOURCES boundary-segment-test.cc COMPILE_ONLY) else() dune_add_test(SOURCES foamgrid-test.cc) dune_add_test(SOURCES local-refine-test.cc) + dune_add_test(SOURCES boundary-segment-test.cc) endif() dune_add_test(SOURCES global-refine-test.cc) diff --git a/dune/foamgrid/test/boundary-segment-test.cc b/dune/foamgrid/test/boundary-segment-test.cc new file mode 100644 index 0000000000000000000000000000000000000000..877b0c4cc2ba201a3ea7d07a3c483576e652c91c --- /dev/null +++ b/dune/foamgrid/test/boundary-segment-test.cc @@ -0,0 +1,107 @@ +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +template +void checkBoundarySegments(const Grid& grid, + const Dune::GridFactory& factory, + const std::vector& boundaryMarkers, + int boundaryIntersectionsExpected, + int insertedBoundariesExpected, + int numIntersectionsExpected, + int boundaryIdExpected = 42) +{ + std::cout << " Checking boundary segment indices" << std::endl; + + int numIntersections = 0; + int boundaryIntersections = 0; + int insertedBoundaries = 0; + + for ( const auto& element : elements(grid.leafGridView()) ) + { + for ( const auto& intersection : intersections(grid.leafGridView(), element) ) + { + numIntersections++; + if ( intersection.boundary() ) + { + boundaryIntersections++; + if ( factory.wasInserted( intersection ) ) + { + insertedBoundaries++; + const int bndsegIdx = intersection.boundarySegmentIndex(); + const int bndId = boundaryMarkers[bndsegIdx]; + if (bndId != boundaryIdExpected) + DUNE_THROW(Dune::GridError, "Wrong boundary ID. Expected " << boundaryIdExpected << " got " << bndId); + } + } + } + } + + // sanity checks + if (boundaryIntersections != boundaryIntersectionsExpected) + DUNE_THROW(Dune::GridError, "Wrong number of boundary intersections. Expected " << boundaryIntersectionsExpected << " got " << boundaryIntersections); + if (insertedBoundaries != insertedBoundariesExpected) + DUNE_THROW(Dune::GridError, "Wrong number of inserted boundaries. Expected " << insertedBoundariesExpected << " got " << insertedBoundaries); + if (numIntersections != numIntersectionsExpected) + DUNE_THROW(Dune::GridError, "Wrong number of intersections. Expected " << numIntersectionsExpected << " got " << numIntersections); +} + +int main (int argc, char *argv[]) try +{ + using namespace Dune; + + Dune::MPIHelper::instance(argc, argv); + + // paths to gmsh test files + const std::string dune_foamgrid_path = std::string(DUNE_FOAMGRID_EXAMPLE_GRIDS_PATH) + "gmsh/"; + + { + std::cout << "\n################################################\n"; + std::cout << "Checking FoamGrid<1, 2> (1d in 2d grid)\n"; + std::cout << "################################################\n\n"; + + std::cout << " Creating grid" << std::endl; + using Grid = FoamGrid<1, 2>; + GridFactory factory; + + std::vector boundaryMarkers, elementMarkers; + GmshReader::read(factory, dune_foamgrid_path + "line1d2dbseg.msh", boundaryMarkers, elementMarkers, /*verbose*/ true, /*insertBoundarySegments*/ true); + auto grid = std::shared_ptr(factory.createGrid()); + + // check the boundary segments + checkBoundarySegments(*grid, factory, boundaryMarkers, 2, 1, 8); + } + + { + std::cout << "\n################################################\n"; + std::cout << "Checking FoamGrid<2, 3> (2d in 3d grid)\n"; + std::cout << "################################################\n\n"; + + std::cout << " Creating grid" << std::endl; + using Grid = FoamGrid<2, 3>; + GridFactory factory; + + std::vector boundaryMarkers, elementMarkers; + GmshReader::read(factory, dune_foamgrid_path + "junction2d3dbseg.msh", boundaryMarkers, elementMarkers, /*verbose*/ true, /*insertBoundarySegments*/ true); + auto grid = std::shared_ptr(factory.createGrid()); + + // check the boundary segments + checkBoundarySegments(*grid, factory, boundaryMarkers, 6, 3, 12); + } +} +// ////////////////////////////////// +// Error handler +// ///////////////////////////////// +catch (const Dune::Exception& e) { + std::cout << e << std::endl; + return 1; +}