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