dune-pdelab  2.4.1
hangingnodemanager.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 
4 #ifndef HANGINGNODEMANAGER_HH
5 #define HANGINGNODEMANAGER_HH
6 
7 #include<dune/grid/common/grid.hh>
8 #include<dune/grid/common/mcmgmapper.hh>
9 #include<dune/common/float_cmp.hh>
10 
11 #include"../common/geometrywrapper.hh"
12 
13 namespace Dune {
14  namespace PDELab {
15 
16  template<class Grid, class BoundaryFunction>
18  {
19  private:
20 #ifdef DEBUG
21  enum{ verbosity = 2 };
22 #else
23  enum{ verbosity = 0 };
24 #endif
25  typedef typename Grid::LeafIndexSet::IndexType IndexType;
26 
27  private:
28  class NodeInfo
29  {
30  public:
31  // Minimum level of elements containing this node
32  unsigned short minimum_level;
33 
34  // Maximum level of elements containing this node
35  unsigned short maximum_level;
36 
37  // Minimum level of elements touching this node
38  unsigned short minimum_touching_level;
39 
40  bool is_boundary;
41 
42  void addLevel(unsigned short level)
43  {
44  minimum_level = std::min(minimum_level,level);
45  maximum_level = std::max(maximum_level,level);
46  }
47 
48  void addTouchingLevel(unsigned short level)
49  {
50  minimum_touching_level = std::min(minimum_touching_level,level);
51  }
52 
53  NodeInfo() : minimum_level(std::numeric_limits<unsigned short>::max()),
54  maximum_level(0),
55  minimum_touching_level(std::numeric_limits<unsigned short>::max()),
56  is_boundary(false)
57  {}
58  };
59 
60  std::vector<NodeInfo> node_info;
61 
62  public:
63 
64  class NodeState
65  {
66  friend class HangingNodeManager;
67  private:
68  bool is_boundary;
69  bool is_hanging;
70 
71  public:
72 
73  inline bool isBoundary() const { return is_boundary; }
74  inline bool isHanging() const { return is_hanging; }
75 
76  NodeState() :is_boundary(false),
77  is_hanging(false)
78  {}
79  };
80 
81 
82  typedef typename Grid::LeafGridView GridView;
83  enum {dim = GridView::dimension};
84  typedef typename GridView::template Codim<0>::EntityPointer CellEntityPointer;
85  typedef typename GridView::template Codim<0>::Entity Cell;
86 
87  typedef typename GridView::template Codim<dim>::EntityPointer VertexEntityPointer;
88  typedef typename GridView::template Codim<0>::Iterator Iterator;
89  typedef typename GridView::IntersectionIterator IntersectionIterator;
90  typedef typename GridView::Grid::ctype ctype;
91  typedef typename Dune::FieldVector<ctype,dim> Point;
92  typedef typename Dune::FieldVector<ctype,dim-1> FacePoint;
93 
94  typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,
95  MCMGElementLayout> CellMapper;
96 
97  Grid & grid;
98  const BoundaryFunction & boundaryFunction;
99  CellMapper cell_mapper;
100 
101  public:
102 
103  void analyzeView()
104  {
105  cell_mapper.update();
106  const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
107 
108  node_info = std::vector<NodeInfo>(indexSet.size(dim));
109 
110  GridView gv = grid.leafGridView();
111 
112  // loop over all codim<0> leaf elements of the partially refined grid
113  for(const auto& cell : elements(gv)) {
114 
115  const Dune::ReferenceElement<double,dim> &
116  reference_element =
117  Dune::ReferenceElements<double,dim>::general(cell.geometry().type());
118 
119 
120  // level of this element
121  const unsigned short level = cell.level();
122 
123  // number of vertices in this element
124  const IndexType v_size = reference_element.size(dim);
125 
126  // update minimum_level and maximum_level for vertices in this
127  // cell
128  // loop over all vertices of the element
129  for(IndexType i=0; i<v_size; ++i){
130  const IndexType v_globalindex = indexSet.subIndex(cell,i,dim);
131  NodeInfo& ni = node_info[v_globalindex];
132  ni.addLevel(level);
133 
134  if(verbosity>10){
135  // This will produce a lot of output on the screen!
136  std::cout << " cell-id=" << cell_mapper.index(cell);
137  std::cout << " level=" << level;
138  std::cout << " v_size=" << v_size;
139  std::cout << " v_globalindex = " << v_globalindex;
140  std::cout << " maximum_level = " << ni.maximum_level;
141  std::cout << " minimum_level = " << ni.minimum_level;
142  std::cout << std::endl;
143  }
144 
145  }
146 
147  // Now we still have to update minimum_touching_level for this
148  // cell
149 
150  typedef typename IntersectionIterator::Intersection Intersection;
151 
152  // Loop over faces
153  unsigned int intersection_index = 0;
154  for(const auto& intersection : intersections(gv,cell)) {
155  ++intersection_index;
156 
157  const Dune::ReferenceElement<double,dim-1> &
158  reference_face_element =
159  Dune::ReferenceElements<double,dim-1>::general(intersection.geometry().type());
160 
161  const int eLocalIndex = intersection.indexInInside();
162  const int e_level = intersection.inside().level();
163 
164  // number of vertices on the face
165  const int e_v_size = reference_element.size(eLocalIndex,1,dim);
166 
167  if(intersection.boundary()) {
168 
169  // loop over vertices on the face
170  for(int i=0; i<e_v_size;++i){
171  const int e_v_index = reference_element.subEntity(eLocalIndex,1,i,dim);
172  const IndexType v_globalindex = indexSet.subIndex(cell,e_v_index,dim);
173 
174  const FacePoint facelocal_position = reference_face_element.position(i,dim-1);
175 
176  /*
177  typename BoundaryFunction::Traits::RangeType boundary_value;
178  boundaryFunction.evaluate(IntersectionGeometry<Intersection>(*fit,intersection_index),
179  facelocal_position,boundary_value);
180  if(boundary_value)
181  node_info[v_globalindex].is_boundary=true;
182  else
183  node_info[v_globalindex].is_boundary=false;
184  */
185 
186  NodeInfo& ni = node_info[v_globalindex];
187  ni.is_boundary = boundaryFunction.isDirichlet(IntersectionGeometry<Intersection>(intersection,intersection_index),facelocal_position);
188  ni.addTouchingLevel(e_level);
189  }
190 
191  // We are done here - the remaining tests are only required for neighbor intersections
192  continue;
193  }
194 
195  const int f_level = intersection.outside().level();
196 
197  // a conforming face has no hanging nodes
198  if(intersection.conforming())
199  continue;
200 
201  // so far no support for initially non conforming grids
202  assert(e_level != f_level);
203 
204  // this check needs to be performed on the element containing
205  // the hanging node only
206  if(e_level < f_level)
207  continue;
208 
209  // loop over vertices on the face
210  for(int i=0; i<e_v_size;++i){
211  const int e_v_index = reference_element.subEntity(eLocalIndex,1,i,dim);
212  const IndexType v_globalindex = indexSet.subIndex(cell,e_v_index,dim);
213 
214  node_info[v_globalindex].addTouchingLevel(f_level);
215  }
216 
217  } // end of loop over faces
218 
219  } // end loop over codim<0> leaf elements
220  }
221 
222  HangingNodeManager(Grid & _grid, const BoundaryFunction & _boundaryFunction)
223  : grid(_grid),
224  boundaryFunction(_boundaryFunction),
225  cell_mapper(grid.leafGridView())
226  { analyzeView(); }
227 
228  const std::vector<NodeState> hangingNodes(const Cell& e) const
229  {
230  const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
231  std::vector<NodeState> is_hanging;
232 
233  const Dune::ReferenceElement<double,dim> &
234  reference_element =
235  Dune::ReferenceElements<double,dim>::general(e.geometry().type());
236 
237  // number of vertices in this element
238  const IndexType v_size = reference_element.size(dim);
239 
240  // make sure the return array is big enough
241  is_hanging.resize(v_size);
242 
243  // update minimum_level and maximum_level for vertices in this
244  // cell
245  // loop over vertices of the element
246  for(IndexType i=0; i<v_size; ++i){
247  const IndexType v_globalindex = indexSet.subIndex(e,i,dim);
248 
249  // here we make use of the fact that a node is hanging if and
250  // only if it touches a cell of a level smaller than the
251  // smallest level of all element containing the node
252  const NodeInfo & v_info = node_info[v_globalindex];
253  if(v_info.minimum_touching_level < v_info.minimum_level){
254  is_hanging[i].is_hanging = true;
255  is_hanging[i].is_boundary = v_info.is_boundary;
256 #ifndef NDEBUG
257  if(verbosity>1){
258  const Point & local = reference_element.position(i,dim);
259  const Point global = e.geometry().global(local);
260  if(verbosity)
261  std::cout << "Found hanging node with id " << v_globalindex << " at " << global << std::endl;
262  }
263 #endif
264  }
265  else{
266  is_hanging[i].is_hanging = false;
267  is_hanging[i].is_boundary = v_info.is_boundary;
268  }
269  }
270 
271  return is_hanging;
272  }
273 
275  {
276  if(verbosity)
277  std::cout << "Begin isolation of hanging nodes" << std::endl;
278 
279  const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
280 
281  size_t iterations(0);
282 
283  bool reiterate;
284 
285  // Iterate until the isolation condition is achieved.
286  do{
287  size_t refinements(0);
288  reiterate = false;
289 
290  GridView gv = grid.leafGridView();
291 
292  // loop over all codim<0> leaf elements of the partially refined grid
293  for(const auto& cell : elements(gv)) {
294 
295  const Dune::ReferenceElement<double,dim> &
296  reference_element =
297  Dune::ReferenceElements<double,dim>::general(cell.geometry().type());
298 
299  //std::cout << "cell center = " << it->geometry().center() << std::endl;
300 
301  // get the refinement level of the element
302  const unsigned short level = cell.level();
303 
304  //std::cout << "level = " << level << std::endl;
305 
306  // number of vertices in this element
307  const IndexType v_size = reference_element.size(dim);
308 
309  // update minimum_level and maximum_level for vertices in this
310  // cell
311  // loop over vertices of the element
312  for(IndexType i=0; i<v_size; ++i){
313 
314  const IndexType v_globalindex = indexSet.subIndex(cell,i,dim);
315 
316  const NodeInfo & v_info = node_info[v_globalindex];
317 
318  //std::cout << "maximum_level = " << v_info.maximum_level << std::endl;
319 
320  const unsigned short level_diff = v_info.maximum_level - level;
321  if( level_diff > 1){
322  grid.mark(1, cell); // Mark this element for an extra refinement if it has a hanging node belonging to a neighbouring element of a refinement level + 2 or more
323  reiterate = true; // Once an element has to be refined, the procedure needs to be repeated!
324  refinements++; // Count the number of refinements.
325 
326  if(verbosity>10){
327  // This will produce a lot of output on the screen!
328  std::cout << " cell-id=" << cell_mapper.index(cell);
329  std::cout << " level=" << level;
330  std::cout << " v_size=" << v_size;
331  std::cout << " v_globalindex = " << v_globalindex;
332  std::cout << std::endl;
333  std::cout << "Refining element nr " << cell_mapper.index(cell)
334  << " to isolate hanging nodes. Level diff = "
335  << v_info.maximum_level << " - " << level<< std::endl;
336  }
337  break;
338  }
339  } // end of loop over vertices
340 
341 
342  if( cell.geometry().type().isSimplex() ){
343  //
344  // SPECIAL CASE for SIMPLICES:
345  // Add extra check to find out "neighbouring" elements of level +2 or more
346  // which share only a hanging node, but do not share an intersection
347  // with the current element.
348  //
349  if( !reiterate ){
350 
351  //std::cout << "Extra check for SIMPLICES:" << std::endl;
352 
353  unsigned int intersection_index = 0;
354 
355  bool bJumpOut = false;
356 
357  // Loop over faces
358  for(const auto& intersection : intersections(gv,cell)) {
359  ++intersection_index;
360 
361  // only internal faces need to be taken care of
362  if(!intersection.boundary()) {
363 
364  const int e_level = intersection.inside().level();
365  const int f_level = intersection.outside().level();
366 
367  if( f_level > e_level ){
368 
369  // We have to locate the potential hanging node
370  // for which we do the extra Check.
371 
372  // get element-local index of the intersection
373  const int eLocalIndex = intersection.indexInInside();
374 
375  // Number of vertices on the face:
376  // A face(=edge) in a triangle has two vertices.
377  // A face(=triangle) in a tetrahedron has three vertices.
378  // const int e_v_size = reference_element.size(eLocalIndex,1,dim);
379 
380  int nEdgeCenters = 0;
381  if( dim == 2 ){
382  // 2D-case: We need to check later for each vertex of the
383  // neigbouring element if it lies on the center of the element edge.
384  // Take care: fit->geometry().center() might return the face
385  // center of a refined neighbouring element!
386  // But we need the face center of the coarse face of the
387  // current element. Therefore loop over vertices on the face
388  // to calculate the proper face center for the coarse face!
389  nEdgeCenters = 1;
390  }
391  else{
392  // 3D-case: We need to check later for each vertex of the
393  // neigbouring element if it lies on the center of one of
394  // the 3 edges of the element face.
395  nEdgeCenters = 3;
396  }
397  std::vector<Dune::FieldVector<ctype,dim> >
398  edgecenter( nEdgeCenters, Dune::FieldVector<ctype,dim>(0) );
399  //std::cout << " edgecenter = " << edgecenter << std::endl;
400 
401  // loop over center of the face (2d) or center of the edges of the face(3d)
402  for(int counter=0; counter<nEdgeCenters; ++counter){
403 
404  int cornerIndex1 = counter % dim;
405  int cornerIndex2 = (counter+1) % dim;
406 
407  const int e_v_index_1 =
408  reference_element.subEntity(eLocalIndex,1,cornerIndex1,dim);
409 
410  const int e_v_index_2 =
411  reference_element.subEntity(eLocalIndex,1,cornerIndex2,dim);
412 
413  auto vertex1 = cell.template subEntity<dim>(e_v_index_1);
414 
415  auto vertex2 = cell.template subEntity<dim>(e_v_index_2);
416 
417  edgecenter[counter] += vertex1.geometry().center();
418  edgecenter[counter] += vertex2.geometry().center();
419  edgecenter[counter] /= ctype( 2 );
420  //std::cout << " edgecenter = " << edgecenter << std::endl;
421 
422 
423  //
424  // check for the neighbouring element now...
425  //
426  const Dune::ReferenceElement<double,dim> &
427  nb_reference_element =
428  Dune::ReferenceElements<double,dim>::general( intersection.outside().geometry().type() );
429 
430  // number of vertices in that neigbouring element
431  const IndexType nb_v_size = nb_reference_element.size(dim);
432 
433  // loop over vertices of the neigbouring element
434  for(IndexType i=0; i<nb_v_size; ++i){
435 
436  auto nb_vertex = intersection.outside().template subEntity<dim>(i);
437 
438  bool doExtraCheck = false;
439 
440  Dune::FieldVector<ctype,dim> center_diff ( edgecenter[counter] );
441 
442  center_diff -= nb_vertex.geometry().center();
443 
444  //std::cout << "nb_vertex = " << nb_vertex->geometry().center() << std::endl;
445 
446  if( center_diff.two_norm2() < 5e-12 ){
447  doExtraCheck = true;
448  }
449 
450 
451  if( doExtraCheck ){
452 
453  //std::cout << "doExtraCheck for node at "
454  // << nb_vertex->geometry().center() << std::endl;
455 
456  const IndexType nb_v_globalindex = indexSet.index(nb_vertex);
457 
458  const NodeInfo & nb_v_info = node_info[nb_v_globalindex];
459 
460  const unsigned short level_diff = nb_v_info.maximum_level - level;
461 
462  if( level_diff > 1){
463  bJumpOut = true;
464  grid.mark(1, cell); // Mark this element for an extra refinement if it has a hanging node belonging to a neighbouring element of a refinement level + 2 or more
465  reiterate = true; // Once an element has to be refined, the procedure needs to be repeated!
466  refinements++; // Count the number of refinements.
467 
468  if(verbosity>10){
469  // This will produce a lot of output on the screen!
470  std::cout << " cell-id=" << cell_mapper.index(cell);
471  std::cout << " level=" << level;
472  std::cout << " v_size=" << v_size;
473  std::cout << std::endl;
474  std::cout << " Extra refining for element nr " << cell_mapper.index(cell)
475  << " to isolate hanging nodes. Level diff = "
476  << nb_v_info.maximum_level << " - " << level<< std::endl;
477  }
478  break;
479 
480  } // end if level_diff > 1
481 
482  } // end if( doExtraCheck )
483  if( bJumpOut ) break;
484  } // end of loop over vertices of the neigbouring element
485  if( bJumpOut ) break;
486  } // end counter loop
487 
488  } // end if( f_level > e_level )
489 
490  } // end if not boundary
491  if( bJumpOut ) break;
492  } // end of loop over faces of the element
493 
494  } // end if(!reiterate)
495 
496  } // end if geometry().type().isSimplex()
497 
498  } // end of loop over all codim<0> leaf elements
499 
500 
501  if(reiterate){
502  if(verbosity)
503  std::cout << "Re-adapt for isolation of hanging nodes..." << std::endl;
504  grid.preAdapt();
505  grid.adapt();
506  grid.postAdapt();
507  analyzeView();
508  }
509 
510  iterations++;
511  if(verbosity)
512  std::cout << "In iteration " << iterations << " there were "
513  << refinements << " grid cells to be refined additionally to isolate hanging nodes." << std::endl;
514  }while(reiterate);
515  }
516 
517  }; // end class HangingNodeManager
518 
519  } // end namespace PDELab
520 } // end namespace Dune
521 #endif
Definition: hangingnodemanager.hh:17
Grid::LeafGridView GridView
Definition: hangingnodemanager.hh:82
void analyzeView()
Definition: hangingnodemanager.hh:103
Dune::FieldVector< ctype, dim > Point
Definition: hangingnodemanager.hh:91
GridView::IntersectionIterator IntersectionIterator
Definition: hangingnodemanager.hh:89
const std::vector< NodeState > hangingNodes(const Cell &e) const
Definition: hangingnodemanager.hh:228
Grid & grid
Definition: hangingnodemanager.hh:97
Definition: hangingnodemanager.hh:64
GridView::template Codim< 0 >::EntityPointer CellEntityPointer
Definition: hangingnodemanager.hh:84
Definition: adaptivity.hh:27
GridView::template Codim< dim >::EntityPointer VertexEntityPointer
Definition: hangingnodemanager.hh:87
bool isBoundary() const
Definition: hangingnodemanager.hh:73
Dune::MultipleCodimMultipleGeomTypeMapper< GridView, MCMGElementLayout > CellMapper
Definition: hangingnodemanager.hh:95
HangingNodeManager(Grid &_grid, const BoundaryFunction &_boundaryFunction)
Definition: hangingnodemanager.hh:222
const Entity & e
Definition: localfunctionspace.hh:111
bool isHanging() const
Definition: hangingnodemanager.hh:74
NodeState()
Definition: hangingnodemanager.hh:76
GridView::template Codim< 0 >::Entity Cell
Definition: hangingnodemanager.hh:85
GridView::Grid::ctype ctype
Definition: hangingnodemanager.hh:90
const BoundaryFunction & boundaryFunction
Definition: hangingnodemanager.hh:98
Wrap intersection.
Definition: geometrywrapper.hh:56
GridView::template Codim< 0 >::Iterator Iterator
Definition: hangingnodemanager.hh:88
Dune::FieldVector< ctype, dim-1 > FacePoint
Definition: hangingnodemanager.hh:92
Definition: hangingnodemanager.hh:83
void adaptToIsolatedHangingNodes()
Definition: hangingnodemanager.hh:274
CellMapper cell_mapper
Definition: hangingnodemanager.hh:99