dune-alugrid  2.8-git
bisectioncompatibility.hh
Go to the documentation of this file.
1 #ifndef ALUGRID_COMPATIBILITY_CHECK_HH_INCLUDED
2 #define ALUGRID_COMPATIBILITY_CHECK_HH_INCLUDED
3 
4 #include <iostream>
5 #include <array>
6 #include <map>
7 #include <list>
8 #include <vector>
9 #include <algorithm>
10 #include <assert.h>
11 #include <random>
12 
13 #include <dune/common/timer.hh>
14 
16 {
17  static int& variant ()
18  {
19  static int var = 0; // default is every vertex in set V0
20  return var;
21  }
22 
23  static int& threshold ()
24  {
25  static int var = 2; // default is 2 for threshold
26  return var;
27  }
28 
29  static int& useAnnouncedEdge ()
30  {
31  static int var = 1 ; // default is use announced edges
32  return var;
33  }
34 };
35 
36 //Class to correct the element orientation to make bisection work in 3d
37 // It provides different algorithms to orientate a grid.
38 // Also implements checks for compatibility.
39 template <class VertexVector>
41 {
43 public:
44  // type of vertex coordinates stored inside the factory
45  typedef VertexVector VertexVectorType;
46 
47  typedef std::array<unsigned int, 3> FaceType;
48  typedef std::vector< unsigned int > ElementType;
49  typedef std::array<unsigned int, 2> EdgeType;
50  typedef std::map< FaceType, EdgeType > FaceMapType;
51  typedef std::pair< FaceType, EdgeType > FaceElementType;
52 
53 protected:
55 
56  //the elements to be renumbered
57  std::vector<ElementType> elements_;
58  std::vector<bool> elementOrientation_;
59  //the neighbouring structure
61  // the number of vertices
62  const size_t nVertices_;
63  //A tag vector for vertices to
64  //decide whether they are in V_1 or v_0
65  std::vector<bool> containedInV0_;
66  //the element types
67  std::vector<int> types_;
68  //true if stevenson notation is used
69  //false for ALBERTA
71 
72  //the 2 nodes of the refinement edge
73  EdgeType type0nodes_; // = stevensonRefinement_ ? 0,3 : 0,1 ;
74  //the faces opposite of type 0 nodes
76  //The interior node of a type 1 element
77  unsigned int type1node_; // = stevensonRefinement_ ? 1 : 2;
78  //the face opposite of the interior node
79  unsigned int type1face_; // = 3 - type1node_ ;
80 
81  // 0 = put all vertices in V0,
82  // 1 = longest edge,
83  // 2 = least adjacent elements,
84  // 3 = random ordering
88 
89 public:
90  //constructor taking elements
91  //assumes standard orientation elemIndex % 2
93  const std::vector<ElementType>& elements,
94  const bool stevenson)
95  : vertices_( vertices ),
96  elements_( elements ),
97  elementOrientation_(elements_.size(), true),
98  nVertices_( vertices_.size() ),
100  types_(elements_.size(), 0),
101  stevensonRefinement_(stevenson),
104  type1node_( stevensonRefinement_ ? 1 : 2 ),
105  type1face_( 3 - type1node_ )
106  {
107  //build the information about neighbours
108  Dune::Timer timer;
109  buildNeighbors();
110  std::cout << "Build neighbors took " << timer.elapsed() << " sec." << std::endl;
111  }
112 
113  //check for strong compatibility
115  {
116  int result = 0;
117  bool verbose = false;
118  unsigned int bndFaces = 0;
119  std::vector<int> nonCompatFacesAtVertex(nVertices_, 0 );
120  for(auto&& face : neighbours_)
121  {
122  if( face.second[0] == face.second[1] )
123  {
124  bndFaces++;
125  }
126  else if(! checkStrongCompatibility(face, verbose))
127  {
128  ++result;
129  for(size_t i =0; i < 3; ++i)
130  {
131  ++(nonCompatFacesAtVertex[face.first[i]]);
132  }
133  }
134  }
135  std::cout << "NotStrongCompatibleMacroFaces" << " InnerFaces " << " TotalFaces " << "Maximum/Vertex " << " Minimum/Vertex " << std::endl;
136  std::cout << result << " " << neighbours_.size() - bndFaces << " " << neighbours_.size() << " " << *(std::max_element(nonCompatFacesAtVertex.begin(), nonCompatFacesAtVertex.end())) << " " << *(std::min_element(nonCompatFacesAtVertex.begin(), nonCompatFacesAtVertex.end())) << std::endl << std::endl;
137  return result;
138  }
139 
140  //check grid for compatibility
141  //i.e. walk over all faces and check their compatibility.
143  {
144  for(auto&& face : neighbours_)
145  {
146  if(!checkFaceCompatibility(face)) return false;
147  }
148  return true;
149  }
150 
152  {
153  //set types to 0, and switch vertices 2,3 for elemIndex % 2
154  applyStandardOrientation();
155  bool result = compatibilityCheck();
156  //set types to 0, and switch vertices 2,3 for elemIndex % 2
157  applyStandardOrientation();
158  return result;
159  }
160 
161  //print the neighbouring structure
162  void printNB()
163  {
164  for(auto&& face : neighbours_)
165  {
166  std::cout << face.first[0] << "," << face.first[1] << "," << face.first[2] << " : " << face.second[0] << ", " << face.second[1] << std::endl;
167  }
168  }
169 
170  //print an element with orientation and all refinement edges
171  void printElement(int index)
172  {
173  ElementType el = elements_[index];
174  EdgeType edge;
175  std::cout << "[" << el[0] ;
176  for(int i=1; i < 4; ++i)
177  std::cout << ","<< el[i] ;
178  std::cout << "] Refinement Edges: ";
179  for(int i=0; i< 4; ++i)
180  {
181  getRefinementEdge(el, i, edge, types_[index]);
182  std::cout << "[" << edge[0] << "," << edge[1] << "] ";
183  }
184  std::cout << " Type: " << types_[ index ] << " V1: ";
185  for (size_t i = 0; i < 4; ++i) {
186  if( ! containedInV0_[ el [ i ] ] )
187  std::cout << el[i] << "," ;
188  }
189  std::cout << std::endl;
190  }
191 
192  //an algorithm using only elements of type 0
193  //it works by sorting the vertices in a global ordering
194  //and tries to make as many reflected neighbors as possible.
196  {
197  // convert ordering of refinement edge to stevenson
199 
200  //calculate the sets V0 and V1
201  calculateV0( variant_, threshold_ );
202  const bool useAnnounced = useAnnouncedEdge_;
203 
204  // all elements are type 0
205  std::fill( types_.begin(), types_.end(), 0 );
206 
207  std::list<std::pair<FaceType, EdgeType> > activeFaceList; // use std::find to find
208  std::vector<bool> doneElements(elements_.size(), false);
209 
210  std::vector< std::pair< double, std::pair< int,int > > > vertexOrder( nVertices_ , std::make_pair(-1.0, std::make_pair(-1,-2) ) );
211  const double eps = std::numeric_limits< double >::epsilon() * double(nVertices_) * 10.0;
212 
213 
214  const unsigned int l1 = -1;
215  //create the vertex priority List
216  const int numberOfElements = elements_.size();
217  Dune::Timer timer;
218  for(int counter = 0; counter < numberOfElements ; ++counter)
219  {
220  FaceElementType faceElement = std::make_pair( FaceType{l1,l1,l1}, EdgeType{l1,l1} );
221  if(counter == 0)
222  {
223  FaceType face;
224  const ElementType& el0 = elements_[0];
225  getFace(el0, type0faces_[0], face);
226  faceElement = FaceElementType(std::make_pair( face , EdgeType{0,0} ) );
227 
228  //orientate E_0 (add vertices to vertexPriorityList)
229  for(unsigned int i=0 ; i < 4 ; ++i)
230  {
231  int vtx = el0[ i ];
232  vertexOrder[ vtx ].first = i+1;
233  // previous vertex index or -1
234  vertexOrder[ vtx ].second.first = i > 0 ? el0[ i-1 ] : -1;
235  vertexOrder[ vtx ].second.second = i < 3 ? el0[ i+1 ] : -2;
236  }
237  }
238  else
239  {
240  assert( ! activeFaceList.empty() );
241  //take element at first face from activeFaceList
242  auto it = activeFaceList.begin();
243  int el1 = it->second[ 0 ];
244  int el2 = it->second[ 1 ];
245 
246  while( doneElements[ el1 ] && doneElements[ el2 ] )
247  {
248  activeFaceList.erase( it++ );
249  /*
250  if( it == activeFaceList.end() )
251  {
252  FaceType face;
253  const ElementType& el0 = elements_[0];
254  getFace(el0, type0faces_[0], face);
255  faceElement = FaceElementType(std::make_pair( face , EdgeType( {0,0} ) ) );
256 
257  //orientate E_0 (add vertices to vertexPriorityList)
258  for(unsigned int i=0 ; i < 4 ; ++i)
259  {
260  int vtx = el0[ i ];
261  vertexOrder[ vtx ].first = i+1;
262  // previous vertex index or -1
263  vertexOrder[ vtx ].second.first = i > 0 ? el0[ i-1 ] : -1;
264  vertexOrder[ vtx ].second.second = i < 3 ? el0[ i+1 ] : -2;
265  }
266  }
267  */
268 
269  assert( it != activeFaceList.end() );
270  el1 = it->second[ 0 ];
271  el2 = it->second[ 1 ];
272  }
273 
274  faceElement = *it;
275  }
276 
277  //el has been worked on
278  int elIndex = faceElement.second[0];
279 
280  //neigh is to be worked on
281  int neighIndex = faceElement.second[1];
282  if(!doneElements[elIndex])
283  {
284  assert(doneElements[neighIndex] || counter == 0 );
285  std::swap(elIndex, neighIndex);
286  }
287 
288  assert(!doneElements[neighIndex] || counter == 0);
289  doneElements[neighIndex] = true;
290  ElementType el = elements_[elIndex];
291  ElementType & neigh = elements_[neighIndex];
292  unsigned int faceInEl = getFaceIndex(el, faceElement.first);
293  unsigned int nodeInEl = 3 - faceInEl;
294  unsigned int faceInNeigh = getFaceIndex(neigh, faceElement.first);
295  unsigned int nodeInNeigh = 3 - faceInNeigh;
296 
297  //helper element that contains neigh
298  // in the ordering of vertexPriorityList
299  ElementType newNeigh(el);
300 
301  //insertion of new vertex before nodeInEl
302  if( vertexOrder[ neigh [ nodeInNeigh ] ].first < 0 )
303  {
304  int vxIndex = el[ nodeInEl ];
305 
306  //this takes care that the children will be reflected neighbors
307  //if nodeInNeigh = 3 && nodeInEl = 0 insert after el[3]
308  if( useAnnounced && (nodeInEl == type0nodes_[0] && nodeInNeigh == type0nodes_[1] ) )
309  {
310  auto& vxPair = vertexOrder[ el[ nodeInNeigh ] ];
311  // got to next vertex
312  vxIndex = vxPair.second.second;
313  if( vxIndex == -2 )
314  {
315  vertexOrder[ neigh [ nodeInNeigh ] ] = std::make_pair( vxPair.first+1.0, std::make_pair( el[ nodeInNeigh ], -2 ) );
316  vxPair.second.second = neigh [ nodeInNeigh ];
317  }
318  }
319 
320  if( vxIndex >= 0 )
321  {
322  //if nodeInNeigh = 0 && nodeInEl = 3 insert before el[0]
323  if ( useAnnounced && ( nodeInEl == type0nodes_[1] && nodeInNeigh == type0nodes_[0] ) )
324  {
325  vxIndex = el[ nodeInNeigh ];
326  }
327 
328  auto& vxPair = vertexOrder[ vxIndex ];
329  assert( vxPair.first >= 0 );
330 
331  // new vertex weight is average between previous and this one
332  const int prevIdx = vxPair.second.first ;
333  double prevOrder = ( prevIdx == -1 ) ? 0.0 : vertexOrder[ prevIdx ].first;
334  double newOrder = 0.5 * (vxPair.first + prevOrder);
335 
336  vertexOrder[ neigh [ nodeInNeigh ] ] = std::make_pair( newOrder, std::make_pair( prevIdx, vxIndex ) );
337  if( prevIdx >= 0 )
338  vertexOrder[ prevIdx ].second.second = neigh [ nodeInNeigh ];
339  vxPair.second.first = neigh [ nodeInNeigh ];
340  assert( vxPair.first > newOrder );
341 
342  if( (vxPair.first - newOrder) < eps )
343  {
344 #ifndef NDEBUG
345  Dune::Timer restimer;
346  std::cout << "Rescale vertex order weights." << std::endl;
347 #endif
348  std::map< double, int > newVertexOrder;
349  const int size = vertexOrder.size();
350  for( int j=0; j<size; ++j )
351  {
352  if( vertexOrder[ j ].first > 0.0 )
353  {
354  newVertexOrder.insert( std::make_pair( vertexOrder[ j ].first, j ) );
355  }
356  }
357 
358  int count = 1;
359  for( const auto& vx: newVertexOrder )
360  {
361  assert( vx.second >=0 && vx.second < size );
362  vertexOrder[ vx.second ].first = count++ ;
363  }
364 
365  for( int j=0; j<size; ++j )
366  {
367  if( vertexOrder[ j ].first > 0.0 && vertexOrder[ j ].second.first >= 0 )
368  {
369  if( vertexOrder[ j ].first <= vertexOrder[ vertexOrder[ j ].second.first ].first )
370  std::abort();
371  }
372  }
373 #ifndef NDEBUG
374  std::cout << "Rescale done, time = " << restimer.elapsed() << std::endl;
375 #endif
376  }
377  }
378  }
379 
380  {
381  // sort vertices in ascending order
382  std::map< double, int > vx;
383  for( int j=0; j<4; ++j )
384  {
385  assert( vertexOrder[ neigh[ j ] ].first > 0 );
386  vx.insert( std::make_pair( vertexOrder[ neigh[ j ] ].first, j ) );
387  }
388 
389  int count = 0;
390  for( const auto& vxItem : vx )
391  {
392  newNeigh[ count++ ] = neigh[ vxItem.second ] ;
393  }
394  }
395 
396  //adjust type of newNeigh
397  unsigned int type = 0;
398  bool contained[4];
399  for(unsigned int i =0; i < 4; ++i)
400  {
401  contained[i] = containedInV0_[newNeigh[i]];
402  }
403  //if all are in V0 or all are V1
404  //we do nothing, set type to 0 and keep the order
405  if( ( contained[0] && contained[1] && contained[2] && contained[3] ) || ( !contained[0] && !contained[1] && !contained[2] && !contained[3] ) )
406  {
407  //do nothing
408  ;
409  }
410  else
411  {
412  ElementType V0Part(newNeigh);
413  ElementType V1Part(newNeigh);
414  for(unsigned int i = 0 ; i < 4; ++i)
415  {
416  if( contained[ i ] )
417  {
418  auto it = std::find(V1Part.begin(),V1Part.end(),newNeigh[i]);
419  V1Part.erase( it );
420  }
421  else
422  {
423  auto it = std::find(V0Part.begin(),V0Part.end(),newNeigh[i]);
424  V0Part.erase( it );
425  ++type;
426  }
427  }
428  for(unsigned int i = 0; i < 4; ++i)
429  {
430  if(i == 0)
431  newNeigh[ i ] = V0Part[ i ];
432  else if( i <= type )
433  newNeigh[ i ] = V1Part[ i - 1 ] ;
434  else if( i > type)
435  newNeigh[ i ] = V0Part[ i - type];
436  }
437  }
438  types_[neighIndex] = type % 3;
439 
440  //reorientate neigh using the helper element newNeigh
441  //we use swaps to be able to track the elementOrientation_
442  bool neighOrientation = elementOrientation_[neighIndex];
443  for(unsigned int i =0 ; i < 3; ++i)
444  {
445  if( newNeigh[i] != neigh[i] )
446  {
447  auto neighIt = std::find(neigh.begin() + i,neigh.end(),newNeigh[i]);
448  std::swap(*neighIt,neigh[i]);
449  neighOrientation = ! neighOrientation;
450  }
451  }
452  elementOrientation_[neighIndex] = neighOrientation;
453  //add and remove faces from activeFaceList
454  for(unsigned int i = 0; i < 4 ; ++i)
455  {
456  getFace(neigh, i, faceElement);
457  //do nothing for boundary
458  if(faceElement.second[0] == faceElement.second[1])
459  continue;
460 
461  //if face does not contain ref Edge
462  if(i != type0faces_[0] && i != type0faces_[1] )
463  activeFaceList.push_back(faceElement);
464  else
465  activeFaceList.push_front(faceElement);
466  }
467 
468 #ifndef NDEBUG
469  const int onePercent = numberOfElements / 100 ;
470  if( onePercent > 0 && counter % onePercent == 0 )
471  {
472  std::cout << "Done: element " << counter << " of " << numberOfElements << " time used = " << timer.elapsed() << std::endl;
473  timer.reset();
474  }
475 #endif
476  }// end elements ?
477 
478  return compatibilityCheck();
479  }
480 
481 
482  //An algorithm using only elements of type 1
484  {
485  // convert to stevenson ordering
487 
488  //set all types to 1
489  std::fill( types_.begin(), types_.end(), 1 );
490 
491  //the currently active Faces.
492  //and the free faces that can still be adjusted at the end.
493  FaceMapType activeFaces, freeFaces;
494  //the finished elements. The number indicates the fixed node
495  //if it is -1, the element has not been touched yet.
496  std::vector<int> nodePriority;
497  nodePriority.resize(nVertices_ , -1);
498  int currNodePriority =nVertices_;
499 
500  const unsigned int numberOfElements = elements_.size();
501  //walk over all elements
502  for(unsigned int elIndex =0 ; elIndex < numberOfElements; ++elIndex)
503  {
504  //if no node is constrained and no face is active, fix one (e.g. smallest index)
505  ElementType & el = elements_[elIndex];
506  int priorityNode = -1;
507  FaceType face;
508  int freeNode = -1;
509  for(int i = 0; i < 4; ++i)
510  {
511  int tmpPrio = nodePriority[el[i]];
512  //if a node has positive priority
513  if( tmpPrio > -1 )
514  {
515  if( priorityNode < 0 || tmpPrio > nodePriority[el[priorityNode]] )
516  priorityNode = i;
517  }
518  getFace(el,3 - i,face );
519  //if we have a free face, the opposite node is good to be fixed
520  if(freeFaces.find(face) != freeFaces.end())
521  freeNode = i;
522  }
523  if(priorityNode > -1)
524  {
525  fixNode(el, priorityNode);
526  }
527  else if(freeNode > -1)
528  {
529  nodePriority[el[freeNode]] = currNodePriority;
530  fixNode(el, freeNode);
531  --currNodePriority;
532  }
533  else //fix a random node
534  {
535  nodePriority[el[type1node_]] = currNodePriority;
536  fixNode(el, type1node_);
537  --currNodePriority;
538  }
539 
540  FaceElementType faceElement;
541  //walk over all faces
542  //add face 2 to freeFaces - if already free -great, match and keep, if active and not free, match and erase
543  getFace(el,type1face_, faceElement);
544  getFace(el,type1face_, face);
545 
546  const auto freeFacesEnd = freeFaces.end();
547  const auto freeFaceIt = freeFaces.find(face);
548  if( freeFaceIt != freeFacesEnd)
549  {
550  while(!checkFaceCompatibility(faceElement))
551  {
552  rotate(el);
553  }
554  freeFaceIt->second[1] = elIndex;
555  }
556  else if(activeFaces.find(face) != activeFaces.end())
557  {
558  while(!checkFaceCompatibility(faceElement))
559  {
560  rotate(el);
561  }
562  activeFaces.erase(face);
563  }
564  else
565  {
566  freeFaces.insert({{face,{elIndex,elIndex}}});
567  }
568  //add others to activeFaces - if already there, delete, if already free, match and erase
569  //for(int i=0; i<4; ++i ) //auto&& i : {0,1,2,3})
570  for(auto&& i : {0,1,2,3})
571  {
572  if (i == type1face_) continue;
573 
574  getFace(el,i,face);
575  getFace(el,i,faceElement);
576 
577  unsigned int neighborIndex = faceElement.second[0] == elIndex ? faceElement.second[1] : faceElement.second[0];
578  if(freeFaces.find(face) != freeFacesEnd)
579  {
580  while(!checkFaceCompatibility(faceElement))
581  {
582  rotate(elements_[neighborIndex]);
583  }
584  }
585  else if(activeFaces.find(face) != activeFaces.end())
586  {
587  if(!checkFaceCompatibility(faceElement))
588  {
589  checkFaceCompatibility(faceElement,true) ;
590  return false;
591  }
592  activeFaces.erase(face);
593  }
594  else
595  {
596  activeFaces.insert({{face,{elIndex,elIndex}}});
597  }
598  }
599  }
600 
601  //now postprocessing of freeFaces. possibly - not really necessary, has to be thought about
602  //useful for parallelization .
603  /*
604  for(auto&& face : freeFaces)
605  {
606  unsigned int elementIndex = face.second[0];
607  unsigned int neighborIndex = face.second[1];
608  //give refinement edge positive priority
609  //and non-refinement edge negative priority less than -1 and counting
610  }
611  */
612  return true;
613  }
614 
615  void returnElements(std::vector<ElementType> & elements,
616  std::vector<bool>& elementOrientation,
617  std::vector<int>& types,
618  const bool stevenson = false )
619  {
620  if( stevenson )
621  {
623  }
624  else
625  {
627  }
628 
629  //This needs to happen, before the boundaryIds are
630  //recreated in the GridFactory
631  const int nElements = elements_.size();
632  for( int el = 0; el<nElements; ++el )
633  {
634  // in ALU only elements with negative orientation can be inserted
635  if( ! elementOrientation_[ el ] )
636  {
637  // the refinement edge is 0 -- 1, so we can swap 2 and 3
638  std::swap( elements_[ el ][ 2 ], elements_[ el ][ 3 ] );
639  }
640  }
641 
642  elements = elements_;
643  elementOrientation = elementOrientation_;
644  types = types_;
645  }
646 
648  {
650  {
651  swapRefinementType();
652  }
653  }
654 
656  {
657  if( ! stevensonRefinement_ )
658  {
659  swapRefinementType();
660  }
661  }
662 
663 private:
664  void swapRefinementType()
665  {
666  const int nElements = elements_.size();
667  for( int el = 0; el<nElements; ++el )
668  {
670  std::swap(elements_[ el ][ 1 ], elements_[ el ][ 3 ]);
671  }
672 
673  // swap refinement flags
677  type1node_ = stevensonRefinement_ ? 1 : 2 ;
678  type1face_ = ( 3 - type1node_ );
679  }
680 
681  //switch vertices 2,3 for all elements with elemIndex % 2
682  void applyStandardOrientation ()
683  {
684  int i = 0;
685  for(auto & element : elements_ )
686  {
687  if ( i % 2 == 0 )
688  {
689  std::swap(element[2],element[3]);
691  }
692  ++i;
693  }
694  types_.resize(elements_.size(), 0);
695  }
696 
697  //check face for compatibility
698  bool checkFaceCompatibility(std::pair<FaceType, EdgeType> face, bool verbose = false)
699  {
700  EdgeType edge1,edge2;
701  int elIndex = face.second[0];
702  int neighIndex = face.second[1];
703  if(elIndex != neighIndex)
704  {
705  getRefinementEdge(elements_[elIndex], face.first, edge1, types_[elIndex]);
706  getRefinementEdge(elements_[neighIndex], face.first, edge2, types_[neighIndex]);
707  if(edge1 != edge2)
708  {
709  if (verbose)
710  {
711  /* std::cerr << "Face: " << face.first[0] << ", " << face.first[1] << ", " << face.first[2]
712  << " has refinement edge: " << edge1[0] << ", " << edge1[1] <<
713  " from one side and: " << edge2[0] << ", " << edge2[1] <<
714  " from the other." << std::endl; */
715  printElement(elIndex);
716  printElement(neighIndex);
717  }
718  return false;
719  }
720  }
721  return true;
722  }
723 
724  //check face for strong compatibility
725  bool checkStrongCompatibility(std::pair<FaceType, EdgeType> face, bool verbose = false)
726  {
727  int elIndex = face.second[0];
728  int neighIndex = face.second[1];
729  bool result = false;
730  //if we are not boundary
731  if(elIndex != neighIndex)
732  {
733  int elType = types_[elIndex];
734  int neighType = types_[neighIndex];
735  unsigned int elFaceIndex = getFaceIndex(elements_[elIndex], face.first);
736  unsigned int elNodeIndex = 3 - elFaceIndex;
737  unsigned int neighFaceIndex = getFaceIndex(elements_[neighIndex], face.first);
738  unsigned int neighNodeIndex = 3 - neighFaceIndex;
739  if(elType == neighType)
740  {
741  //if the local face indices coincide, they are reflected neighbors
742  // if the refinement edge is not contained in the shared face, we
743  // have reflected neighbors of the children, if the face is in the same direction
744  // and the other edge of the refinement edge is the missing one.
745  if( elNodeIndex == neighNodeIndex ||
746  (elNodeIndex == type0nodes_[0] && neighNodeIndex == type0nodes_[1]) ||
747  (elNodeIndex == type0nodes_[1] && neighNodeIndex == type0nodes_[0]) )
748  { result = true; }
749  }
750  else
751  {
752  if( elType % 3 == (neighType +1) %3 )
753  {
754  std::swap(elType, neighType);
755  std::swap(elNodeIndex, neighNodeIndex);
756  std::swap(elFaceIndex, neighFaceIndex);
757  std::swap(elIndex, neighIndex);
758  }
759  switch (elType) {
760  case 0:
761  assert(neighType == 1);
762  if( ( neighNodeIndex == type1node_ ) && ( elNodeIndex == type0nodes_[0] || elNodeIndex == type0nodes_[1] ) )
763  { result = true; }
764  break;
765  case 1:
766  assert(neighType == 2);
767  if( ( neighNodeIndex == type1node_ ) && ( elNodeIndex == type0nodes_[0] || elNodeIndex == type0nodes_[1] ) )
768  { result = true; }
769  break;
770  case 2:
771  assert(neighType == 0);
772  if( ( neighNodeIndex == type1node_ ) && ( elNodeIndex == type0nodes_[0] || elNodeIndex == type0nodes_[1] ) )
773  { result = true; }
774  break;
775  default:
776  std::cerr << "No other types than 0,1,2 implemented. Aborting" << std::endl;
777  abort();
778  }
779  }
780  }
781  else //boundary faces
782  {
783  result =true;
784  }
785  if (verbose && result == false )
786  {
787  /* std::cerr << "Face: " << face.first[0] << ", " << face.first[1] << ", " << face.first[2]
788  << " has refinement edge: " << edge1[0] << ", " << edge1[1] <<
789  " from one side and: " << edge2[0] << ", " << edge2[1] <<
790  " from the other." << std::endl; */
791  printElement(elIndex);
792  printElement(neighIndex);
793  }
794  return result;
795  }
796 
797  void fixNode(ElementType& el, int node)
798  {
799  if(!(node == type1node_))
800  {
801  //swap the node at the right position
802  std::swap(el[node],el[type1node_]);
803  //also swap two other nodes to keep the volume positive
804  //2 and 0 are never type1node_
805  std::swap(el[(type1node_+1)%4],el[(type1node_+2)%4]);
806  }
807  }
808 
809  //The rotations that keep the type 1 node fixed
810  void rotate(ElementType& el)
811  {
812  std::swap(el[(type1node_+1)%4],el[(type1node_+2)%4]);
813  std::swap(el[(type1node_+1)%4],el[(type1node_+3)%4]);
814  }
815 
816  //get the sorted face of an element to a corresponding index
817  //the index coincides with the missing vertex
818  void getFace(ElementType el, int faceIndex, FaceType & face )
819  {
820  switch(faceIndex)
821  {
822  case 3 :
823  face = {el[1],el[2],el[3]};
824  break;
825  case 2 :
826  face = {el[0],el[2],el[3]};
827  break;
828  case 1 :
829  face = {el[0],el[1],el[3]};
830  break;
831  case 0 :
832  face = {el[0],el[1],el[2]};
833  break;
834  default :
835  std::cerr << "index " << faceIndex << " NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
836  std::abort();
837  }
838  std::sort(face.begin(), face.end());
839  return;
840  }
841 
842  void getFace(ElementType el, int faceIndex, FaceElementType & faceElement)
843  {
844  FaceType face;
845  getFace(el, faceIndex, face);
846  faceElement = *neighbours_.find(face);
847  }
848 
849  //get the sorted initial refinement edge of a face of an element
850  //this has to be adjusted, when using stevensonRefinement
851  //orientation switches indices 2 and 3 in the internal ordering
852  void getRefinementEdge(ElementType el, int faceIndex, EdgeType & edge, int type )
853  {
854  if(type == 0)
855  {
857  {
858  switch(faceIndex)
859  {
860  case 0 :
861  edge = {el[0],el[2]};
862  break;
863  case 1 :
864  edge = {el[0],el[3]};
865  break;
866  case 2 :
867  edge = {el[0],el[3]};
868  break;
869  case 3 :
870  edge = {el[1],el[3]};
871  break;
872  default :
873  std::cerr << "index " << faceIndex << " NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
874  std::abort();
875  }
876  }
877  else //ALBERTA Refinement
878  {
879  switch(faceIndex)
880  {
881  case 0 :
882  edge = {el[0],el[1]};
883  break;
884  case 1 :
885  edge = {el[0],el[1]};
886  break;
887  case 2 :
888  edge = {el[0],el[2]};
889  break;
890  case 3 :
891  edge = {el[1],el[3]};
892  break;
893  default :
894  std::cerr << "index " << faceIndex << " NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
895  std::abort();
896  }
897  }
898  }
899  else if(type == 1 || type == 2)
900  {
902  {
903  switch(faceIndex)
904  {
905  case 0 :
906  edge = {el[0],el[2]};
907  break;
908  case 1 :
909  edge = {el[0],el[3]};
910  break;
911  case 2 :
912  edge = {el[0],el[3]};
913  break;
914  case 3 :
915  edge = {el[2],el[3]};
916  break;
917  default :
918  std::cerr << "index " << faceIndex << " NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
919  std::abort();
920  }
921  }
922  else //ALBERTA Refinement
923  {
924  switch(faceIndex)
925  {
926  case 0 :
927  edge = {el[0],el[1]};
928  break;
929  case 1 :
930  edge = {el[0],el[1]};
931  break;
932  case 2 :
933  edge = {el[0],el[3]};
934  break;
935  case 3 :
936  edge = {el[1],el[3]};
937  break;
938  default :
939  std::cerr << "index " << faceIndex << " NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
940  std::abort();
941  }
942  }
943  }
944  else
945  {
946  std::cerr << "no other types than 0, 1, 2 implemented." << std::endl;
947  std::abort();
948  }
949  std::sort(edge.begin(),edge.end());
950  return;
951  }
952 
953  //get the sorted initial refinement edge of a face of an element
954  void getRefinementEdge(ElementType el, FaceType face, EdgeType & edge, int type )
955  {
956  getRefinementEdge(el, getFaceIndex(el, face), edge, type);
957  return;
958  }
959 
960  //get the index of a face in an elements
961  //this could be improved by exploiting that faces are sorted
962  int getFaceIndex(ElementType el, FaceType face)
963  {
964  for(int i =0; i<4 ; ++i)
965  {
966  if(!( el[i] == face[0] || el[i] == face[1] || el[i] == face[2] ) )
967  return 3 - i ;
968  }
969  return -1;
970  }
971 
972  //build the structure containing the neighbors
973  //consists of a face and the two indices belonging to
974  //the elements that share said face
975  //boundary faces have two times the same index
976  //this is executed in the constructor
977  void buildNeighbors()
978  {
979  // clear existing structures
980  neighbours_.clear();
981 
982  FaceType face;
983  EdgeType indexPair;
984 
985  unsigned int index = 0;
986  const auto nend = neighbours_.end();
987  for(auto&& el : elements_)
988  {
989  for(int i = 0; i< 4; ++i)
990  {
991  getFace(el, i, face);
992  auto faceInList = neighbours_.find(face);
993  if(faceInList == nend)
994  {
995  indexPair = {index, index};
996  neighbours_.insert(std::make_pair (face, indexPair ) );
997  }
998  else
999  {
1000  assert(faceInList != neighbours_.end());
1001  faceInList->second[1] = index;
1002  }
1003  }
1004  ++index;
1005  }
1006  }
1007 
1016  void calculateV0(int variante, int threshold)
1017  {
1018  //For now, everything is in V0
1019  switch (variante) {
1020  case 1:
1021  {
1022  std::fill(containedInV0_.begin(),containedInV0_.end(), false);
1023  std::vector<int> numberOfAdjacentRefEdges(nVertices_, 0);
1024  //we assume that the edges have been sorted and
1025  //the refinement edge is, where it belongs
1026  for(auto&& el : elements_)
1027  {
1028  numberOfAdjacentRefEdges [ el [ type0nodes_[ 0 ] ] ] ++;
1029  numberOfAdjacentRefEdges [ el [ type0nodes_[ 1 ] ] ] ++;
1030  }
1031  for(unsigned int i = 0; i <nVertices_ ; ++i)
1032  {
1033  if(numberOfAdjacentRefEdges[ i ] >= threshold )
1034  {
1035  containedInV0_[ i ] = true;
1036  }
1037  }
1038  }
1039  break;
1040  case 2:
1041  {
1042  std::vector<int> numberOfAdjacentElements(nVertices_, 0);
1043  std::vector<bool> vertexOnBoundary(nVertices_, false);
1044  for(auto&& neigh : neighbours_)
1045  {
1046  //We want to treat boundary vertices differently
1047  if(neigh.second[1] == neigh.second[0])
1048  {
1049  for(unsigned int i =0; i < 3 ; ++i)
1050  {
1051  vertexOnBoundary[ neigh.first [ i ] ] =true;
1052  }
1053  }
1054  }
1055  //calculate numberOfAdjacentElements
1056  for(auto&& el : elements_)
1057  {
1058  for(unsigned int i =0; i < 4; ++i)
1059  {
1060  numberOfAdjacentElements[ el [ i ] ] ++;
1061  }
1062  }
1063  for(unsigned int i = 0; i <nVertices_ ; ++i)
1064  {
1065  double bound = vertexOnBoundary[ i ] ? threshold / 2. : threshold ;
1066  if(numberOfAdjacentElements[ i ] < bound )
1067  {
1068  containedInV0_[ i ] = false;
1069  }
1070  }
1071 
1072  }
1073  break;
1074  case 3:
1075  {
1076  std::default_random_engine generator;
1077  std::uniform_int_distribution<int> distribution(1,6);
1078  for(unsigned int i = 0; i < nVertices_; ++i)
1079  {
1080  int roll = distribution(generator); // generates number in the range 1..6
1081  if(roll < 3)
1082  {
1083  containedInV0_[ i ] = false;
1084  }
1085  }
1086  }
1087  break;
1088  default: ;
1089  }
1090  int sizeOfV1 = 0;
1091  int sizeOfV0 = 0;
1092  for(unsigned int i =0 ; i < nVertices_; ++i)
1093  {
1094  if( containedInV0_ [ i ] )
1095  {
1096  ++sizeOfV0;
1097  }
1098  else
1099  {
1100  ++sizeOfV1;
1101  }
1102  }
1103 #ifndef NDEBUG
1104  std::cout << "#V0 #V1" << std::endl << sizeOfV0 << " " << sizeOfV1 << std::endl;
1105 #endif
1106  }
1107 
1108 }; //class bisectioncompatibility
1109 
1110 
1111 
1112 #endif //ALUGRID_COMPATIBILITY_CHECK_HH_INCLUDED
Definition: bisectioncompatibility.hh:16
static int & variant()
Definition: bisectioncompatibility.hh:17
static int & threshold()
Definition: bisectioncompatibility.hh:23
static int & useAnnouncedEdge()
Definition: bisectioncompatibility.hh:29
Definition: bisectioncompatibility.hh:41
bool type1Algorithm()
Definition: bisectioncompatibility.hh:483
void printElement(int index)
Definition: bisectioncompatibility.hh:171
BisectionCompatibility(const VertexVectorType &vertices, const std::vector< ElementType > &elements, const bool stevenson)
Definition: bisectioncompatibility.hh:92
std::vector< bool > elementOrientation_
Definition: bisectioncompatibility.hh:58
void printNB()
Definition: bisectioncompatibility.hh:162
std::vector< int > types_
Definition: bisectioncompatibility.hh:67
unsigned int type1face_
Definition: bisectioncompatibility.hh:79
FaceMapType neighbours_
Definition: bisectioncompatibility.hh:60
std::array< unsigned int, 3 > FaceType
Definition: bisectioncompatibility.hh:47
const int variant_
Definition: bisectioncompatibility.hh:85
int stronglyCompatibleFaces()
Definition: bisectioncompatibility.hh:114
const VertexVectorType & vertices_
Definition: bisectioncompatibility.hh:54
const size_t nVertices_
Definition: bisectioncompatibility.hh:62
const int threshold_
Definition: bisectioncompatibility.hh:86
unsigned int type1node_
Definition: bisectioncompatibility.hh:77
VertexVector VertexVectorType
Definition: bisectioncompatibility.hh:45
std::vector< unsigned int > ElementType
Definition: bisectioncompatibility.hh:48
const bool useAnnouncedEdge_
Definition: bisectioncompatibility.hh:87
EdgeType type0nodes_
Definition: bisectioncompatibility.hh:73
std::array< unsigned int, 2 > EdgeType
Definition: bisectioncompatibility.hh:49
void alberta2Stevenson()
Definition: bisectioncompatibility.hh:655
std::pair< FaceType, EdgeType > FaceElementType
Definition: bisectioncompatibility.hh:51
void returnElements(std::vector< ElementType > &elements, std::vector< bool > &elementOrientation, std::vector< int > &types, const bool stevenson=false)
Definition: bisectioncompatibility.hh:615
std::vector< ElementType > elements_
Definition: bisectioncompatibility.hh:57
void stevenson2Alberta()
Definition: bisectioncompatibility.hh:647
std::map< FaceType, EdgeType > FaceMapType
Definition: bisectioncompatibility.hh:50
bool make6CompatibilityCheck()
Definition: bisectioncompatibility.hh:151
EdgeType type0faces_
Definition: bisectioncompatibility.hh:75
std::vector< bool > containedInV0_
Definition: bisectioncompatibility.hh:65
bool stevensonRefinement_
Definition: bisectioncompatibility.hh:70
bool type0Algorithm()
Definition: bisectioncompatibility.hh:195
bool compatibilityCheck()
Definition: bisectioncompatibility.hh:142