dune-grid-glue  2.9
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:
3 // SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
10 #ifndef DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
11 #define DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
12 
13 
14 #include <iostream>
15 #include <iomanip>
16 #include <vector>
17 #include <stack>
18 #include <set>
19 #include <utility>
20 #include <map>
21 #include <memory>
22 #include <algorithm>
23 
24 #include <dune/common/fvector.hh>
25 #include <dune/common/bitsetvector.hh>
26 #include <dune/common/stdstreams.hh>
27 #include <dune/common/timer.hh>
28 
29 #include <dune/geometry/referenceelements.hh>
30 #include <dune/grid/common/grid.hh>
31 
35 
36 namespace Dune {
37 namespace GridGlue {
38 
55 template<class T, int grid1Dim, int grid2Dim, int dimworld>
57  : public Merger<T,grid1Dim,grid2Dim,dimworld>
58 {
60 
61 public:
62 
63  /* 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 */
64 
66  typedef T ctype;
67 
69  using Grid1Coords = typename Base::Grid1Coords;
70 
72  using Grid2Coords = typename Base::Grid2Coords;
73 
75  using WorldCoords = typename Base::WorldCoords;
76 
78 
79 protected:
80 
85 
86  bool valid = false;
87 
91  {}
92 
93  virtual ~StandardMerge() = default;
94 
99  virtual void computeIntersections(const Dune::GeometryType& grid1ElementType,
100  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
101  std::bitset<(1<<grid1Dim)>& neighborIntersects1,
102  unsigned int grid1Index,
103  const Dune::GeometryType& grid2ElementType,
104  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
105  std::bitset<(1<<grid2Dim)>& neighborIntersects2,
106  unsigned int grid2Index,
107  std::vector<SimplicialIntersection>& intersections) = 0;
108 
112  bool computeIntersection(unsigned int candidate0, unsigned int candidate1,
113  const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
114  const std::vector<Dune::GeometryType>& grid1_element_types,
115  std::bitset<(1<<grid1Dim)>& neighborIntersects1,
116  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
117  const std::vector<Dune::GeometryType>& grid2_element_types,
118  std::bitset<(1<<grid2Dim)>& neighborIntersects2,
119  bool insert = true);
120 
121  /* M E M B E R V A R I A B L E S */
122 
123  std::shared_ptr<IntersectionListProvider> intersectionListProvider_;
124  std::shared_ptr<IntersectionList> intersectionList_;
125 
127  std::vector<std::vector<unsigned int> > grid1ElementCorners_;
128  std::vector<std::vector<unsigned int> > grid2ElementCorners_;
129 
130  std::vector<std::vector<int> > elementNeighbors1_;
131  std::vector<std::vector<int> > elementNeighbors2_;
132 
133 public:
134 
135  /* 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 */
136 
140  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
141  const std::vector<unsigned int>& grid1_elements,
142  const std::vector<Dune::GeometryType>& grid1_element_types,
143  const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
144  const std::vector<unsigned int>& grid2_elements,
145  const std::vector<Dune::GeometryType>& grid2_element_types) override;
146 
147 
148  /* P R O B I N G T H E M E R G E D G R I D */
149 
150  void clear() override
151  {
152  // Delete old internal data, from a possible previous run
153  intersectionListProvider_->clear();
154  purge(grid1ElementCorners_);
155  purge(grid2ElementCorners_);
156 
157  valid = false;
158  }
159 
160  std::shared_ptr<IntersectionList> intersectionList() const final
161  {
162  assert(valid);
163  return intersectionList_;
164  }
165 
166  void enableFallback(bool fallback)
167  {
168  m_enableFallback = fallback;
169  }
170 
171  void enableBruteForce(bool bruteForce)
172  {
173  m_enableBruteForce = bruteForce;
174  }
175 
176 private:
180  bool m_enableFallback = false;
181 
185  bool m_enableBruteForce = false;
186 
187  auto& intersections()
188  { return intersectionListProvider_->intersections(); }
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 
203  void generateSeed(std::vector<int>& seeds,
204  Dune::BitSetVector<1>& isHandled2,
205  std::stack<unsigned>& candidates2,
206  const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
207  const std::vector<Dune::GeometryType>& grid1_element_types,
208  const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
209  const std::vector<Dune::GeometryType>& grid2_element_types);
210 
214  int insertIntersections(unsigned int candidate1, unsigned int candidate2,std::vector<SimplicialIntersection>& intersections);
215 
219  int bruteForceSearch(int candidate1,
220  const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
221  const std::vector<Dune::GeometryType>& grid1_element_types,
222  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
223  const std::vector<Dune::GeometryType>& grid2_element_types);
224 
228  std::pair<bool, unsigned int>
229  intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
230  SimplicialIntersection& intersection);
231 
235  template <int gridDim>
236  void computeNeighborsPerElement(const std::vector<Dune::GeometryType>& gridElementTypes,
237  const std::vector<std::vector<unsigned int> >& gridElementCorners,
238  std::vector<std::vector<int> >& elementNeighbors);
239 
240  void buildAdvancingFront(
241  const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
242  const std::vector<unsigned int>& grid1_elements,
243  const std::vector<Dune::GeometryType>& grid1_element_types,
244  const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
245  const std::vector<unsigned int>& grid2_elements,
246  const std::vector<Dune::GeometryType>& grid2_element_types
247  );
248 
249  void buildBruteForce(
250  const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
251  const std::vector<unsigned int>& grid1_elements,
252  const std::vector<Dune::GeometryType>& grid1_element_types,
253  const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
254  const std::vector<unsigned int>& grid2_elements,
255  const std::vector<Dune::GeometryType>& grid2_element_types
256  );
257 };
258 
259 
260 /* IMPLEMENTATION */
261 
262 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
263 bool StandardMerge<T,grid1Dim,grid2Dim,dimworld>::computeIntersection(unsigned int candidate0, unsigned int candidate1,
264  const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
265  const std::vector<Dune::GeometryType>& grid1_element_types,
266  std::bitset<(1<<grid1Dim)>& neighborIntersects1,
267  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
268  const std::vector<Dune::GeometryType>& grid2_element_types,
269  std::bitset<(1<<grid2Dim)>& neighborIntersects2,
270  bool insert)
271 {
272  // Select vertices of the grid1 element
273  int grid1NumVertices = grid1ElementCorners_[candidate0].size();
274  std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
275  for (int i=0; i<grid1NumVertices; i++)
276  grid1ElementCorners[i] = grid1Coords[grid1ElementCorners_[candidate0][i]];
277 
278  // Select vertices of the grid2 element
279  int grid2NumVertices = grid2ElementCorners_[candidate1].size();
280  std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
281  for (int i=0; i<grid2NumVertices; i++)
282  grid2ElementCorners[i] = grid2Coords[grid2ElementCorners_[candidate1][i]];
283 
284  // ///////////////////////////////////////////////////////
285  // Compute the intersection between the two elements
286  // ///////////////////////////////////////////////////////
287 
288  std::vector<SimplicialIntersection> intersections(0);
289 
290  // compute the intersections
291  computeIntersections(grid1_element_types[candidate0], grid1ElementCorners,
292  neighborIntersects1, candidate0,
293  grid2_element_types[candidate1], grid2ElementCorners,
294  neighborIntersects2, candidate1,
295  intersections);
296 
297  // insert intersections if needed
298  if(insert && !intersections.empty())
299  insertIntersections(candidate0,candidate1,intersections);
300 
301  // Have we found an intersection?
302  return !intersections.empty() || neighborIntersects1.any() || neighborIntersects2.any();
303 
304 }
305 
306 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
308  const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
309  const std::vector<Dune::GeometryType>& grid1_element_types,
310  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
311  const std::vector<Dune::GeometryType>& grid2_element_types)
312 {
313  std::bitset<(1<<grid1Dim)> neighborIntersects1;
314  std::bitset<(1<<grid2Dim)> neighborIntersects2;
315  for (std::size_t i=0; i<grid1_element_types.size(); i++) {
316 
317  bool intersectionFound = computeIntersection(i, candidate1,
318  grid1Coords, grid1_element_types, neighborIntersects1,
319  grid2Coords, grid2_element_types, neighborIntersects2,
320  false);
321 
322  // if there is an intersection, i is our new seed candidate on the grid1 side
323  if (intersectionFound)
324  return i;
325 
326  }
327 
328  return -1;
329 }
330 
331 
332 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
333 template<int gridDim>
334 void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
335 computeNeighborsPerElement(const std::vector<Dune::GeometryType>& gridElementTypes,
336  const std::vector<std::vector<unsigned int> >& gridElementCorners,
337  std::vector<std::vector<int> >& elementNeighbors)
338 {
339  typedef std::vector<unsigned int> FaceType;
340  typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
341 
343  // First: grid 1
345  FaceSetType faces;
346  elementNeighbors.resize(gridElementTypes.size());
347 
348  for (size_t i=0; i<gridElementTypes.size(); i++)
349  elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
350 
351  for (size_t i=0; i<gridElementTypes.size(); i++) { //iterate over all elements
352  const auto& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
353 
354  for (size_t j=0; j<(size_t)refElement.size(1); j++) { // iterate over all faces of the element
355 
356  FaceType face;
357  // extract element face
358  for (size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
359  face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
360 
361  // sort the face vertices to get rid of twists and other permutations
362  std::sort(face.begin(), face.end());
363 
364  typename FaceSetType::iterator faceHandle = faces.find(face);
365 
366  if (faceHandle == faces.end()) {
367 
368  // face has not been visited before
369  faces.insert(std::make_pair(face, std::make_pair(i,j)));
370 
371  } else {
372 
373  // face has been visited before: store the mutual neighbor information
374  elementNeighbors[i][j] = faceHandle->second.first;
375  elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
376 
377  faces.erase(faceHandle);
378 
379  }
380 
381  }
382 
383  }
384 }
385 
386 // /////////////////////////////////////////////////////////////////////
387 // Compute the intersection of all pairs of elements
388 // Linear algorithm by Gander and Japhet, Proc. of DD18
389 // /////////////////////////////////////////////////////////////////////
390 
391 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
392 void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
393  const std::vector<unsigned int>& grid1_elements,
394  const std::vector<Dune::GeometryType>& grid1_element_types,
395  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
396  const std::vector<unsigned int>& grid2_elements,
397  const std::vector<Dune::GeometryType>& grid2_element_types
398  )
399 {
400 
401  std::cout << "StandardMerge building merged grid..." << std::endl;
402  Dune::Timer watch;
403 
404  clear();
405  // clear global intersection list
406  intersectionListProvider_->clear();
407  this->counter = 0;
408 
409  // /////////////////////////////////////////////////////////////////////
410  // Copy element corners into a data structure with block-structure.
411  // This is not as efficient but a lot easier to use.
412  // We may think about efficiency later.
413  // /////////////////////////////////////////////////////////////////////
414 
415  // first the grid1 side
416  grid1ElementCorners_.resize(grid1_element_types.size());
417 
418  unsigned int grid1CornerCounter = 0;
419 
420  for (std::size_t i=0; i<grid1_element_types.size(); i++) {
421 
422  // Select vertices of the grid1 element
423  int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
424  grid1ElementCorners_[i].resize(numVertices);
425  for (int j=0; j<numVertices; j++)
426  grid1ElementCorners_[i][j] = grid1_elements[grid1CornerCounter++];
427 
428  }
429 
430  // then the grid2 side
431  grid2ElementCorners_.resize(grid2_element_types.size());
432 
433  unsigned int grid2CornerCounter = 0;
434 
435  for (std::size_t i=0; i<grid2_element_types.size(); i++) {
436 
437  // Select vertices of the grid2 element
438  int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
439  grid2ElementCorners_[i].resize(numVertices);
440  for (int j=0; j<numVertices; j++)
441  grid2ElementCorners_[i][j] = grid2_elements[grid2CornerCounter++];
442 
443  }
444 
446  // Compute the face neighbors for each element
448 
449  computeNeighborsPerElement<grid1Dim>(grid1_element_types, grid1ElementCorners_, elementNeighbors1_);
450  computeNeighborsPerElement<grid2Dim>(grid2_element_types, grid2ElementCorners_, elementNeighbors2_);
451 
452  std::cout << "setup took " << watch.elapsed() << " seconds." << std::endl;
453 
454  if (m_enableBruteForce)
455  buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
456  else
457  buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
458 
459  valid = true;
460  std::cout << "intersection construction took " << watch.elapsed() << " seconds." << std::endl;
461 }
462 
463 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
465  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 {
474  // Data structures for the advancing-front algorithm
476 
477  std::stack<unsigned int> candidates1;
478  std::stack<unsigned int> candidates2;
479 
480  std::vector<int> seeds(grid2_element_types.size(), -1);
481 
482  // /////////////////////////////////////////////////////////////////////
483  // Do a brute-force search to find one pair of intersecting elements
484  // to start the advancing-front type algorithm with.
485  // /////////////////////////////////////////////////////////////////////
486 
487  // Set flag if element has been handled
488  Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
489 
490  // Set flag if the element has been entered in the queue
491  Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
492 
493  generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
494 
495  // /////////////////////////////////////////////////////////////////////
496  // Main loop
497  // /////////////////////////////////////////////////////////////////////
498 
499  std::set<unsigned int> isHandled1;
500 
501  std::set<unsigned int> isCandidate1;
502 
503  while (!candidates2.empty()) {
504 
505  // Get the next element on the grid2 side
506  unsigned int currentCandidate2 = candidates2.top();
507  int seed = seeds[currentCandidate2];
508  assert(seed >= 0);
509 
510  candidates2.pop();
511  isHandled2[currentCandidate2] = true;
512 
513  // Start advancing front algorithm on the grid1 side from the 'seed' element that
514  // we stored along with the current grid2 element
515  candidates1.push(seed);
516 
517  isHandled1.clear();
518  isCandidate1.clear();
519 
520  while (!candidates1.empty()) {
521 
522  unsigned int currentCandidate1 = candidates1.top();
523  candidates1.pop();
524  isHandled1.insert(currentCandidate1);
525 
526  // Test whether there is an intersection between currentCandidate0 and currentCandidate1
527  std::bitset<(1<<grid1Dim)> neighborIntersects1;
528  std::bitset<(1<<grid2Dim)> neighborIntersects2;
529  bool intersectionFound = computeIntersection(currentCandidate1, currentCandidate2,
530  grid1Coords,grid1_element_types, neighborIntersects1,
531  grid2Coords,grid2_element_types, neighborIntersects2);
532 
533  for (size_t i=0; i<neighborIntersects2.size(); i++)
534  if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
535  seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
536 
537  // add neighbors of candidate0 to the list of elements to be checked
538  if (intersectionFound) {
539 
540  for (size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
541 
542  int neighbor = elementNeighbors1_[currentCandidate1][i];
543 
544  if (neighbor == -1) // do nothing at the grid boundary
545  continue;
546 
547  if (isHandled1.find(neighbor) == isHandled1.end()
548  && isCandidate1.find(neighbor) == isCandidate1.end()) {
549  candidates1.push(neighbor);
550  isCandidate1.insert(neighbor);
551  }
552 
553  }
554 
555  }
556 
557  }
558 
559  // We have now found all intersections of elements in the grid1 side with currentCandidate2
560  // Now we add all neighbors of currentCandidate2 that have not been treated yet as new
561  // candidates.
562 
563  // Do we have an unhandled neighbor with a seed?
564  bool seedFound = !candidates2.empty();
565  for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
566 
567  int neighbor = elementNeighbors2_[currentCandidate2][i];
568 
569  if (neighbor == -1) // do nothing at the grid boundary
570  continue;
571 
572  // Add all unhandled intersecting neighbors to the queue
573  if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
574 
575  isCandidate2[neighbor][0] = true;
576  candidates2.push(neighbor);
577  seedFound = true;
578  }
579  }
580 
581  if (seedFound || !m_enableFallback)
582  continue;
583 
584  // There is no neighbor with a seed, so we need to be a bit more aggressive...
585  // get all neighbors of currentCandidate2, but not currentCandidate2 itself
586  for (size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
587 
588  int neighbor = elementNeighbors2_[currentCandidate2][i];
589 
590  if (neighbor == -1) // do nothing at the grid boundary
591  continue;
592 
593  if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
594 
595  // Get a seed element for the new grid2 element
596  // Look for an element on the grid1 side that intersects the new grid2 element.
597  int seed = -1;
598 
599  // Look among the ones that have been tested during the last iteration.
600  for (typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
601  seedIt != isHandled1.end(); ++seedIt) {
602 
603  std::bitset<(1<<grid1Dim)> neighborIntersects1;
604  std::bitset<(1<<grid2Dim)> neighborIntersects2;
605  bool intersectionFound = computeIntersection(*seedIt, neighbor,
606  grid1Coords, grid1_element_types, neighborIntersects1,
607  grid2Coords, grid2_element_types, neighborIntersects2,
608  false);
609 
610  // if the intersection is nonempty, *seedIt is our new seed candidate on the grid1 side
611  if (intersectionFound) {
612  seed = *seedIt;
613  Dune::dwarn << "Algorithm entered first fallback method and found a new seed in the build algorithm." <<
614  "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
615  break;
616  }
617 
618  }
619 
620  if (seed < 0) {
621  // The fast method didn't find a grid1 element that intersects with
622  // the new grid2 candidate. We have to do a brute-force search.
623  seed = bruteForceSearch(neighbor,
624  grid1Coords,grid1_element_types,
625  grid2Coords,grid2_element_types);
626  Dune::dwarn << "Algorithm entered second fallback method. This probably should not happen." << std::endl;
627 
628  }
629 
630  // We have tried all we could: the candidate is 'handled' now
631  isCandidate2[neighbor] = true;
632 
633  // still no seed? Then the new grid2 candidate isn't overlapped by anything
634  if (seed < 0)
635  continue;
636 
637  // we have a seed now
638  candidates2.push(neighbor);
639  seeds[neighbor] = seed;
640  seedFound = true;
641 
642  }
643 
644  }
645 
646  /* Do a brute-force search if there is still no seed:
647  * There might still be a disconnected region out there.
648  */
649  if (!seedFound && candidates2.empty()) {
650  generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
651  }
652  }
653 }
654 
655 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
656 void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildBruteForce(
657  const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
658  const std::vector<unsigned int>& grid1_elements,
659  const std::vector<Dune::GeometryType>& grid1_element_types,
660  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
661  const std::vector<unsigned int>& grid2_elements,
662  const std::vector<Dune::GeometryType>& grid2_element_types
663  )
664 {
665  std::bitset<(1<<grid1Dim)> neighborIntersects1;
666  std::bitset<(1<<grid2Dim)> neighborIntersects2;
667 
668  for (unsigned i = 0; i < grid1_element_types.size(); ++i) {
669  for (unsigned j = 0; j < grid2_element_types.size(); ++j) {
670  (void) computeIntersection(i, j,
671  grid1Coords, grid1_element_types, neighborIntersects1,
672  grid2Coords, grid2_element_types, neighborIntersects2);
673  }
674  }
675 }
676 
677 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
678 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)
679 {
680  for (std::size_t j=0; j<grid2_element_types.size(); j++) {
681 
682  if (seeds[j] > 0 || isHandled2[j][0])
683  continue;
684 
685  int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
686 
687  if (seed >= 0) {
688  candidates2.push(j); // the candidate and a seed for the candidate
689  seeds[j] = seed;
690  break;
691  } else // If the brute force search did not find any intersection we can skip this element
692  isHandled2[j] = true;
693  }
694 }
695 
696 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
697 int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::insertIntersections(unsigned int candidate1, unsigned int candidate2,
698  std::vector<SimplicialIntersection>& intersections)
699 {
700  typedef typename std::vector<SimplicialIntersection>::size_type size_t;
701  int count = 0;
702 
703  for (size_t i = 0; i < intersections.size(); ++i) {
704  // get the intersection index of the current intersection from intersections in this->intersections
705  bool found;
706  unsigned int index;
707  std::tie(found, index) = intersectionIndex(candidate1,candidate2,intersections[i]);
708 
709  if (found && index >= this->intersections().size()) { //the intersection is not yet contained in this->intersections
710  this->intersections().push_back(intersections[i]); // insert
711 
712  ++count;
713  } else if (found) {
714  auto& intersection = this->intersections()[index];
715 
716  // insert each grid1 element and local representation of intersections[i] with parent candidate1
717  for (size_t j = 0; j < intersections[i].parents0.size(); ++j) {
718  intersection.parents0.push_back(candidate1);
719  intersection.corners0.push_back(intersections[i].corners0[j]);
720  }
721 
722  // insert each grid2 element and local representation of intersections[i] with parent candidate2
723  for (size_t j = 0; j < intersections[i].parents1.size(); ++j) {
724  intersection.parents1.push_back(candidate2);
725  intersection.corners1.push_back(intersections[i].corners1[j]);
726  }
727 
728  ++count;
729  } else {
730  Dune::dwarn << "Computed the same intersection twice!" << std::endl;
731  }
732  }
733  return count;
734 }
735 
736 template<typename T, int grid1Dim, int grid2Dim, int dimworld>
737 std::pair<bool, unsigned int>
738 StandardMerge<T,grid1Dim,grid2Dim,dimworld>::intersectionIndex(unsigned int grid1Index, unsigned int grid2Index,
739  SimplicialIntersection& intersection) {
740 
741 
742  // return index in intersections_ if at least one local representation of a Simplicial Intersection (SI)
743  // of intersections_ is equal to the local representation of one element in intersections
744 
745  std::size_t n_intersections = this->intersections().size();
746  if (grid1Dim == grid2Dim)
747  return {true, n_intersections};
748 
749  T eps = 1e-10;
750 
751  for (std::size_t i = 0; i < n_intersections; ++i) {
752 
753  // compare the local representation of the subelements of the SI
754  for (std::size_t ei = 0; ei < this->intersections()[i].parents0.size(); ++ei) // merger subelement
755  {
756  if (this->intersections()[i].parents0[ei] == grid1Index)
757  {
758  for (std::size_t er = 0; er < intersection.parents0.size(); ++er) // list subelement
759  {
760  bool found_all = true;
761  // compare the local coordinate representations
762  for (std::size_t ci = 0; ci < this->intersections()[i].corners0[ei].size(); ++ci)
763  {
764  Dune::FieldVector<T,grid1Dim> ni = this->intersections()[i].corners0[ei][ci];
765  bool found_ni = false;
766  for (std::size_t cr = 0; cr < intersection.corners0[er].size(); ++cr)
767  {
768  Dune::FieldVector<T,grid1Dim> nr = intersection.corners0[er][cr];
769 
770  found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
771  if (found_ni)
772  break;
773  }
774  found_all = found_all && found_ni;
775 
776  if (!found_ni)
777  break;
778  }
779 
780  if (found_all && (this->intersections()[i].parents1[ei] != grid2Index))
781  return {true, i};
782  else if (found_all)
783  return {false, 0};
784  }
785  }
786  }
787 
788  // compare the local representation of the subelements of the SI
789  for (std::size_t ei = 0; ei < this->intersections()[i].parents1.size(); ++ei) // merger subelement
790  {
791  if (this->intersections()[i].parents1[ei] == grid2Index)
792  {
793  for (std::size_t er = 0; er < intersection.parents1.size(); ++er) // list subelement
794  {
795  bool found_all = true;
796  // compare the local coordinate representations
797  for (std::size_t ci = 0; ci < this->intersections()[i].corners1[ei].size(); ++ci)
798  {
799  Dune::FieldVector<T,grid2Dim> ni = this->intersections()[i].corners1[ei][ci];
800  bool found_ni = false;
801  for (std::size_t cr = 0; cr < intersection.corners1[er].size(); ++cr)
802  {
803  Dune::FieldVector<T,grid2Dim> nr = intersection.corners1[er][cr];
804  found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
805 
806  if (found_ni)
807  break;
808  }
809  found_all = found_all && found_ni;
810 
811  if (!found_ni)
812  break;
813  }
814 
815  if (found_all && (this->intersections()[i].parents0[ei] != grid1Index))
816  return {true, i};
817  else if (found_all)
818  return {false, 0};
819  }
820  }
821  }
822  }
823 
824  return {true, n_intersections};
825 }
826 
827 #define DECL extern
828 #define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \
829  DECL template \
830  void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \
831  const std::vector<unsigned int>& grid1_elements, \
832  const std::vector<Dune::GeometryType>& grid1_element_types, \
833  const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \
834  const std::vector<unsigned int>& grid2_elements, \
835  const std::vector<Dune::GeometryType>& grid2_element_types \
836  )
837 
838 STANDARD_MERGE_INSTANTIATE(double,1,1,1);
839 STANDARD_MERGE_INSTANTIATE(double,2,2,2);
840 STANDARD_MERGE_INSTANTIATE(double,3,3,3);
841 #undef STANDARD_MERGE_INSTANTIATE
842 #undef DECL
843 
844 } /* namespace GridGlue */
845 } /* namespace Dune */
846 
847 #endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
Definition: gridglue.hh:37
STANDARD_MERGE_INSTANTIATE(double, 1, 1, 1)
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Definition: intersectionlist.hh:207
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:27
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: merger.hh:37
Dune::GridGlue::IntersectionList< Grid1Coords, Grid2Coords > IntersectionList
Definition: merger.hh:39
Dune::FieldVector< T, grid1Dim > Grid1Coords
the local coordinate type for the grid1 coordinates
Definition: merger.hh:31
Dune::FieldVector< T, grid2Dim > Grid2Coords
the local coordinate type for the grid2 coordinates
Definition: merger.hh:34
Common base class for many merger implementations: produce pairs of entities that may intersect.
Definition: standardmerge.hh:58
virtual void computeIntersections(const Dune::GeometryType &grid1ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid1ElementCorners, std::bitset<(1<< grid1Dim)> &neighborIntersects1, unsigned int grid1Index, const Dune::GeometryType &grid2ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid2ElementCorners, std::bitset<(1<< grid2Dim)> &neighborIntersects2, unsigned int grid2Index, std::vector< SimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
std::shared_ptr< IntersectionList > intersectionList_
Definition: standardmerge.hh:124
typename Base::Grid1Coords Grid1Coords
Type used for local coordinates on the grid1 side.
Definition: standardmerge.hh:69
void enableFallback(bool fallback)
Definition: standardmerge.hh:166
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:131
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:130
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:66
void clear() override
Definition: standardmerge.hh:150
typename Base::IntersectionList IntersectionList
Definition: standardmerge.hh:77
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:263
void enableBruteForce(bool bruteForce)
Definition: standardmerge.hh:171
std::shared_ptr< IntersectionList > intersectionList() const final
Definition: standardmerge.hh:160
std::shared_ptr< IntersectionListProvider > intersectionListProvider_
Definition: standardmerge.hh:123
void build(const std::vector< Dune::FieldVector< T, dimworld > > &grid1_Coords, const std::vector< unsigned int > &grid1_elements, const std::vector< Dune::GeometryType > &grid1_element_types, const std::vector< Dune::FieldVector< T, dimworld > > &grid2_coords, const std::vector< unsigned int > &grid2_elements, const std::vector< Dune::GeometryType > &grid2_element_types) override
builds the merged grid
Definition: standardmerge.hh:392
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:127
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:128
typename Base::Grid2Coords Grid2Coords
Type used for local coordinates on the grid2 side.
Definition: standardmerge.hh:72
virtual ~StandardMerge()=default
SimplicialIntersection RemoteSimplicialIntersection
Definition: standardmerge.hh:84
bool valid
Definition: standardmerge.hh:86
StandardMerge()
Definition: standardmerge.hh:88
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition: standardmerge.hh:83
typename Base::WorldCoords WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:75