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