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