dune-grid-glue  2.7.0
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 
34 namespace Dune {
35 namespace GridGlue {
36 
53 template<class T, int grid1Dim, int grid2Dim, int dimworld>
55  : public Merger<T,grid1Dim,grid2Dim,dimworld>
56 {
58 
59 public:
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 
67  using Grid1Coords = typename Base::Grid1Coords;
68 
70  using Grid2Coords = typename Base::Grid2Coords;
71 
73  using WorldCoords = typename Base::WorldCoords;
74 
76 
77 protected:
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 
131 public:
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
151  intersectionListProvider_->clear();
152  purge(grid1ElementCorners_);
153  purge(grid2ElementCorners_);
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 
174 private:
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 
260 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
261 bool 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,
293  intersections);
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 
304 template<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 
330 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
331 template<int gridDim>
332 void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
333 computeNeighborsPerElement(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 
389 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
390 void 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 
461 template<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 
653 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
654 void 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 
675 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
676 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)
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 
694 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
695 int 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 
734 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
735 std::pair<bool, unsigned int>
736 StandardMerge<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 
836 STANDARD_MERGE_INSTANTIATE(double,1,1,1);
837 STANDARD_MERGE_INSTANTIATE(double,2,2,2);
838 STANDARD_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
Dune::GridGlue::StandardMerge::elementNeighbors2_
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:129
Dune
Definition: gridglue.hh:35
Dune::GridGlue::IntersectionList
Definition: intersectionlist.hh:132
Dune::GridGlue::StandardMerge::intersectionList
std::shared_ptr< IntersectionList > intersectionList() const final
Definition: standardmerge.hh:158
Dune::GridGlue::StandardMerge::~StandardMerge
virtual ~StandardMerge()=default
Dune::GridGlue::StandardMerge
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: standardmerge.hh:56
Dune::GridGlue::StandardMerge::grid1ElementCorners_
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:125
intersectionlist.hh
Dune::GridGlue::StandardMerge::enableFallback
void enableFallback(bool fallback)
Definition: standardmerge.hh:164
Dune::GridGlue::Merger< double, grid1Dim, grid2Dim, dimworld >::Grid1Coords
Dune::FieldVector< double, grid1Dim > Grid1Coords
the local coordinate type for the grid1 coordinates
Definition: merger.hh:29
Dune::GridGlue::SimplicialIntersectionListProvider::SimplicialIntersection
Definition: intersectionlist.hh:220
merger.hh
Dune::GridGlue::intersections
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Dune::GridGlue::StandardMerge::intersectionList_
std::shared_ptr< IntersectionList > intersectionList_
Definition: standardmerge.hh:122
Dune::GridGlue::StandardMerge::computeIntersections
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.
Dune::GridGlue::StandardMerge::grid2ElementCorners_
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:126
Dune::GridGlue::Merger< double, grid1Dim, grid2Dim, dimworld >::Grid2Coords
Dune::FieldVector< double, grid2Dim > Grid2Coords
the local coordinate type for the grid2 coordinates
Definition: merger.hh:32
Dune::GridGlue::StandardMerge::intersectionListProvider_
std::shared_ptr< IntersectionListProvider > intersectionListProvider_
Definition: standardmerge.hh:121
Dune::GridGlue::StandardMerge< double, dim, dim, dimworld >::SimplicialIntersection
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition: standardmerge.hh:81
Dune::GridGlue::StandardMerge::elementNeighbors1_
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:128
computeintersection.hh
Dune::GridGlue::Merger< T, grid1Dim, grid2Dim, dimworld >::IntersectionList
Dune::GridGlue::IntersectionList< Grid1Coords, Grid2Coords > IntersectionList
Definition: merger.hh:37
Dune::GridGlue::StandardMerge::valid
bool valid
Definition: standardmerge.hh:84
Dune::GridGlue::STANDARD_MERGE_INSTANTIATE
STANDARD_MERGE_INSTANTIATE(double, 1, 1, 1)
Dune::GridGlue::StandardMerge::StandardMerge
StandardMerge()
Definition: standardmerge.hh:86
Dune::GridGlue::StandardMerge::ctype
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:64
Dune::GridGlue::Merger< double, grid1Dim, grid2Dim, dimworld >::WorldCoords
Dune::FieldVector< double, dimworld > WorldCoords
the coordinate type used in this interface
Definition: merger.hh:35
Dune::GridGlue::StandardMerge::enableBruteForce
void enableBruteForce(bool bruteForce)
Definition: standardmerge.hh:169
Dune::GridGlue::Merger
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:25
Dune::GridGlue::StandardMerge::computeIntersection
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
Dune::GridGlue::StandardMerge< double, dim, dim, dimworld >::RemoteSimplicialIntersection
SimplicialIntersection RemoteSimplicialIntersection
Definition: standardmerge.hh:82
Dune::GridGlue::StandardMerge::clear
void clear() override
Definition: standardmerge.hh:148
Dune::GridGlue::SimplicialIntersectionListProvider
Definition: intersectionlist.hh:205
Dune::GridGlue::StandardMerge::build
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