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