dune-grid-glue 2.8-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
33
34namespace Dune {
35namespace GridGlue {
36
53template<class T, int grid1Dim, int grid2Dim, int dimworld>
55 : public Merger<T,grid1Dim,grid2Dim,dimworld>
56{
58
59public:
60
61 /* 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 */
62
64 typedef T ctype;
65
68
71
74
76
77protected:
78
83
84 bool valid = false;
85
89 {}
90
91 virtual ~StandardMerge() = default;
92
97 virtual void computeIntersections(const Dune::GeometryType& grid1ElementType,
98 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
99 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
100 unsigned int grid1Index,
101 const Dune::GeometryType& grid2ElementType,
102 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
103 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
104 unsigned int grid2Index,
105 std::vector<SimplicialIntersection>& intersections) = 0;
106
110 bool computeIntersection(unsigned int candidate0, unsigned int candidate1,
111 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
112 const std::vector<Dune::GeometryType>& grid1_element_types,
113 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
114 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
115 const std::vector<Dune::GeometryType>& grid2_element_types,
116 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
117 bool insert = true);
118
119 /* M E M B E R V A R I A B L E S */
120
121 std::shared_ptr<IntersectionListProvider> intersectionListProvider_;
122 std::shared_ptr<IntersectionList> intersectionList_;
123
125 std::vector<std::vector<unsigned int> > grid1ElementCorners_;
126 std::vector<std::vector<unsigned int> > grid2ElementCorners_;
127
128 std::vector<std::vector<int> > elementNeighbors1_;
129 std::vector<std::vector<int> > elementNeighbors2_;
130
131public:
132
133 /* 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 */
134
138 void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
139 const std::vector<unsigned int>& grid1_elements,
140 const std::vector<Dune::GeometryType>& grid1_element_types,
141 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
142 const std::vector<unsigned int>& grid2_elements,
143 const std::vector<Dune::GeometryType>& grid2_element_types) override;
144
145
146 /* P R O B I N G T H E M E R G E D G R I D */
147
148 void clear() override
149 {
150 // Delete old internal data, from a possible previous run
154
155 valid = false;
156 }
157
158 std::shared_ptr<IntersectionList> intersectionList() const final
159 {
160 assert(valid);
161 return intersectionList_;
162 }
163
164 void enableFallback(bool fallback)
165 {
166 m_enableFallback = fallback;
167 }
168
169 void enableBruteForce(bool bruteForce)
170 {
171 m_enableBruteForce = bruteForce;
172 }
173
174private:
178 bool m_enableFallback = false;
179
183 bool m_enableBruteForce = false;
184
185 auto& intersections()
186 { return intersectionListProvider_->intersections(); }
187
189 template<typename V>
190 static void purge(V & v)
191 {
192 v.clear();
193 V v2(v);
194 v.swap(v2);
195 }
196
201 void generateSeed(std::vector<int>& seeds,
202 Dune::BitSetVector<1>& isHandled2,
203 std::stack<unsigned>& candidates2,
204 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
205 const std::vector<Dune::GeometryType>& grid1_element_types,
206 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
207 const std::vector<Dune::GeometryType>& grid2_element_types);
208
212 int insertIntersections(unsigned int candidate1, unsigned int candidate2,std::vector<SimplicialIntersection>& intersections);
213
217 int bruteForceSearch(int candidate1,
218 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
219 const std::vector<Dune::GeometryType>& grid1_element_types,
220 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
221 const std::vector<Dune::GeometryType>& grid2_element_types);
222
226 std::pair<bool, unsigned int>
227 intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
228 SimplicialIntersection& intersection);
229
233 template <int gridDim>
234 void computeNeighborsPerElement(const std::vector<Dune::GeometryType>& gridElementTypes,
235 const std::vector<std::vector<unsigned int> >& gridElementCorners,
236 std::vector<std::vector<int> >& elementNeighbors);
237
238 void buildAdvancingFront(
239 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
240 const std::vector<unsigned int>& grid1_elements,
241 const std::vector<Dune::GeometryType>& grid1_element_types,
242 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
243 const std::vector<unsigned int>& grid2_elements,
244 const std::vector<Dune::GeometryType>& grid2_element_types
245 );
246
247 void buildBruteForce(
248 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
249 const std::vector<unsigned int>& grid1_elements,
250 const std::vector<Dune::GeometryType>& grid1_element_types,
251 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
252 const std::vector<unsigned int>& grid2_elements,
253 const std::vector<Dune::GeometryType>& grid2_element_types
254 );
255};
256
257
258/* IMPLEMENTATION */
259
260template<typename T, int grid1Dim, int grid2Dim, int dimworld>
261bool StandardMerge<T,grid1Dim,grid2Dim,dimworld>::computeIntersection(unsigned int candidate0, unsigned int candidate1,
262 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
263 const std::vector<Dune::GeometryType>& grid1_element_types,
264 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
265 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
266 const std::vector<Dune::GeometryType>& grid2_element_types,
267 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
268 bool insert)
269{
270 // Select vertices of the grid1 element
271 int grid1NumVertices = grid1ElementCorners_[candidate0].size();
272 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
273 for (int i=0; i<grid1NumVertices; i++)
274 grid1ElementCorners[i] = grid1Coords[grid1ElementCorners_[candidate0][i]];
275
276 // Select vertices of the grid2 element
277 int grid2NumVertices = grid2ElementCorners_[candidate1].size();
278 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
279 for (int i=0; i<grid2NumVertices; i++)
280 grid2ElementCorners[i] = grid2Coords[grid2ElementCorners_[candidate1][i]];
281
282 // ///////////////////////////////////////////////////////
283 // Compute the intersection between the two elements
284 // ///////////////////////////////////////////////////////
285
286 std::vector<SimplicialIntersection> intersections(0);
287
288 // compute the intersections
289 computeIntersections(grid1_element_types[candidate0], grid1ElementCorners,
290 neighborIntersects1, candidate0,
291 grid2_element_types[candidate1], grid2ElementCorners,
292 neighborIntersects2, candidate1,
294
295 // insert intersections if needed
296 if(insert && !intersections.empty())
297 insertIntersections(candidate0,candidate1,intersections);
298
299 // Have we found an intersection?
300 return !intersections.empty() || neighborIntersects1.any() || neighborIntersects2.any();
301
302}
303
304template<typename T, int grid1Dim, int grid2Dim, int dimworld>
306 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
307 const std::vector<Dune::GeometryType>& grid1_element_types,
308 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
309 const std::vector<Dune::GeometryType>& grid2_element_types)
310{
311 std::bitset<(1<<grid1Dim)> neighborIntersects1;
312 std::bitset<(1<<grid2Dim)> neighborIntersects2;
313 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
314
315 bool intersectionFound = computeIntersection(i, candidate1,
316 grid1Coords, grid1_element_types, neighborIntersects1,
317 grid2Coords, grid2_element_types, neighborIntersects2,
318 false);
319
320 // if there is an intersection, i is our new seed candidate on the grid1 side
321 if (intersectionFound)
322 return i;
323
324 }
325
326 return -1;
327}
328
329
330template<typename T, int grid1Dim, int grid2Dim, int dimworld>
331template<int gridDim>
332void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
333computeNeighborsPerElement(const std::vector<Dune::GeometryType>& gridElementTypes,
334 const std::vector<std::vector<unsigned int> >& gridElementCorners,
335 std::vector<std::vector<int> >& elementNeighbors)
336{
337 typedef std::vector<unsigned int> FaceType;
338 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
339
341 // First: grid 1
343 FaceSetType faces;
344 elementNeighbors.resize(gridElementTypes.size());
345
346 for (size_t i=0; i<gridElementTypes.size(); i++)
347 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
348
349 for (size_t i=0; i<gridElementTypes.size(); i++) { //iterate over all elements
350 const auto& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
351
352 for (size_t j=0; j<(size_t)refElement.size(1); j++) { // iterate over all faces of the element
353
354 FaceType face;
355 // extract element face
356 for (size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
357 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
358
359 // sort the face vertices to get rid of twists and other permutations
360 std::sort(face.begin(), face.end());
361
362 typename FaceSetType::iterator faceHandle = faces.find(face);
363
364 if (faceHandle == faces.end()) {
365
366 // face has not been visited before
367 faces.insert(std::make_pair(face, std::make_pair(i,j)));
368
369 } else {
370
371 // face has been visited before: store the mutual neighbor information
372 elementNeighbors[i][j] = faceHandle->second.first;
373 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
374
375 faces.erase(faceHandle);
376
377 }
378
379 }
380
381 }
382}
383
384// /////////////////////////////////////////////////////////////////////
385// Compute the intersection of all pairs of elements
386// Linear algorithm by Gander and Japhet, Proc. of DD18
387// /////////////////////////////////////////////////////////////////////
388
389template<typename T, int grid1Dim, int grid2Dim, int dimworld>
390void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
391 const std::vector<unsigned int>& grid1_elements,
392 const std::vector<Dune::GeometryType>& grid1_element_types,
393 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
394 const std::vector<unsigned int>& grid2_elements,
395 const std::vector<Dune::GeometryType>& grid2_element_types
396 )
397{
398
399 std::cout << "StandardMerge building merged grid..." << std::endl;
400 Dune::Timer watch;
401
402 clear();
403 // clear global intersection list
404 intersectionListProvider_->clear();
405 this->counter = 0;
406
407 // /////////////////////////////////////////////////////////////////////
408 // Copy element corners into a data structure with block-structure.
409 // This is not as efficient but a lot easier to use.
410 // We may think about efficiency later.
411 // /////////////////////////////////////////////////////////////////////
412
413 // first the grid1 side
414 grid1ElementCorners_.resize(grid1_element_types.size());
415
416 unsigned int grid1CornerCounter = 0;
417
418 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
419
420 // Select vertices of the grid1 element
421 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
422 grid1ElementCorners_[i].resize(numVertices);
423 for (int j=0; j<numVertices; j++)
424 grid1ElementCorners_[i][j] = grid1_elements[grid1CornerCounter++];
425
426 }
427
428 // then the grid2 side
429 grid2ElementCorners_.resize(grid2_element_types.size());
430
431 unsigned int grid2CornerCounter = 0;
432
433 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
434
435 // Select vertices of the grid2 element
436 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
437 grid2ElementCorners_[i].resize(numVertices);
438 for (int j=0; j<numVertices; j++)
439 grid2ElementCorners_[i][j] = grid2_elements[grid2CornerCounter++];
440
441 }
442
444 // Compute the face neighbors for each element
446
447 computeNeighborsPerElement<grid1Dim>(grid1_element_types, grid1ElementCorners_, elementNeighbors1_);
448 computeNeighborsPerElement<grid2Dim>(grid2_element_types, grid2ElementCorners_, elementNeighbors2_);
449
450 std::cout << "setup took " << watch.elapsed() << " seconds." << std::endl;
451
452 if (m_enableBruteForce)
453 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
454 else
455 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
456
457 valid = true;
458 std::cout << "intersection construction took " << watch.elapsed() << " seconds." << std::endl;
459}
460
461template<typename T, int grid1Dim, int grid2Dim, int dimworld>
463 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
464 const std::vector<unsigned int>& grid1_elements,
465 const std::vector<Dune::GeometryType>& grid1_element_types,
466 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
467 const std::vector<unsigned int>& grid2_elements,
468 const std::vector<Dune::GeometryType>& grid2_element_types
469 )
470{
472 // Data structures for the advancing-front algorithm
474
475 std::stack<unsigned int> candidates1;
476 std::stack<unsigned int> candidates2;
477
478 std::vector<int> seeds(grid2_element_types.size(), -1);
479
480 // /////////////////////////////////////////////////////////////////////
481 // Do a brute-force search to find one pair of intersecting elements
482 // to start the advancing-front type algorithm with.
483 // /////////////////////////////////////////////////////////////////////
484
485 // Set flag if element has been handled
486 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
487
488 // Set flag if the element has been entered in the queue
489 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
490
491 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
492
493 // /////////////////////////////////////////////////////////////////////
494 // Main loop
495 // /////////////////////////////////////////////////////////////////////
496
497 std::set<unsigned int> isHandled1;
498
499 std::set<unsigned int> isCandidate1;
500
501 while (!candidates2.empty()) {
502
503 // Get the next element on the grid2 side
504 unsigned int currentCandidate2 = candidates2.top();
505 int seed = seeds[currentCandidate2];
506 assert(seed >= 0);
507
508 candidates2.pop();
509 isHandled2[currentCandidate2] = true;
510
511 // Start advancing front algorithm on the grid1 side from the 'seed' element that
512 // we stored along with the current grid2 element
513 candidates1.push(seed);
514
515 isHandled1.clear();
516 isCandidate1.clear();
517
518 while (!candidates1.empty()) {
519
520 unsigned int currentCandidate1 = candidates1.top();
521 candidates1.pop();
522 isHandled1.insert(currentCandidate1);
523
524 // Test whether there is an intersection between currentCandidate0 and currentCandidate1
525 std::bitset<(1<<grid1Dim)> neighborIntersects1;
526 std::bitset<(1<<grid2Dim)> neighborIntersects2;
527 bool intersectionFound = computeIntersection(currentCandidate1, currentCandidate2,
528 grid1Coords,grid1_element_types, neighborIntersects1,
529 grid2Coords,grid2_element_types, neighborIntersects2);
530
531 for (size_t i=0; i<neighborIntersects2.size(); i++)
532 if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
533 seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
534
535 // add neighbors of candidate0 to the list of elements to be checked
536 if (intersectionFound) {
537
538 for (size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
539
540 int neighbor = elementNeighbors1_[currentCandidate1][i];
541
542 if (neighbor == -1) // do nothing at the grid boundary
543 continue;
544
545 if (isHandled1.find(neighbor) == isHandled1.end()
546 && isCandidate1.find(neighbor) == isCandidate1.end()) {
547 candidates1.push(neighbor);
548 isCandidate1.insert(neighbor);
549 }
550
551 }
552
553 }
554
555 }
556
557 // We have now found all intersections of elements in the grid1 side with currentCandidate2
558 // Now we add all neighbors of currentCandidate2 that have not been treated yet as new
559 // candidates.
560
561 // Do we have an unhandled neighbor with a seed?
562 bool seedFound = !candidates2.empty();
563 for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
564
565 int neighbor = elementNeighbors2_[currentCandidate2][i];
566
567 if (neighbor == -1) // do nothing at the grid boundary
568 continue;
569
570 // Add all unhandled intersecting neighbors to the queue
571 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
572
573 isCandidate2[neighbor][0] = true;
574 candidates2.push(neighbor);
575 seedFound = true;
576 }
577 }
578
579 if (seedFound || !m_enableFallback)
580 continue;
581
582 // There is no neighbor with a seed, so we need to be a bit more aggressive...
583 // get all neighbors of currentCandidate2, but not currentCandidate2 itself
584 for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
585
586 int neighbor = elementNeighbors2_[currentCandidate2][i];
587
588 if (neighbor == -1) // do nothing at the grid boundary
589 continue;
590
591 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
592
593 // Get a seed element for the new grid2 element
594 // Look for an element on the grid1 side that intersects the new grid2 element.
595 int seed = -1;
596
597 // Look among the ones that have been tested during the last iteration.
598 for (typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
599 seedIt != isHandled1.end(); ++seedIt) {
600
601 std::bitset<(1<<grid1Dim)> neighborIntersects1;
602 std::bitset<(1<<grid2Dim)> neighborIntersects2;
603 bool intersectionFound = computeIntersection(*seedIt, neighbor,
604 grid1Coords, grid1_element_types, neighborIntersects1,
605 grid2Coords, grid2_element_types, neighborIntersects2,
606 false);
607
608 // if the intersection is nonempty, *seedIt is our new seed candidate on the grid1 side
609 if (intersectionFound) {
610 seed = *seedIt;
611 Dune::dwarn << "Algorithm entered first fallback method and found a new seed in the build algorithm." <<
612 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
613 break;
614 }
615
616 }
617
618 if (seed < 0) {
619 // The fast method didn't find a grid1 element that intersects with
620 // the new grid2 candidate. We have to do a brute-force search.
621 seed = bruteForceSearch(neighbor,
622 grid1Coords,grid1_element_types,
623 grid2Coords,grid2_element_types);
624 Dune::dwarn << "Algorithm entered second fallback method. This probably should not happen." << std::endl;
625
626 }
627
628 // We have tried all we could: the candidate is 'handled' now
629 isCandidate2[neighbor] = true;
630
631 // still no seed? Then the new grid2 candidate isn't overlapped by anything
632 if (seed < 0)
633 continue;
634
635 // we have a seed now
636 candidates2.push(neighbor);
637 seeds[neighbor] = seed;
638 seedFound = true;
639
640 }
641
642 }
643
644 /* Do a brute-force search if there is still no seed:
645 * There might still be a disconnected region out there.
646 */
647 if (!seedFound && candidates2.empty()) {
648 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
649 }
650 }
651}
652
653template<typename T, int grid1Dim, int grid2Dim, int dimworld>
654void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildBruteForce(
655 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
656 const std::vector<unsigned int>& grid1_elements,
657 const std::vector<Dune::GeometryType>& grid1_element_types,
658 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
659 const std::vector<unsigned int>& grid2_elements,
660 const std::vector<Dune::GeometryType>& grid2_element_types
661 )
662{
663 std::bitset<(1<<grid1Dim)> neighborIntersects1;
664 std::bitset<(1<<grid2Dim)> neighborIntersects2;
665
666 for (unsigned i = 0; i < grid1_element_types.size(); ++i) {
667 for (unsigned j = 0; j < grid2_element_types.size(); ++j) {
668 (void) computeIntersection(i, j,
669 grid1Coords, grid1_element_types, neighborIntersects1,
670 grid2Coords, grid2_element_types, neighborIntersects2);
671 }
672 }
673}
674
675template<typename T, int grid1Dim, int grid2Dim, int dimworld>
676void 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)
677{
678 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
679
680 if (seeds[j] > 0 || isHandled2[j][0])
681 continue;
682
683 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
684
685 if (seed >= 0) {
686 candidates2.push(j); // the candidate and a seed for the candidate
687 seeds[j] = seed;
688 break;
689 } else // If the brute force search did not find any intersection we can skip this element
690 isHandled2[j] = true;
691 }
692}
693
694template<typename T, int grid1Dim, int grid2Dim, int dimworld>
695int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::insertIntersections(unsigned int candidate1, unsigned int candidate2,
696 std::vector<SimplicialIntersection>& intersections)
697{
698 typedef typename std::vector<SimplicialIntersection>::size_type size_t;
699 int count = 0;
700
701 for (size_t i = 0; i < intersections.size(); ++i) {
702 // get the intersection index of the current intersection from intersections in this->intersections
703 bool found;
704 unsigned int index;
705 std::tie(found, index) = intersectionIndex(candidate1,candidate2,intersections[i]);
706
707 if (found && index >= this->intersections().size()) { //the intersection is not yet contained in this->intersections
708 this->intersections().push_back(intersections[i]); // insert
709
710 ++count;
711 } else if (found) {
712 auto& intersection = this->intersections()[index];
713
714 // insert each grid1 element and local representation of intersections[i] with parent candidate1
715 for (size_t j = 0; j < intersections[i].parents0.size(); ++j) {
716 intersection.parents0.push_back(candidate1);
717 intersection.corners0.push_back(intersections[i].corners0[j]);
718 }
719
720 // insert each grid2 element and local representation of intersections[i] with parent candidate2
721 for (size_t j = 0; j < intersections[i].parents1.size(); ++j) {
722 intersection.parents1.push_back(candidate2);
723 intersection.corners1.push_back(intersections[i].corners1[j]);
724 }
725
726 ++count;
727 } else {
728 Dune::dwarn << "Computed the same intersection twice!" << std::endl;
729 }
730 }
731 return count;
732}
733
734template<typename T, int grid1Dim, int grid2Dim, int dimworld>
735std::pair<bool, unsigned int>
736StandardMerge<T,grid1Dim,grid2Dim,dimworld>::intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
737 SimplicialIntersection& intersection) {
738
739
740 // return index in intersections_ if at least one local representation of a Simplicial Intersection (SI)
741 // of intersections_ is equal to the local representation of one element in intersections
742
743 std::size_t n_intersections = this->intersections().size();
744 if (grid1Dim == grid2Dim)
745 return {true, n_intersections};
746
747 T eps = 1e-10;
748
749 for (std::size_t i = 0; i < n_intersections; ++i) {
750
751 // compare the local representation of the subelements of the SI
752 for (std::size_t ei = 0; ei < this->intersections()[i].parents0.size(); ++ei) // merger subelement
753 {
754 if (this->intersections()[i].parents0[ei] == grid1Index)
755 {
756 for (std::size_t er = 0; er < intersection.parents0.size(); ++er) // list subelement
757 {
758 bool found_all = true;
759 // compare the local coordinate representations
760 for (std::size_t ci = 0; ci < this->intersections()[i].corners0[ei].size(); ++ci)
761 {
762 Dune::FieldVector<T,grid1Dim> ni = this->intersections()[i].corners0[ei][ci];
763 bool found_ni = false;
764 for (std::size_t cr = 0; cr < intersection.corners0[er].size(); ++cr)
765 {
766 Dune::FieldVector<T,grid1Dim> nr = intersection.corners0[er][cr];
767
768 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
769 if (found_ni)
770 break;
771 }
772 found_all = found_all && found_ni;
773
774 if (!found_ni)
775 break;
776 }
777
778 if (found_all && (this->intersections()[i].parents1[ei] != grid2Index))
779 return {true, i};
780 else if (found_all)
781 return {false, 0};
782 }
783 }
784 }
785
786 // compare the local representation of the subelements of the SI
787 for (std::size_t ei = 0; ei < this->intersections()[i].parents1.size(); ++ei) // merger subelement
788 {
789 if (this->intersections()[i].parents1[ei] == grid2Index)
790 {
791 for (std::size_t er = 0; er < intersection.parents1.size(); ++er) // list subelement
792 {
793 bool found_all = true;
794 // compare the local coordinate representations
795 for (std::size_t ci = 0; ci < this->intersections()[i].corners1[ei].size(); ++ci)
796 {
797 Dune::FieldVector<T,grid2Dim> ni = this->intersections()[i].corners1[ei][ci];
798 bool found_ni = false;
799 for (std::size_t cr = 0; cr < intersection.corners1[er].size(); ++cr)
800 {
801 Dune::FieldVector<T,grid2Dim> nr = intersection.corners1[er][cr];
802 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
803
804 if (found_ni)
805 break;
806 }
807 found_all = found_all && found_ni;
808
809 if (!found_ni)
810 break;
811 }
812
813 if (found_all && (this->intersections()[i].parents0[ei] != grid1Index))
814 return {true, i};
815 else if (found_all)
816 return {false, 0};
817 }
818 }
819 }
820 }
821
822 return {true, n_intersections};
823}
824
825#define DECL extern
826#define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \
827 DECL template \
828 void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \
829 const std::vector<unsigned int>& grid1_elements, \
830 const std::vector<Dune::GeometryType>& grid1_element_types, \
831 const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \
832 const std::vector<unsigned int>& grid2_elements, \
833 const std::vector<Dune::GeometryType>& grid2_element_types \
834 )
835
836STANDARD_MERGE_INSTANTIATE(double,1,1,1);
837STANDARD_MERGE_INSTANTIATE(double,2,2,2);
838STANDARD_MERGE_INSTANTIATE(double,3,3,3);
839#undef STANDARD_MERGE_INSTANTIATE
840#undef DECL
841
842} /* namespace GridGlue */
843} /* namespace Dune */
844
845#endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
Definition: gridglue.hh:35
STANDARD_MERGE_INSTANTIATE(double, 1, 1, 1)
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Definition: intersectionlist.hh:205
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:25
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: merger.hh:35
Dune::GridGlue::IntersectionList< Grid1Coords, Grid2Coords > IntersectionList
Definition: merger.hh:37
Dune::FieldVector< T, grid1Dim > Grid1Coords
the local coordinate type for the grid1 coordinates
Definition: merger.hh:29
Dune::FieldVector< T, grid2Dim > Grid2Coords
the local coordinate type for the grid2 coordinates
Definition: merger.hh:32
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: standardmerge.hh:56
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.
std::shared_ptr< IntersectionList > intersectionList_
Definition: standardmerge.hh:122
typename Base::Grid1Coords Grid1Coords
Type used for local coordinates on the grid1 side.
Definition: standardmerge.hh:67
void enableFallback(bool fallback)
Definition: standardmerge.hh:164
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:129
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:128
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:64
void clear() override
Definition: standardmerge.hh:148
typename Base::IntersectionList IntersectionList
Definition: standardmerge.hh:75
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:261
void enableBruteForce(bool bruteForce)
Definition: standardmerge.hh:169
std::shared_ptr< IntersectionListProvider > intersectionListProvider_
Definition: standardmerge.hh:121
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:390
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:125
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:126
typename Base::Grid2Coords Grid2Coords
Type used for local coordinates on the grid2 side.
Definition: standardmerge.hh:70
std::shared_ptr< IntersectionList > intersectionList() const final
Definition: standardmerge.hh:158
virtual ~StandardMerge()=default
SimplicialIntersection RemoteSimplicialIntersection
Definition: standardmerge.hh:82
bool valid
Definition: standardmerge.hh:84
StandardMerge()
Definition: standardmerge.hh:86
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition: standardmerge.hh:81
typename Base::WorldCoords WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:73