dune-grid-glue 2.5-git
standardmerge.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
8#ifndef DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
9#define DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
10
11
12#include <iostream>
13#include <iomanip>
14#include <vector>
15#include <stack>
16#include <set>
17#include <utility>
18#include <map>
19#include <memory>
20#include <algorithm>
21
22#include <dune/common/fvector.hh>
23#include <dune/common/bitsetvector.hh>
24#include <dune/common/stdstreams.hh>
25#include <dune/common/timer.hh>
26
27#include <dune/geometry/referenceelements.hh>
28#include <dune/grid/common/grid.hh>
29
32
33namespace Dune {
34namespace GridGlue {
35
52template<class T, int grid1Dim, int grid2Dim, int dimworld>
54 : public Merger<T,grid1Dim,grid2Dim,dimworld>
55{
56
57public:
58
59 /* E X P O R T E D T Y P E S A N D C O N S T A N T S */
60
62 typedef T ctype;
63
66
69
71 typedef Dune::FieldVector<T, dimworld> WorldCoords;
72
73protected:
74
75 bool valid;
76
77 StandardMerge() : valid(false) {}
78
80 {
82 enum {intersectionDim = (grid1Dim<grid2Dim) ? grid1Dim : grid2Dim};
83
86
89 {
90 grid1Entities_.resize(1);
91 grid2Entities_.resize(1);
92 grid1Local_.resize(1);
93 grid2Local_.resize(1);
94 }
95
97 RemoteSimplicialIntersection(int grid1Entity, int grid2Entity)
98 {
99 grid1Entities_.resize(1);
100 grid2Entities_.resize(1);
101 grid1Local_.resize(1);
102 grid2Local_.resize(1);
103
104 grid1Entities_[0] = grid1Entity;
105 grid2Entities_[0] = grid2Entity;
106 }
107
108 // Local coordinates in the grid1 entity
109 std::vector<std::array<Dune::FieldVector<T,grid1Dim>, nVertices> > grid1Local_;
110
111 // Local coordinates in the grid2 entity
112 std::vector<std::array<Dune::FieldVector<T,grid2Dim>, nVertices> > grid2Local_;
113
114 //
115 std::vector<unsigned int> grid1Entities_;
116
117 std::vector<unsigned int> grid2Entities_;
118
119 };
120
125 virtual void computeIntersections(const Dune::GeometryType& grid1ElementType,
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;
134
138 bool computeIntersection(unsigned int candidate0, unsigned int candidate1,
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,
145 bool insert = true);
146
147 /* M E M B E R V A R I A B L E S */
148
150 std::vector<RemoteSimplicialIntersection> intersections_;
151
153 std::vector<std::vector<unsigned int> > grid1ElementCorners_;
154 std::vector<std::vector<unsigned int> > grid2ElementCorners_;
155
156 std::vector<std::vector<int> > elementNeighbors1_;
157 std::vector<std::vector<int> > elementNeighbors2_;
158
159public:
160
161 /* C O N C E P T I M P L E M E N T I N G I N T E R F A C E */
162
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
172 );
173
174
175 /* Q U E S T I O N I N G T H E M E R G E D G R I D */
176
179 unsigned int nSimplices() const;
180
181 void clear()
182 {
183 // Delete old internal data, from a possible previous run
184 purge(intersections_);
187
188 valid = false;
189 }
190
191 void enableFallback(bool fallback)
192 {
193 m_enableFallback = fallback;
194 }
195
196 void enableBruteForce(bool bruteForce)
197 {
198 m_enableBruteForce = bruteForce;
199 }
200
201private:
205 bool m_enableFallback = false;
206
210 bool m_enableBruteForce = false;
211
213 template<typename V>
214 static void purge(V & v)
215 {
216 v.clear();
217 V v2(v);
218 v.swap(v2);
219 }
220
221 /* M A P P I N G O N I N D E X B A S I S */
222
228 unsigned int grid1Parents(unsigned int idx) const;
229
235 unsigned int grid2Parents(unsigned int idx) const;
236
242 unsigned int grid1Parent(unsigned int idx, unsigned int parId = 0) const;
243
249 unsigned int grid2Parent(unsigned int idx, unsigned int parId = 0) const;
250
251
252 /* G E O M E T R I C A L I N F O R M A T I O N */
253
261 Grid1Coords grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
262
270 Grid2Coords grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
271
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);
283
287 int insertIntersections(unsigned int candidate1, unsigned int candidate2,std::vector<RemoteSimplicialIntersection>& intersections);
288
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);
297
301 std::pair<bool, unsigned int>
302 intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
303 RemoteSimplicialIntersection& intersection);
304
308 template <int gridDim>
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);
312
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
320 );
321
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
329 );
330};
331
332
333/* IMPLEMENTATION */
334
335template<typename T, int grid1Dim, int grid2Dim, int dimworld>
336bool StandardMerge<T,grid1Dim,grid2Dim,dimworld>::computeIntersection(unsigned int candidate0, unsigned int candidate1,
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,
343 bool insert)
344{
345 // Select vertices of the grid1 element
346 int grid1NumVertices = grid1ElementCorners_[candidate0].size();
347 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
348 for (int i=0; i<grid1NumVertices; i++)
349 grid1ElementCorners[i] = grid1Coords[grid1ElementCorners_[candidate0][i]];
350
351 // Select vertices of the grid2 element
352 int grid2NumVertices = grid2ElementCorners_[candidate1].size();
353 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
354 for (int i=0; i<grid2NumVertices; i++)
355 grid2ElementCorners[i] = grid2Coords[grid2ElementCorners_[candidate1][i]];
356
357 // ///////////////////////////////////////////////////////
358 // Compute the intersection between the two elements
359 // ///////////////////////////////////////////////////////
360
361 std::vector<RemoteSimplicialIntersection> intersections(0);
362
363 // compute the intersections
364 computeIntersections(grid1_element_types[candidate0], grid1ElementCorners,
365 neighborIntersects1, candidate0,
366 grid2_element_types[candidate1], grid2ElementCorners,
367 neighborIntersects2, candidate1,
369
370 // insert intersections if needed
371 if(insert && intersections.size() > 0)
372 insertIntersections(candidate0,candidate1,intersections);
373
374 // Have we found an intersection?
375 return (intersections.size() > 0 || neighborIntersects1.any() || neighborIntersects2.any());
376
377}
378
379template<typename T, int grid1Dim, int grid2Dim, 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)
385{
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++) {
389
390 bool intersectionFound = computeIntersection(i, candidate1,
391 grid1Coords, grid1_element_types, neighborIntersects1,
392 grid2Coords, grid2_element_types, neighborIntersects2,
393 false);
394
395 // if there is an intersection, i is our new seed candidate on the grid1 side
396 if (intersectionFound)
397 return i;
398
399 }
400
401 return -1;
402}
403
404
405template<typename T, int grid1Dim, int grid2Dim, int dimworld>
406template<int gridDim>
407void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
408computeNeighborsPerElement(const std::vector<Dune::GeometryType>& gridElementTypes,
409 const std::vector<std::vector<unsigned int> >& gridElementCorners,
410 std::vector<std::vector<int> >& elementNeighbors)
411{
412 typedef std::vector<unsigned int> FaceType;
413 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
414
416 // First: grid 1
418 FaceSetType faces;
419 elementNeighbors.resize(gridElementTypes.size());
420
421 for (size_t i=0; i<gridElementTypes.size(); i++)
422 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
423
424 for (size_t i=0; i<gridElementTypes.size(); i++) { //iterate over all elements
425 const Dune::ReferenceElement<T,gridDim>& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
426
427 for (size_t j=0; j<(size_t)refElement.size(1); j++) { // iterate over all faces of the element
428
429 FaceType face;
430 // extract element face
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)]);
433
434 // sort the face vertices to get rid of twists and other permutations
435 std::sort(face.begin(), face.end());
436
437 typename FaceSetType::iterator faceHandle = faces.find(face);
438
439 if (faceHandle == faces.end()) {
440
441 // face has not been visited before
442 faces.insert(std::make_pair(face, std::make_pair(i,j)));
443
444 } else {
445
446 // face has been visited before: store the mutual neighbor information
447 elementNeighbors[i][j] = faceHandle->second.first;
448 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
449
450 faces.erase(faceHandle);
451
452 }
453
454 }
455
456 }
457}
458
459// /////////////////////////////////////////////////////////////////////
460// Compute the intersection of all pairs of elements
461// Linear algorithm by Gander and Japhet, Proc. of DD18
462// /////////////////////////////////////////////////////////////////////
463
464template<typename T, int grid1Dim, int grid2Dim, int dimworld>
465void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
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
471 )
472{
473
474 std::cout << "StandardMerge building merged grid..." << std::endl;
475 Dune::Timer watch;
476
477 clear();
478 // clear global intersection list
479 intersections_.clear();
480 this->counter = 0;
481
482 // /////////////////////////////////////////////////////////////////////
483 // Copy element corners into a data structure with block-structure.
484 // This is not as efficient but a lot easier to use.
485 // We may think about efficiency later.
486 // /////////////////////////////////////////////////////////////////////
487
488 // first the grid1 side
489 grid1ElementCorners_.resize(grid1_element_types.size());
490
491 unsigned int grid1CornerCounter = 0;
492
493 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
494
495 // Select vertices of the grid1 element
496 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
497 grid1ElementCorners_[i].resize(numVertices);
498 for (int j=0; j<numVertices; j++)
499 grid1ElementCorners_[i][j] = grid1_elements[grid1CornerCounter++];
500
501 }
502
503 // then the grid2 side
504 grid2ElementCorners_.resize(grid2_element_types.size());
505
506 unsigned int grid2CornerCounter = 0;
507
508 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
509
510 // Select vertices of the grid2 element
511 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
512 grid2ElementCorners_[i].resize(numVertices);
513 for (int j=0; j<numVertices; j++)
514 grid2ElementCorners_[i][j] = grid2_elements[grid2CornerCounter++];
515
516 }
517
519 // Compute the face neighbors for each element
521
522 computeNeighborsPerElement<grid1Dim>(grid1_element_types, grid1ElementCorners_, elementNeighbors1_);
523 computeNeighborsPerElement<grid2Dim>(grid2_element_types, grid2ElementCorners_, elementNeighbors2_);
524
525 std::cout << "setup took " << watch.elapsed() << " seconds." << std::endl;
526
527 if (m_enableBruteForce)
528 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
529 else
530 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
531
532 valid = true;
533 std::cout << "intersection construction took " << watch.elapsed() << " seconds." << std::endl;
534}
535
536template<typename T, int grid1Dim, int grid2Dim, 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
544 )
545{
547 // Data structures for the advancing-front algorithm
549
550 std::stack<unsigned int> candidates1;
551 std::stack<unsigned int> candidates2;
552
553 std::vector<int> seeds(grid2_element_types.size(), -1);
554
555 // /////////////////////////////////////////////////////////////////////
556 // Do a brute-force search to find one pair of intersecting elements
557 // to start the advancing-front type algorithm with.
558 // /////////////////////////////////////////////////////////////////////
559
560 // Set flag if element has been handled
561 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
562
563 // Set flag if the element has been entered in the queue
564 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
565
566 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
567
568 // /////////////////////////////////////////////////////////////////////
569 // Main loop
570 // /////////////////////////////////////////////////////////////////////
571
572 std::set<unsigned int> isHandled1;
573
574 std::set<unsigned int> isCandidate1;
575
576 while (!candidates2.empty()) {
577
578 // Get the next element on the grid2 side
579 unsigned int currentCandidate2 = candidates2.top();
580 int seed = seeds[currentCandidate2];
581 assert(seed >= 0);
582
583 candidates2.pop();
584 isHandled2[currentCandidate2] = true;
585
586 // Start advancing front algorithm on the grid1 side from the 'seed' element that
587 // we stored along with the current grid2 element
588 candidates1.push(seed);
589
590 isHandled1.clear();
591 isCandidate1.clear();
592
593 while (!candidates1.empty()) {
594
595 unsigned int currentCandidate1 = candidates1.top();
596 candidates1.pop();
597 isHandled1.insert(currentCandidate1);
598
599 // Test whether there is an intersection between currentCandidate0 and currentCandidate1
600 std::bitset<(1<<grid1Dim)> neighborIntersects1;
601 std::bitset<(1<<grid2Dim)> neighborIntersects2;
602 bool intersectionFound = computeIntersection(currentCandidate1, currentCandidate2,
603 grid1Coords,grid1_element_types, neighborIntersects1,
604 grid2Coords,grid2_element_types, neighborIntersects2);
605
606 for (size_t i=0; i<neighborIntersects2.size(); i++)
607 if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
608 seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
609
610 // add neighbors of candidate0 to the list of elements to be checked
611 if (intersectionFound) {
612
613 for (size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
614
615 int neighbor = elementNeighbors1_[currentCandidate1][i];
616
617 if (neighbor == -1) // do nothing at the grid boundary
618 continue;
619
620 if (isHandled1.find(neighbor) == isHandled1.end()
621 && isCandidate1.find(neighbor) == isCandidate1.end()) {
622 candidates1.push(neighbor);
623 isCandidate1.insert(neighbor);
624 }
625
626 }
627
628 }
629
630 }
631
632 // We have now found all intersections of elements in the grid1 side with currentCandidate2
633 // Now we add all neighbors of currentCandidate2 that have not been treated yet as new
634 // candidates.
635
636 // Do we have an unhandled neighbor with a seed?
637 bool seedFound = !candidates2.empty();
638 for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
639
640 int neighbor = elementNeighbors2_[currentCandidate2][i];
641
642 if (neighbor == -1) // do nothing at the grid boundary
643 continue;
644
645 // Add all unhandled intersecting neighbors to the queue
646 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
647
648 isCandidate2[neighbor][0] = true;
649 candidates2.push(neighbor);
650 seedFound = true;
651 }
652 }
653
654 if (seedFound || !m_enableFallback)
655 continue;
656
657 // There is no neighbor with a seed, so we need to be a bit more aggressive...
658 // get all neighbors of currentCandidate2, but not currentCandidate2 itself
659 for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
660
661 int neighbor = elementNeighbors2_[currentCandidate2][i];
662
663 if (neighbor == -1) // do nothing at the grid boundary
664 continue;
665
666 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
667
668 // Get a seed element for the new grid2 element
669 // Look for an element on the grid1 side that intersects the new grid2 element.
670 int seed = -1;
671
672 // Look among the ones that have been tested during the last iteration.
673 for (typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
674 seedIt != isHandled1.end(); ++seedIt) {
675
676 std::bitset<(1<<grid1Dim)> neighborIntersects1;
677 std::bitset<(1<<grid2Dim)> neighborIntersects2;
678 bool intersectionFound = computeIntersection(*seedIt, neighbor,
679 grid1Coords, grid1_element_types, neighborIntersects1,
680 grid2Coords, grid2_element_types, neighborIntersects2,
681 false);
682
683 // if the intersection is nonempty, *seedIt is our new seed candidate on the grid1 side
684 if (intersectionFound) {
685 seed = *seedIt;
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;
688 break;
689 }
690
691 }
692
693 if (seed < 0) {
694 // The fast method didn't find a grid1 element that intersects with
695 // the new grid2 candidate. We have to do a brute-force search.
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;
700
701 }
702
703 // We have tried all we could: the candidate is 'handled' now
704 isCandidate2[neighbor] = true;
705
706 // still no seed? Then the new grid2 candidate isn't overlapped by anything
707 if (seed < 0)
708 continue;
709
710 // we have a seed now
711 candidates2.push(neighbor);
712 seeds[neighbor] = seed;
713 seedFound = true;
714
715 }
716
717 }
718
719 /* Do a brute-force search if there is still no seed:
720 * There might still be a disconnected region out there.
721 */
722 if (!seedFound && candidates2.empty()) {
723 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
724 }
725 }
726}
727
728template<typename T, int grid1Dim, int grid2Dim, int dimworld>
729void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildBruteForce(
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
736 )
737{
738 std::bitset<(1<<grid1Dim)> neighborIntersects1;
739 std::bitset<(1<<grid2Dim)> neighborIntersects2;
740
741 for (unsigned i = 0; i < grid1_element_types.size(); ++i) {
742 for (unsigned j = 0; j < grid2_element_types.size(); ++j) {
743 (void) computeIntersection(i, j,
744 grid1Coords, grid1_element_types, neighborIntersects1,
745 grid2Coords, grid2_element_types, neighborIntersects2);
746 }
747 }
748}
749
750template<typename T, int grid1Dim, int grid2Dim, int dimworld>
752{
753 assert(valid);
754 return intersections_.size();
755}
756
757template<typename T, int grid1Dim, int grid2Dim, int dimworld>
758inline unsigned int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid1Parents(unsigned int idx) const
759{
760 assert(valid);
761 return (intersections_[idx].grid1Entities_).size();
762}
763
764template<typename T, int grid1Dim, int grid2Dim, int dimworld>
765inline unsigned int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid2Parents(unsigned int idx) const
766{
767 assert(valid);
768 return (intersections_[idx].grid2Entities_).size();
769}
770
771
772template<typename T, int grid1Dim, int grid2Dim, int dimworld>
773inline unsigned int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid1Parent(unsigned int idx, unsigned int parId) const
774{
775 assert(valid);
776 return intersections_[idx].grid1Entities_[parId];
777}
778
779
780template<typename T, int grid1Dim, int grid2Dim, int dimworld>
781inline unsigned int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid2Parent(unsigned int idx, unsigned int parId) const
782{
783 assert(valid);
784 return intersections_[idx].grid2Entities_[parId];
785}
786
787
788template<typename T, int grid1Dim, int grid2Dim, int dimworld>
789typename StandardMerge<T,grid1Dim,grid2Dim,dimworld>::Grid1Coords StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
790{
791 assert(valid);
792 return intersections_[idx].grid1Local_[parId][corner];
793}
794
795
796template<typename T, int grid1Dim, int grid2Dim, int dimworld>
797typename StandardMerge<T,grid1Dim,grid2Dim,dimworld>::Grid2Coords StandardMerge<T,grid1Dim,grid2Dim,dimworld>::grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
798{
799 assert(valid);
800 return intersections_[idx].grid2Local_[parId][corner];
801}
802
803template<typename T, int grid1Dim, int grid2Dim, int dimworld>
804void 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)
805{
806 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
807
808 if (seeds[j] > 0 || isHandled2[j][0])
809 continue;
810
811 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
812
813 if (seed >= 0) {
814 candidates2.push(j); // the candidate and a seed for the candidate
815 seeds[j] = seed;
816 break;
817 } else // If the brute force search did not find any intersection we can skip this element
818 isHandled2[j] = true;
819 }
820}
821
822template<typename T, int grid1Dim, int grid2Dim, int dimworld>
823int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::insertIntersections(unsigned int candidate1, unsigned int candidate2,
824 std::vector<RemoteSimplicialIntersection>& intersections)
825{
826 typedef typename std::vector<RemoteSimplicialIntersection>::size_type size_t;
827 int count = 0;
828
829 for (size_t i = 0; i < intersections.size(); ++i) {
830 // get the intersection index of the current intersection from intersections in this->intersections
831 bool found;
832 unsigned int index;
833 std::tie(found, index) = intersectionIndex(candidate1,candidate2,intersections[i]);
834
835 if (found && index >= this->intersections_.size()) { //the intersection is not yet contained in this->intersections
836 this->intersections_.push_back(intersections[i]); // insert
837
838 ++count;
839 } else if (found) {
840 // insert each grid1 element and local representation of intersections[i] with parent candidate1
841 for (size_t j = 0; j < intersections[i].grid1Entities_.size(); ++j) {
842 this->intersections_[index].grid1Entities_.push_back(candidate1);
843 this->intersections_[index].grid1Local_.push_back(intersections[i].grid1Local_[j]);
844 }
845
846 // insert each grid2 element and local representation of intersections[i] with parent candidate2
847 for (size_t j = 0; j < intersections[i].grid2Entities_.size(); ++j) {
848 this->intersections_[index].grid2Entities_.push_back(candidate2);
849 this->intersections_[index].grid2Local_.push_back(intersections[i].grid2Local_[j]);
850 }
851
852 ++count;
853 } else {
854 Dune::dwarn << "Computed the same intersection twice!" << std::endl;
855 }
856 }
857 return count;
858}
859
860template<typename T, int grid1Dim, int grid2Dim, int dimworld>
861std::pair<bool, unsigned int>
862StandardMerge<T,grid1Dim,grid2Dim,dimworld>::intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
863 RemoteSimplicialIntersection& intersection) {
864
865
866 // return index in intersections_ if at least one local representation of a Remote Simplicial Intersection (RSI)
867 // of intersections_ is equal to the local representation of one element in intersections
868
869 std::size_t n_intersections = this->intersections_.size();
870 if (grid1Dim == grid2Dim)
871 return {true, n_intersections};
872
873 T eps = 1e-10;
874
875 for (std::size_t i = 0; i < n_intersections; ++i) {
876
877 // compare the local representation of the subelements of the RSI
878 for (std::size_t ei = 0; ei < this->intersections_[i].grid1Entities_.size(); ++ei) // merger subelement
879 {
880 if (this->intersections_[i].grid1Entities_[ei] == grid1Index)
881 {
882 for (std::size_t er = 0; er < intersection.grid1Entities_.size(); ++er) // list subelement
883 {
884 bool found_all = true;
885 // compare the local coordinate representations
886 for (std::size_t ci = 0; ci < this->intersections_[i].grid1Local_[ei].size(); ++ci)
887 {
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)
891 {
892 Dune::FieldVector<T,grid1Dim> nr = intersection.grid1Local_[er][cr];
893
894 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
895 if (found_ni)
896 break;
897 }
898 found_all = found_all && found_ni;
899
900 if (!found_ni)
901 break;
902 }
903
904 if (found_all && (this->intersections_[i].grid2Entities_[ei] != grid2Index))
905 return {true, i};
906 else if (found_all)
907 return {false, 0};
908 }
909 }
910 }
911
912 // compare the local representation of the subelements of the RSI
913 for (std::size_t ei = 0; ei < this->intersections_[i].grid2Entities_.size(); ++ei) // merger subelement
914 {
915 if (this->intersections_[i].grid2Entities_[ei] == grid2Index)
916 {
917 for (std::size_t er = 0; er < intersection.grid2Entities_.size(); ++er) // list subelement
918 {
919 bool found_all = true;
920 // compare the local coordinate representations
921 for (std::size_t ci = 0; ci < this->intersections_[i].grid2Local_[ei].size(); ++ci)
922 {
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)
926 {
927 Dune::FieldVector<T,grid2Dim> nr = intersection.grid2Local_[er][cr];
928 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
929
930 if (found_ni)
931 break;
932 }
933 found_all = found_all && found_ni;
934
935 if (!found_ni)
936 break;
937 }
938
939 if (found_all && (this->intersections_[i].grid1Entities_[ei] != grid1Index))
940 return {true, i};
941 else if (found_all)
942 return {false, 0};
943 }
944 }
945 }
946 }
947
948 return {true, n_intersections};
949}
950
951#define DECL extern
952#define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \
953 DECL template \
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 \
960 )
961
962STANDARD_MERGE_INSTANTIATE(double,1,1,1);
963STANDARD_MERGE_INSTANTIATE(double,2,2,2);
964STANDARD_MERGE_INSTANTIATE(double,3,3,3);
965#undef STANDARD_MERGE_INSTANTIATE
966#undef DECL
967
968} /* namespace GridGlue */
969} /* namespace Dune */
970
971#endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
#define STANDARD_MERGE_INSTANTIATE(T, A, B, C)
Definition standardmerge.cc:11
Definition gridglue.hh:30
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Coordinate corner(unsigned c)
Definition projection_impl.hh:22
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition merger.hh:91
Dune::FieldVector< ctype, grid1Dim > Grid1Coords
the local coordinate type for the grid1 coordinates
Definition merger.hh:100
Dune::FieldVector< ctype, grid2Dim > Grid2Coords
the local coordinate type for the grid2 coordinates
Definition merger.hh:103
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition standardmerge.hh:55
void enableFallback(bool fallback)
Definition standardmerge.hh:191
std::vector< std::vector< int > > elementNeighbors2_
Definition standardmerge.hh:157
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid1Coords Grid1Coords
Type used for local coordinates on the grid1 side.
Definition standardmerge.hh:65
std::vector< std::vector< int > > elementNeighbors1_
Definition standardmerge.hh:156
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid2Coords Grid2Coords
Type used for local coordinates on the grid2 side.
Definition standardmerge.hh:68
T ctype
the numeric type used in this interface
Definition standardmerge.hh:62
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
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)
Definition standardmerge.hh:465
void enableBruteForce(bool bruteForce)
Definition standardmerge.hh:196
void clear()
Definition standardmerge.hh:181
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.
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition standardmerge.hh:153
std::vector< RemoteSimplicialIntersection > intersections_
The computed intersections.
Definition standardmerge.hh:150
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition standardmerge.hh:154
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 valid
Definition standardmerge.hh:75
StandardMerge()
Definition standardmerge.hh:77
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition standardmerge.hh:71
std::vector< unsigned int > grid1Entities_
Definition standardmerge.hh:115
RemoteSimplicialIntersection()
Default constructor.
Definition standardmerge.hh:88
std::vector< std::array< Dune::FieldVector< T, grid1Dim >, nVertices > > grid1Local_
Definition standardmerge.hh:109
std::vector< unsigned int > grid2Entities_
Definition standardmerge.hh:117
RemoteSimplicialIntersection(int grid1Entity, int grid2Entity)
Constructor for two given entity indices.
Definition standardmerge.hh:97
std::vector< std::array< Dune::FieldVector< T, grid2Dim >, nVertices > > grid2Local_
Definition standardmerge.hh:112