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> 52 template<
class T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
54 :
public Merger<T,grid1Dim,grid2Dim,dimworld>
126 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
127 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
128 unsigned int grid1Index,
129 const Dune::GeometryType& grid2ElementType,
130 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
131 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
132 unsigned int grid2Index,
133 std::vector<RemoteSimplicialIntersection>&
intersections) = 0;
139 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
140 const std::vector<Dune::GeometryType>& grid1_element_types,
141 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
142 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
143 const std::vector<Dune::GeometryType>& grid2_element_types,
144 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
166 virtual void build(
const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
167 const std::vector<unsigned int>& grid1_elements,
168 const std::vector<Dune::GeometryType>& grid1_element_types,
169 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
170 const std::vector<unsigned int>& grid2_elements,
171 const std::vector<Dune::GeometryType>& grid2_element_types
184 purge(intersections_);
185 purge(grid1ElementCorners_);
186 purge(grid2ElementCorners_);
193 m_enableFallback = fallback;
198 m_enableBruteForce = bruteForce;
205 bool m_enableFallback =
false;
210 bool m_enableBruteForce =
false;
214 static void purge(V & v)
228 unsigned int grid1Parents(
unsigned int idx)
const;
235 unsigned int grid2Parents(
unsigned int idx)
const;
242 unsigned int grid1Parent(
unsigned int idx,
unsigned int parId = 0)
const;
249 unsigned int grid2Parent(
unsigned int idx,
unsigned int parId = 0)
const;
261 Grid1Coords grid1ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
270 Grid2Coords grid2ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
276 void generateSeed(std::vector<int>& seeds,
277 Dune::BitSetVector<1>& isHandled2,
278 std::stack<unsigned>& candidates2,
279 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
280 const std::vector<Dune::GeometryType>& grid1_element_types,
281 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
282 const std::vector<Dune::GeometryType>& grid2_element_types);
287 int insertIntersections(
unsigned int candidate1,
unsigned int candidate2,std::vector<RemoteSimplicialIntersection>&
intersections);
292 int bruteForceSearch(
int candidate1,
293 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
294 const std::vector<Dune::GeometryType>& grid1_element_types,
295 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
296 const std::vector<Dune::GeometryType>& grid2_element_types);
301 std::pair<bool, unsigned int>
302 intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
308 template <
int gr
idDim>
309 void computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
310 const std::vector<std::vector<unsigned int> >& gridElementCorners,
311 std::vector<std::vector<int> >& elementNeighbors);
313 void buildAdvancingFront(
314 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
315 const std::vector<unsigned int>& grid1_elements,
316 const std::vector<Dune::GeometryType>& grid1_element_types,
317 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
318 const std::vector<unsigned int>& grid2_elements,
319 const std::vector<Dune::GeometryType>& grid2_element_types
322 void buildBruteForce(
323 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
324 const std::vector<unsigned int>& grid1_elements,
325 const std::vector<Dune::GeometryType>& grid1_element_types,
326 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
327 const std::vector<unsigned int>& grid2_elements,
328 const std::vector<Dune::GeometryType>& grid2_element_types
335 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
337 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
338 const std::vector<Dune::GeometryType>& grid1_element_types,
339 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
340 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
341 const std::vector<Dune::GeometryType>& grid2_element_types,
342 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
347 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
348 for (
int i=0; i<grid1NumVertices; i++)
353 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
354 for (
int i=0; i<grid2NumVertices; i++)
365 neighborIntersects1, candidate0,
366 grid2_element_types[candidate1], grid2ElementCorners,
367 neighborIntersects2, candidate1,
371 if(insert && intersections.size() > 0)
372 insertIntersections(candidate0,candidate1,intersections);
375 return (intersections.size() > 0 || neighborIntersects1.any() || neighborIntersects2.any());
379 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
381 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
382 const std::vector<Dune::GeometryType>& grid1_element_types,
383 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
384 const std::vector<Dune::GeometryType>& grid2_element_types)
386 std::bitset<(1<<grid1Dim)> neighborIntersects1;
387 std::bitset<(1<<grid2Dim)> neighborIntersects2;
388 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
391 grid1Coords, grid1_element_types, neighborIntersects1,
392 grid2Coords, grid2_element_types, neighborIntersects2,
396 if (intersectionFound)
405 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
406 template<
int gr
idDim>
409 const std::vector<std::vector<unsigned int> >& gridElementCorners,
410 std::vector<std::vector<int> >& elementNeighbors)
412 typedef std::vector<unsigned int> FaceType;
413 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
419 elementNeighbors.resize(gridElementTypes.size());
421 for (
size_t i=0; i<gridElementTypes.size(); i++)
422 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
424 for (
size_t i=0; i<gridElementTypes.size(); i++) {
425 const Dune::ReferenceElement<T,gridDim>& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
427 for (
size_t j=0; j<(size_t)refElement.size(1); j++) {
431 for (
size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
432 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
435 std::sort(face.begin(), face.end());
437 typename FaceSetType::iterator faceHandle = faces.find(face);
439 if (faceHandle == faces.end()) {
442 faces.insert(std::make_pair(face, std::make_pair(i,j)));
447 elementNeighbors[i][j] = faceHandle->second.first;
448 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
450 faces.erase(faceHandle);
464 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
466 const std::vector<unsigned int>& grid1_elements,
467 const std::vector<Dune::GeometryType>& grid1_element_types,
468 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
469 const std::vector<unsigned int>& grid2_elements,
470 const std::vector<Dune::GeometryType>& grid2_element_types
474 std::cout <<
"StandardMerge building merged grid..." << std::endl;
491 unsigned int grid1CornerCounter = 0;
493 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
496 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
498 for (
int j=0; j<numVertices; j++)
506 unsigned int grid2CornerCounter = 0;
508 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
511 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
513 for (
int j=0; j<numVertices; j++)
525 std::cout <<
"setup took " << watch.elapsed() <<
" seconds." << std::endl;
527 if (m_enableBruteForce)
528 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
530 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
533 std::cout <<
"intersection construction took " << watch.elapsed() <<
" seconds." << std::endl;
536 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
538 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
539 const std::vector<unsigned int>& grid1_elements,
540 const std::vector<Dune::GeometryType>& grid1_element_types,
541 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
542 const std::vector<unsigned int>& grid2_elements,
543 const std::vector<Dune::GeometryType>& grid2_element_types
550 std::stack<unsigned int> candidates1;
551 std::stack<unsigned int> candidates2;
553 std::vector<int> seeds(grid2_element_types.size(), -1);
561 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
564 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
566 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
572 std::set<unsigned int> isHandled1;
574 std::set<unsigned int> isCandidate1;
576 while (!candidates2.empty()) {
579 unsigned int currentCandidate2 = candidates2.top();
580 int seed = seeds[currentCandidate2];
584 isHandled2[currentCandidate2] =
true;
588 candidates1.push(seed);
591 isCandidate1.clear();
593 while (!candidates1.empty()) {
595 unsigned int currentCandidate1 = candidates1.top();
597 isHandled1.insert(currentCandidate1);
600 std::bitset<(1<<grid1Dim)> neighborIntersects1;
601 std::bitset<(1<<grid2Dim)> neighborIntersects2;
603 grid1Coords,grid1_element_types, neighborIntersects1,
604 grid2Coords,grid2_element_types, neighborIntersects2);
606 for (
size_t i=0; i<neighborIntersects2.size(); i++)
611 if (intersectionFound) {
620 if (isHandled1.find(neighbor) == isHandled1.end()
621 && isCandidate1.find(neighbor) == isCandidate1.end()) {
622 candidates1.push(neighbor);
623 isCandidate1.insert(neighbor);
637 bool seedFound = !candidates2.empty();
646 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
648 isCandidate2[neighbor][0] =
true;
649 candidates2.push(neighbor);
654 if (seedFound || !m_enableFallback)
666 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
673 for (
typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
674 seedIt != isHandled1.end(); ++seedIt) {
676 std::bitset<(1<<grid1Dim)> neighborIntersects1;
677 std::bitset<(1<<grid2Dim)> neighborIntersects2;
679 grid1Coords, grid1_element_types, neighborIntersects1,
680 grid2Coords, grid2_element_types, neighborIntersects2,
684 if (intersectionFound) {
686 Dune::dwarn <<
"Algorithm entered first fallback method and found a new seed in the build algorithm." <<
687 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
696 seed = bruteForceSearch(neighbor,
697 grid1Coords,grid1_element_types,
698 grid2Coords,grid2_element_types);
699 Dune::dwarn <<
"Algorithm entered second fallback method. This probably should not happen." << std::endl;
704 isCandidate2[neighbor] =
true;
711 candidates2.push(neighbor);
712 seeds[neighbor] = seed;
722 if (!seedFound && candidates2.empty()) {
723 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
728 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
730 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
731 const std::vector<unsigned int>& grid1_elements,
732 const std::vector<Dune::GeometryType>& grid1_element_types,
733 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
734 const std::vector<unsigned int>& grid2_elements,
735 const std::vector<Dune::GeometryType>& grid2_element_types
738 std::bitset<(1<<grid1Dim)> neighborIntersects1;
739 std::bitset<(1<<grid2Dim)> neighborIntersects2;
741 for (
unsigned i = 0; i < grid1_element_types.size(); ++i) {
742 for (
unsigned j = 0; j < grid2_element_types.size(); ++j) {
744 grid1Coords, grid1_element_types, neighborIntersects1,
745 grid2Coords, grid2_element_types, neighborIntersects2);
750 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
757 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
764 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
772 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
780 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
788 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
796 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
803 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
804 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)
806 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
808 if (seeds[j] > 0 || isHandled2[j][0])
811 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
818 isHandled2[j] =
true;
822 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
826 typedef typename std::vector<RemoteSimplicialIntersection>::size_type size_t;
829 for (
size_t i = 0; i < intersections.size(); ++i) {
833 std::tie(found, index) = intersectionIndex(candidate1,candidate2,intersections[i]);
841 for (
size_t j = 0; j < intersections[i].grid1Entities_.size(); ++j) {
847 for (
size_t j = 0; j < intersections[i].grid2Entities_.size(); ++j) {
854 Dune::dwarn <<
"Computed the same intersection twice!" << std::endl;
860 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
861 std::pair<bool, unsigned int>
870 if (grid1Dim == grid2Dim)
871 return {
true, n_intersections};
875 for (std::size_t i = 0; i < n_intersections; ++i) {
878 for (std::size_t ei = 0; ei < this->
intersections_[i].grid1Entities_.size(); ++ei)
882 for (std::size_t er = 0; er < intersection.grid1Entities_.size(); ++er)
884 bool found_all =
true;
886 for (std::size_t ci = 0; ci < this->
intersections_[i].grid1Local_[ei].size(); ++ci)
888 Dune::FieldVector<T,grid1Dim> ni = this->
intersections_[i].grid1Local_[ei][ci];
889 bool found_ni =
false;
890 for (std::size_t cr = 0; cr < intersection.grid1Local_[er].size(); ++cr)
892 Dune::FieldVector<T,grid1Dim> nr = intersection.grid1Local_[er][cr];
894 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
898 found_all = found_all && found_ni;
913 for (std::size_t ei = 0; ei < this->
intersections_[i].grid2Entities_.size(); ++ei)
917 for (std::size_t er = 0; er < intersection.grid2Entities_.size(); ++er)
919 bool found_all =
true;
921 for (std::size_t ci = 0; ci < this->
intersections_[i].grid2Local_[ei].size(); ++ci)
923 Dune::FieldVector<T,grid2Dim> ni = this->
intersections_[i].grid2Local_[ei][ci];
924 bool found_ni =
false;
925 for (std::size_t cr = 0; cr < intersection.grid2Local_[er].size(); ++cr)
927 Dune::FieldVector<T,grid2Dim> nr = intersection.grid2Local_[er][cr];
928 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
933 found_all = found_all && found_ni;
948 return {
true, n_intersections};
952 #define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \ 954 void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \ 955 const std::vector<unsigned int>& grid1_elements, \ 956 const std::vector<Dune::GeometryType>& grid1_element_types, \ 957 const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \ 958 const std::vector<unsigned int>& grid2_elements, \ 959 const std::vector<Dune::GeometryType>& grid2_element_types \ 965 #undef STANDARD_MERGE_INSTANTIATE 971 #endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH unsigned int nSimplices() const
get the number of simplices in the merged grid The indices are then in 0..nSimplices()-1 ...
Definition: standardmerge.hh:751
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:336
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:53
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< RemoteSimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
unsigned int counter
Counts the number of times the computeIntersection method has been called.
Definition: merger.hh:195
RemoteSimplicialIntersection()
Default constructor.
Definition: standardmerge.hh:88
StandardMerge()
Definition: standardmerge.hh:77
STANDARD_MERGE_INSTANTIATE(double, 1, 1, 1)
void enableFallback(bool fallback)
Definition: standardmerge.hh:191
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:153
std::vector< unsigned int > grid2Entities_
Definition: standardmerge.hh:117
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:154
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:157
std::vector< std::array< Dune::FieldVector< T, grid2Dim >, nVertices > > grid2Local_
Definition: standardmerge.hh:112
Definition: standardmerge.hh:85
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:71
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
bool valid
Definition: standardmerge.hh:75
Definition: standardmerge.hh:79
void clear()
Definition: standardmerge.hh:181
RemoteSimplicialIntersection(int grid1Entity, int grid2Entity)
Constructor for two given entity indices.
Definition: standardmerge.hh:97
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:62
std::vector< unsigned int > grid1Entities_
Definition: standardmerge.hh:115
Definition: standardmerge.hh:82
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:16
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:156
virtual 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)
builds the merged grid
Definition: standardmerge.hh:465
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid1Coords Grid1Coords
Type used for local coordinates on the grid1 side.
Definition: standardmerge.hh:65
Coordinate corner(unsigned c)
Definition: projection_impl.hh:22
std::vector< std::array< Dune::FieldVector< T, grid1Dim >, nVertices > > grid1Local_
Definition: standardmerge.hh:109
Definition: gridglue.hh:30
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid2Coords Grid2Coords
Type used for local coordinates on the grid2 side.
Definition: standardmerge.hh:68
void enableBruteForce(bool bruteForce)
Definition: standardmerge.hh:196
std::vector< RemoteSimplicialIntersection > intersections_
The computed intersections.
Definition: standardmerge.hh:150