8 #ifndef DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH 9 #define DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH 22 #include <dune/common/fvector.hh> 23 #include <dune/common/bitsetvector.hh> 24 #include <dune/common/stdstreams.hh> 25 #include <dune/common/timer.hh> 27 #include <dune/geometry/referenceelements.hh> 28 #include <dune/grid/common/grid.hh> 53 template<
class T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
55 :
public Merger<T,grid1Dim,grid2Dim,dimworld>
96 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
97 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
98 unsigned int grid1Index,
99 const Dune::GeometryType& grid2ElementType,
100 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
101 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
102 unsigned int grid2Index,
103 std::vector<SimplicialIntersection>& intersections) = 0;
109 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
110 const std::vector<Dune::GeometryType>& grid1_element_types,
111 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
112 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
113 const std::vector<Dune::GeometryType>& grid2_element_types,
114 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
136 void build(
const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
137 const std::vector<unsigned int>& grid1_elements,
138 const std::vector<Dune::GeometryType>& grid1_element_types,
139 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
140 const std::vector<unsigned int>& grid2_elements,
141 const std::vector<Dune::GeometryType>& grid2_element_types)
override;
149 intersectionListProvider_->clear();
150 purge(grid1ElementCorners_);
151 purge(grid2ElementCorners_);
164 m_enableFallback = fallback;
169 m_enableBruteForce = bruteForce;
176 bool m_enableFallback =
false;
181 bool m_enableBruteForce =
false;
183 auto& intersections()
184 {
return intersectionListProvider_->intersections(); }
188 static void purge(V & v)
199 void generateSeed(std::vector<int>& seeds,
200 Dune::BitSetVector<1>& isHandled2,
201 std::stack<unsigned>& candidates2,
202 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
203 const std::vector<Dune::GeometryType>& grid1_element_types,
204 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
205 const std::vector<Dune::GeometryType>& grid2_element_types);
210 int insertIntersections(
unsigned int candidate1,
unsigned int candidate2,std::vector<SimplicialIntersection>& intersections);
215 int bruteForceSearch(
int candidate1,
216 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
217 const std::vector<Dune::GeometryType>& grid1_element_types,
218 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
219 const std::vector<Dune::GeometryType>& grid2_element_types);
224 std::pair<bool, unsigned int>
225 intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
231 template <
int gr
idDim>
232 void computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
233 const std::vector<std::vector<unsigned int> >& gridElementCorners,
234 std::vector<std::vector<int> >& elementNeighbors);
236 void buildAdvancingFront(
237 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
238 const std::vector<unsigned int>& grid1_elements,
239 const std::vector<Dune::GeometryType>& grid1_element_types,
240 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
241 const std::vector<unsigned int>& grid2_elements,
242 const std::vector<Dune::GeometryType>& grid2_element_types
245 void buildBruteForce(
246 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
247 const std::vector<unsigned int>& grid1_elements,
248 const std::vector<Dune::GeometryType>& grid1_element_types,
249 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
250 const std::vector<unsigned int>& grid2_elements,
251 const std::vector<Dune::GeometryType>& grid2_element_types
258 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
260 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
261 const std::vector<Dune::GeometryType>& grid1_element_types,
262 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
263 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
264 const std::vector<Dune::GeometryType>& grid2_element_types,
265 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
270 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
271 for (
int i=0; i<grid1NumVertices; i++)
276 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
277 for (
int i=0; i<grid2NumVertices; i++)
284 std::vector<SimplicialIntersection> intersections(0);
288 neighborIntersects1, candidate0,
289 grid2_element_types[candidate1], grid2ElementCorners,
290 neighborIntersects2, candidate1,
294 if(insert && !intersections.empty())
295 insertIntersections(candidate0,candidate1,intersections);
298 return !intersections.empty() || neighborIntersects1.any() || neighborIntersects2.any();
302 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
304 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
305 const std::vector<Dune::GeometryType>& grid1_element_types,
306 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
307 const std::vector<Dune::GeometryType>& grid2_element_types)
309 std::bitset<(1<<grid1Dim)> neighborIntersects1;
310 std::bitset<(1<<grid2Dim)> neighborIntersects2;
311 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
314 grid1Coords, grid1_element_types, neighborIntersects1,
315 grid2Coords, grid2_element_types, neighborIntersects2,
319 if (intersectionFound)
328 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
329 template<
int gr
idDim>
332 const std::vector<std::vector<unsigned int> >& gridElementCorners,
333 std::vector<std::vector<int> >& elementNeighbors)
335 typedef std::vector<unsigned int> FaceType;
336 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
342 elementNeighbors.resize(gridElementTypes.size());
344 for (
size_t i=0; i<gridElementTypes.size(); i++)
345 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
347 for (
size_t i=0; i<gridElementTypes.size(); i++) {
348 const auto& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
350 for (
size_t j=0; j<(size_t)refElement.size(1); j++) {
354 for (
size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
355 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
358 std::sort(face.begin(), face.end());
360 typename FaceSetType::iterator faceHandle = faces.find(face);
362 if (faceHandle == faces.end()) {
365 faces.insert(std::make_pair(face, std::make_pair(i,j)));
370 elementNeighbors[i][j] = faceHandle->second.first;
371 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
373 faces.erase(faceHandle);
387 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
389 const std::vector<unsigned int>& grid1_elements,
390 const std::vector<Dune::GeometryType>& grid1_element_types,
391 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
392 const std::vector<unsigned int>& grid2_elements,
393 const std::vector<Dune::GeometryType>& grid2_element_types
397 std::cout <<
"StandardMerge building merged grid..." << std::endl;
414 unsigned int grid1CornerCounter = 0;
416 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
419 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
421 for (
int j=0; j<numVertices; j++)
429 unsigned int grid2CornerCounter = 0;
431 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
434 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
436 for (
int j=0; j<numVertices; j++)
448 std::cout <<
"setup took " << watch.elapsed() <<
" seconds." << std::endl;
450 if (m_enableBruteForce)
451 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
453 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
456 std::cout <<
"intersection construction took " << watch.elapsed() <<
" seconds." << std::endl;
459 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
461 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
462 const std::vector<unsigned int>& grid1_elements,
463 const std::vector<Dune::GeometryType>& grid1_element_types,
464 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
465 const std::vector<unsigned int>& grid2_elements,
466 const std::vector<Dune::GeometryType>& grid2_element_types
473 std::stack<unsigned int> candidates1;
474 std::stack<unsigned int> candidates2;
476 std::vector<int> seeds(grid2_element_types.size(), -1);
484 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
487 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
489 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
495 std::set<unsigned int> isHandled1;
497 std::set<unsigned int> isCandidate1;
499 while (!candidates2.empty()) {
502 unsigned int currentCandidate2 = candidates2.top();
503 int seed = seeds[currentCandidate2];
507 isHandled2[currentCandidate2] =
true;
511 candidates1.push(seed);
514 isCandidate1.clear();
516 while (!candidates1.empty()) {
518 unsigned int currentCandidate1 = candidates1.top();
520 isHandled1.insert(currentCandidate1);
523 std::bitset<(1<<grid1Dim)> neighborIntersects1;
524 std::bitset<(1<<grid2Dim)> neighborIntersects2;
526 grid1Coords,grid1_element_types, neighborIntersects1,
527 grid2Coords,grid2_element_types, neighborIntersects2);
529 for (
size_t i=0; i<neighborIntersects2.size(); i++)
534 if (intersectionFound) {
543 if (isHandled1.find(neighbor) == isHandled1.end()
544 && isCandidate1.find(neighbor) == isCandidate1.end()) {
545 candidates1.push(neighbor);
546 isCandidate1.insert(neighbor);
560 bool seedFound = !candidates2.empty();
569 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
571 isCandidate2[neighbor][0] =
true;
572 candidates2.push(neighbor);
577 if (seedFound || !m_enableFallback)
589 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
596 for (
typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
597 seedIt != isHandled1.end(); ++seedIt) {
599 std::bitset<(1<<grid1Dim)> neighborIntersects1;
600 std::bitset<(1<<grid2Dim)> neighborIntersects2;
602 grid1Coords, grid1_element_types, neighborIntersects1,
603 grid2Coords, grid2_element_types, neighborIntersects2,
607 if (intersectionFound) {
609 Dune::dwarn <<
"Algorithm entered first fallback method and found a new seed in the build algorithm." <<
610 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
619 seed = bruteForceSearch(neighbor,
620 grid1Coords,grid1_element_types,
621 grid2Coords,grid2_element_types);
622 Dune::dwarn <<
"Algorithm entered second fallback method. This probably should not happen." << std::endl;
627 isCandidate2[neighbor] =
true;
634 candidates2.push(neighbor);
635 seeds[neighbor] = seed;
645 if (!seedFound && candidates2.empty()) {
646 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
651 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
653 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
654 const std::vector<unsigned int>& grid1_elements,
655 const std::vector<Dune::GeometryType>& grid1_element_types,
656 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
657 const std::vector<unsigned int>& grid2_elements,
658 const std::vector<Dune::GeometryType>& grid2_element_types
661 std::bitset<(1<<grid1Dim)> neighborIntersects1;
662 std::bitset<(1<<grid2Dim)> neighborIntersects2;
664 for (
unsigned i = 0; i < grid1_element_types.size(); ++i) {
665 for (
unsigned j = 0; j < grid2_element_types.size(); ++j) {
667 grid1Coords, grid1_element_types, neighborIntersects1,
668 grid2Coords, grid2_element_types, neighborIntersects2);
673 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
674 void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::generateSeed(std::vector<int>& seeds, Dune::BitSetVector<1>& isHandled2, std::stack<unsigned>& candidates2,
const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
const std::vector<Dune::GeometryType>& grid1_element_types,
const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
const std::vector<Dune::GeometryType>& grid2_element_types)
676 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
678 if (seeds[j] > 0 || isHandled2[j][0])
681 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
688 isHandled2[j] =
true;
692 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
694 std::vector<SimplicialIntersection>& intersections)
696 typedef typename std::vector<SimplicialIntersection>::size_type size_t;
699 for (
size_t i = 0; i < intersections.size(); ++i) {
703 std::tie(found, index) = intersectionIndex(candidate1,candidate2,intersections[i]);
705 if (found && index >= this->intersections().size()) {
706 this->intersections().push_back(intersections[i]);
710 auto& intersection = this->intersections()[index];
713 for (
size_t j = 0; j < intersections[i].parents0.size(); ++j) {
714 intersection.parents0.push_back(candidate1);
715 intersection.corners0.push_back(intersections[i].corners0[j]);
719 for (
size_t j = 0; j < intersections[i].parents1.size(); ++j) {
720 intersection.parents1.push_back(candidate2);
721 intersection.corners1.push_back(intersections[i].corners1[j]);
726 Dune::dwarn <<
"Computed the same intersection twice!" << std::endl;
732 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
733 std::pair<bool, unsigned int>
741 std::size_t n_intersections = this->intersections().size();
742 if (grid1Dim == grid2Dim)
743 return {
true, n_intersections};
747 for (std::size_t i = 0; i < n_intersections; ++i) {
750 for (std::size_t ei = 0; ei < this->intersections()[i].parents0.size(); ++ei)
752 if (this->intersections()[i].parents0[ei] == grid1Index)
754 for (std::size_t er = 0; er < intersection.parents0.size(); ++er)
756 bool found_all =
true;
758 for (std::size_t ci = 0; ci < this->intersections()[i].corners0[ei].size(); ++ci)
760 Dune::FieldVector<T,grid1Dim> ni = this->intersections()[i].corners0[ei][ci];
761 bool found_ni =
false;
762 for (std::size_t cr = 0; cr < intersection.corners0[er].size(); ++cr)
764 Dune::FieldVector<T,grid1Dim> nr = intersection.corners0[er][cr];
766 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
770 found_all = found_all && found_ni;
776 if (found_all && (this->intersections()[i].parents1[ei] != grid2Index))
785 for (std::size_t ei = 0; ei < this->intersections()[i].parents1.size(); ++ei)
787 if (this->intersections()[i].parents1[ei] == grid2Index)
789 for (std::size_t er = 0; er < intersection.parents1.size(); ++er)
791 bool found_all =
true;
793 for (std::size_t ci = 0; ci < this->intersections()[i].corners1[ei].size(); ++ci)
795 Dune::FieldVector<T,grid2Dim> ni = this->intersections()[i].corners1[ei][ci];
796 bool found_ni =
false;
797 for (std::size_t cr = 0; cr < intersection.corners1[er].size(); ++cr)
799 Dune::FieldVector<T,grid2Dim> nr = intersection.corners1[er][cr];
800 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
805 found_all = found_all && found_ni;
811 if (found_all && (this->intersections()[i].parents0[ei] != grid1Index))
820 return {
true, n_intersections};
824 #define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \ 826 void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \ 827 const std::vector<unsigned int>& grid1_elements, \ 828 const std::vector<Dune::GeometryType>& grid1_element_types, \ 829 const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \ 830 const std::vector<unsigned int>& grid2_elements, \ 831 const std::vector<Dune::GeometryType>& grid2_element_types \ 837 #undef STANDARD_MERGE_INSTANTIATE 843 #endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: merger.hh:35
std::shared_ptr< IntersectionListProvider > intersectionListProvider_
Definition: standardmerge.hh:119
std::shared_ptr< IntersectionList > intersectionList_
Definition: standardmerge.hh:120
SimplicialIntersection RemoteSimplicialIntersection
Definition: standardmerge.hh:82
Definition: intersectionlist.hh:219
Dune::FieldVector< T, grid2Dim > Grid2Coords
the local coordinate type for the grid2 coordinates
Definition: merger.hh:32
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:124
void enableFallback(bool fallback)
Definition: standardmerge.hh:162
STANDARD_MERGE_INSTANTIATE(double, 1, 1, 1)
Dune::FieldVector< T, grid1Dim > Grid1Coords
the local coordinate type for the grid1 coordinates
Definition: merger.hh:29
virtual void computeIntersections(const Dune::GeometryType &grid1ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid1ElementCorners, std::bitset<(1<< grid1Dim)> &neighborIntersects1, unsigned int grid1Index, const Dune::GeometryType &grid2ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid2ElementCorners, std::bitset<(1<< grid2Dim)> &neighborIntersects2, unsigned int grid2Index, std::vector< SimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
Definition: intersectionlist.hh:203
bool valid
Definition: standardmerge.hh:84
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:24
unsigned int counter
Counts the number of times the computeIntersection method has been called.
Definition: merger.hh:112
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition: standardmerge.hh:81
StandardMerge()
Definition: standardmerge.hh:86
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:123
void build(const std::vector< Dune::FieldVector< T, dimworld > > &grid1_Coords, const std::vector< unsigned int > &grid1_elements, const std::vector< Dune::GeometryType > &grid1_element_types, const std::vector< Dune::FieldVector< T, dimworld > > &grid2_coords, const std::vector< unsigned int > &grid2_elements, const std::vector< Dune::GeometryType > &grid2_element_types) override
builds the merged grid
Definition: standardmerge.hh:388
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:54
Dune::GridGlue::IntersectionList< Grid1Coords, Grid2Coords > IntersectionList
Definition: merger.hh:37
bool computeIntersection(unsigned int candidate0, unsigned int candidate1, const std::vector< Dune::FieldVector< T, dimworld > > &grid1Coords, const std::vector< Dune::GeometryType > &grid1_element_types, std::bitset<(1<< grid1Dim)> &neighborIntersects1, const std::vector< Dune::FieldVector< T, dimworld > > &grid2Coords, const std::vector< Dune::GeometryType > &grid2_element_types, std::bitset<(1<< grid2Dim)> &neighborIntersects2, bool insert=true)
Compute the intersection between two overlapping elements.
Definition: standardmerge.hh:259
Definition: intersectionlist.hh:131
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:127
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:126
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:64
void clear() override
Definition: standardmerge.hh:146
Definition: gridglue.hh:35
std::shared_ptr< IntersectionList > intersectionList() const final
Definition: standardmerge.hh:156
void enableBruteForce(bool bruteForce)
Definition: standardmerge.hh:167