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