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