4 #ifndef HANGINGNODEMANAGER_HH 5 #define HANGINGNODEMANAGER_HH 7 #include<dune/grid/common/grid.hh> 8 #include<dune/grid/common/mcmgmapper.hh> 9 #include<dune/common/float_cmp.hh> 11 #include"../common/geometrywrapper.hh" 16 template<
class Gr
id,
class BoundaryFunction>
21 enum{ verbosity = 2 };
23 enum{ verbosity = 0 };
25 typedef typename Grid::LeafIndexSet::IndexType IndexType;
32 unsigned short minimum_level;
35 unsigned short maximum_level;
38 unsigned short minimum_touching_level;
42 void addLevel(
unsigned short level)
44 minimum_level = std::min(minimum_level,level);
45 maximum_level = std::max(maximum_level,level);
48 void addTouchingLevel(
unsigned short level)
50 minimum_touching_level = std::min(minimum_touching_level,level);
53 NodeInfo() : minimum_level(std::numeric_limits<unsigned short>::max()),
55 minimum_touching_level(std::numeric_limits<unsigned short>::max()),
60 std::vector<NodeInfo> node_info;
74 inline bool isHanging()
const {
return is_hanging; }
83 enum {
dim = GridView::dimension};
85 typedef typename GridView::template Codim<0>::Entity
Cell;
88 typedef typename GridView::template Codim<0>::Iterator
Iterator;
90 typedef typename GridView::Grid::ctype
ctype;
91 typedef typename Dune::FieldVector<ctype,dim>
Point;
94 typedef Dune::MultipleCodimMultipleGeomTypeMapper<
GridView,
105 cell_mapper.update();
106 const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
108 node_info = std::vector<NodeInfo>(indexSet.size(
dim));
110 GridView gv = grid.leafGridView();
113 for(
const auto& cell : elements(gv)) {
115 const Dune::ReferenceElement<double,dim> &
117 Dune::ReferenceElements<double,dim>::general(cell.geometry().type());
121 const unsigned short level = cell.level();
124 const IndexType v_size = reference_element.size(
dim);
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];
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;
150 typedef typename IntersectionIterator::Intersection Intersection;
153 unsigned int intersection_index = 0;
154 for(
const auto& intersection : intersections(gv,cell)) {
155 ++intersection_index;
157 const Dune::ReferenceElement<double,
dim-1> &
158 reference_face_element =
159 Dune::ReferenceElements<double,dim-1>::general(intersection.geometry().type());
161 const int eLocalIndex = intersection.indexInInside();
162 const int e_level = intersection.inside().level();
165 const int e_v_size = reference_element.size(eLocalIndex,1,
dim);
167 if(intersection.boundary()) {
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);
174 const FacePoint facelocal_position = reference_face_element.position(i,
dim-1);
186 NodeInfo& ni = node_info[v_globalindex];
188 ni.addTouchingLevel(e_level);
195 const int f_level = intersection.outside().level();
198 if(intersection.conforming())
202 assert(e_level != f_level);
206 if(e_level < f_level)
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);
214 node_info[v_globalindex].addTouchingLevel(f_level);
224 boundaryFunction(_boundaryFunction),
225 cell_mapper(grid.leafGridView())
230 const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
231 std::vector<NodeState> is_hanging;
233 const Dune::ReferenceElement<double,dim> &
235 Dune::ReferenceElements<double,dim>::general(e.geometry().type());
238 const IndexType v_size = reference_element.size(
dim);
241 is_hanging.resize(v_size);
246 for(IndexType i=0; i<v_size; ++i){
247 const IndexType v_globalindex = indexSet.subIndex(e,i,
dim);
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;
258 const Point & local = reference_element.position(i,
dim);
259 const Point global = e.geometry().global(local);
261 std::cout <<
"Found hanging node with id " << v_globalindex <<
" at " << global << std::endl;
266 is_hanging[i].is_hanging =
false;
267 is_hanging[i].is_boundary = v_info.is_boundary;
277 std::cout <<
"Begin isolation of hanging nodes" << std::endl;
279 const typename GridView::IndexSet& indexSet = grid.leafGridView().indexSet();
281 size_t iterations(0);
287 size_t refinements(0);
290 GridView gv = grid.leafGridView();
293 for(
const auto& cell : elements(gv)) {
295 const Dune::ReferenceElement<double,dim> &
297 Dune::ReferenceElements<double,dim>::general(cell.geometry().type());
302 const unsigned short level = cell.level();
307 const IndexType v_size = reference_element.size(
dim);
312 for(IndexType i=0; i<v_size; ++i){
314 const IndexType v_globalindex = indexSet.subIndex(cell,i,
dim);
316 const NodeInfo & v_info = node_info[v_globalindex];
320 const unsigned short level_diff = v_info.maximum_level - level;
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;
342 if( cell.geometry().type().isSimplex() ){
353 unsigned int intersection_index = 0;
355 bool bJumpOut =
false;
358 for(
const auto& intersection : intersections(gv,cell)) {
359 ++intersection_index;
362 if(!intersection.boundary()) {
364 const int e_level = intersection.inside().level();
365 const int f_level = intersection.outside().level();
367 if( f_level > e_level ){
373 const int eLocalIndex = intersection.indexInInside();
380 int nEdgeCenters = 0;
397 std::vector<Dune::FieldVector<ctype,dim> >
398 edgecenter( nEdgeCenters, Dune::FieldVector<ctype,dim>(0) );
402 for(
int counter=0; counter<nEdgeCenters; ++counter){
404 int cornerIndex1 = counter %
dim;
405 int cornerIndex2 = (counter+1) %
dim;
407 const int e_v_index_1 =
408 reference_element.subEntity(eLocalIndex,1,cornerIndex1,
dim);
410 const int e_v_index_2 =
411 reference_element.subEntity(eLocalIndex,1,cornerIndex2,
dim);
413 auto vertex1 = cell.template subEntity<dim>(e_v_index_1);
415 auto vertex2 = cell.template subEntity<dim>(e_v_index_2);
417 edgecenter[counter] += vertex1.geometry().center();
418 edgecenter[counter] += vertex2.geometry().center();
419 edgecenter[counter] /=
ctype( 2 );
426 const Dune::ReferenceElement<double,dim> &
427 nb_reference_element =
428 Dune::ReferenceElements<double,dim>::general( intersection.outside().geometry().type() );
431 const IndexType nb_v_size = nb_reference_element.size(
dim);
434 for(IndexType i=0; i<nb_v_size; ++i){
436 auto nb_vertex = intersection.outside().template subEntity<dim>(i);
438 bool doExtraCheck =
false;
440 Dune::FieldVector<ctype,dim> center_diff ( edgecenter[counter] );
442 center_diff -= nb_vertex.geometry().center();
446 if( center_diff.two_norm2() < 5
e-12 ){
456 const IndexType nb_v_globalindex = indexSet.index(nb_vertex);
458 const NodeInfo & nb_v_info = node_info[nb_v_globalindex];
460 const unsigned short level_diff = nb_v_info.maximum_level - level;
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;
483 if( bJumpOut )
break;
485 if( bJumpOut )
break;
491 if( bJumpOut )
break;
503 std::cout <<
"Re-adapt for isolation of hanging nodes..." << std::endl;
512 std::cout <<
"In iteration " << iterations <<
" there were " 513 << refinements <<
" grid cells to be refined additionally to isolate hanging nodes." << std::endl;
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