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