DUNE PDELab (git)

hangingnodemanager.hh
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
10
11#include"../common/geometrywrapper.hh"
12
13namespace Dune {
14 namespace PDELab {
15
16 template<class Grid, class BoundaryFunction>
17 class HangingNodeManager
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
93
94 Grid & grid;
95 const BoundaryFunction & boundaryFunction;
96 CellMapper cell_mapper;
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
264 void adaptToIsolatedHangingNodes()
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 bool bJumpOut = false;
342
343 // Loop over faces
344 for(const auto& intersection : intersections(gv,cell)) {
345
346 // only internal faces need to be taken care of
347 if(!intersection.boundary()) {
348
349 const int e_level = intersection.inside().level();
350 const int f_level = intersection.outside().level();
351
352 if( f_level > e_level ){
353
354 // We have to locate the potential hanging node
355 // for which we do the extra Check.
356
357 // get element-local index of the intersection
358 const int eLocalIndex = intersection.indexInInside();
359
360 // Number of vertices on the face:
361 // A face(=edge) in a triangle has two vertices.
362 // A face(=triangle) in a tetrahedron has three vertices.
363 // const int e_v_size = reference_element.size(eLocalIndex,1,dim);
364
365 int nEdgeCenters = 0;
366 if( dim == 2 ){
367 // 2D-case: We need to check later for each vertex of the
368 // neigbouring element if it lies on the center of the element edge.
369 // Take care: fit->geometry().center() might return the face
370 // center of a refined neighbouring element!
371 // But we need the face center of the coarse face of the
372 // current element. Therefore loop over vertices on the face
373 // to calculate the proper face center for the coarse face!
374 nEdgeCenters = 1;
375 }
376 else{
377 // 3D-case: We need to check later for each vertex of the
378 // neigbouring element if it lies on the center of one of
379 // the 3 edges of the element face.
380 nEdgeCenters = 3;
381 }
382 std::vector<Dune::FieldVector<ctype,dim> >
383 edgecenter( nEdgeCenters, Dune::FieldVector<ctype,dim>(0) );
384 //std::cout << " edgecenter = " << edgecenter << std::endl;
385
386 // loop over center of the face (2d) or center of the edges of the face(3d)
387 for(int counter=0; counter<nEdgeCenters; ++counter){
388
389 int cornerIndex1 = counter % dim;
390 int cornerIndex2 = (counter+1) % dim;
391
392 const int e_v_index_1 =
393 reference_element.subEntity(eLocalIndex,1,cornerIndex1,dim);
394
395 const int e_v_index_2 =
396 reference_element.subEntity(eLocalIndex,1,cornerIndex2,dim);
397
398 auto vertex1 = cell.template subEntity<dim>(e_v_index_1);
399
400 auto vertex2 = cell.template subEntity<dim>(e_v_index_2);
401
402 edgecenter[counter] += vertex1.geometry().center();
403 edgecenter[counter] += vertex2.geometry().center();
404 edgecenter[counter] /= ctype( 2 );
405 //std::cout << " edgecenter = " << edgecenter << std::endl;
406
407
408 //
409 // check for the neighbouring element now...
410 //
411 auto nb_reference_element = referenceElement(intersection.outside().geometry());
412
413 // number of vertices in that neigbouring element
414 const IndexType nb_v_size = nb_reference_element.size(dim);
415
416 // loop over vertices of the neigbouring element
417 for(IndexType i=0; i<nb_v_size; ++i){
418
419 auto nb_vertex = intersection.outside().template subEntity<dim>(i);
420
421 bool doExtraCheck = false;
422
423 Dune::FieldVector<ctype,dim> center_diff ( edgecenter[counter] );
424
425 center_diff -= nb_vertex.geometry().center();
426
427 //std::cout << "nb_vertex = " << nb_vertex->geometry().center() << std::endl;
428
429 if( center_diff.two_norm2() < 5e-12 ){
430 doExtraCheck = true;
431 }
432
433
434 if( doExtraCheck ){
435
436 //std::cout << "doExtraCheck for node at "
437 // << nb_vertex->geometry().center() << std::endl;
438
439 const IndexType nb_v_globalindex = indexSet.index(nb_vertex);
440
441 const NodeInfo & nb_v_info = node_info[nb_v_globalindex];
442
443 const unsigned short level_diff = nb_v_info.maximum_level - level;
444
445 if( level_diff > 1){
446 bJumpOut = true;
447 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
448 reiterate = true; // Once an element has to be refined, the procedure needs to be repeated!
449 refinements++; // Count the number of refinements.
450
451 if(verbosity>10){
452 // This will produce a lot of output on the screen!
453 std::cout << " cell-id=" << cell_mapper.index(cell);
454 std::cout << " level=" << level;
455 std::cout << " v_size=" << v_size;
456 std::cout << std::endl;
457 std::cout << " Extra refining for element nr " << cell_mapper.index(cell)
458 << " to isolate hanging nodes. Level diff = "
459 << nb_v_info.maximum_level << " - " << level<< std::endl;
460 }
461 break;
462
463 } // end if level_diff > 1
464
465 } // end if( doExtraCheck )
466 if( bJumpOut ) break;
467 } // end of loop over vertices of the neigbouring element
468 if( bJumpOut ) break;
469 } // end counter loop
470
471 } // end if( f_level > e_level )
472
473 } // end if not boundary
474 if( bJumpOut ) break;
475 } // end of loop over faces of the element
476
477 } // end if(!reiterate)
478
479 } // end if geometry().type().isSimplex()
480
481 } // end of loop over all codim<0> leaf elements
482
483
484 if(reiterate){
485 if(verbosity)
486 std::cout << "Re-adapt for isolation of hanging nodes..." << std::endl;
487 grid.preAdapt();
488 grid.adapt();
489 grid.postAdapt();
490 analyzeView();
491 }
492
493 iterations++;
494 if(verbosity)
495 std::cout << "In iteration " << iterations << " there were "
496 << refinements << " grid cells to be refined additionally to isolate hanging nodes." << std::endl;
497 }while(reiterate);
498 }
499
500 }; // end class HangingNodeManager
501
502 } // end namespace PDELab
503} // end namespace Dune
504#endif // DUNE_PDELAB_CONSTRAINTS_HANGINGNODEMANAGER_HH
GridFamily::Traits::LeafGridView LeafGridView
type of view for leaf grid
Definition: grid.hh:400
Dune::Intersection< GridImp, IntersectionImp > Intersection
Type of Intersection this IntersectionIterator points to.
Definition: intersectioniterator.hh:111
void update(const GV &gridView)
Recalculates indices after grid adaptation.
Definition: mcmgmapper.hh:308
Different resources needed by all grid implementations.
Various ways to compare floating-point numbers.
Traits::IntersectionIterator IntersectionIterator
type of the intersection iterator
Definition: gridview.hh:92
Traits::IndexSet IndexSet
type of the index set
Definition: gridview.hh:86
Grid< dim, dimworld, ct, GridFamily >::LeafGridView leafGridView(const Grid< dim, dimworld, ct, GridFamily > &grid)
leaf grid view for the given grid
Definition: grid.hh:805
static constexpr int dimension
The dimension of the grid.
Definition: gridview.hh:134
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
Mapper for multiple codim and multiple geometry types.
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)