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,
106 const typename GridView::IndexSet& indexSet =
grid.leafGridView().indexSet();
108 node_info = std::vector<NodeInfo>(indexSet.size(
dim));
112 Iterator it = gv.template begin<0>();
113 Iterator eit = gv.template end<0>();
118 const Dune::ReferenceElement<double,dim> &
120 Dune::ReferenceElements<double,dim>::general(it->geometry().type());
124 const unsigned short level = it->level();
127 const IndexType v_size = reference_element.size(
dim);
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];
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;
153 unsigned int intersection_index = 0;
156 typedef typename IntersectionIterator::Intersection Intersection;
159 for(;fit!=efit;++fit,++intersection_index){
161 const Dune::ReferenceElement<double,
dim-1> &
162 reference_face_element =
163 Dune::ReferenceElements<double,dim-1>::general(fit->geometry().type());
165 const int eLocalIndex = fit->indexInInside();
166 const int e_level = fit->inside()->level();
169 const int e_v_size = reference_element.size(eLocalIndex,1,
dim);
171 if((*fit).boundary()) {
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);
178 const FacePoint facelocal_position = reference_face_element.position(i,
dim-1);
190 NodeInfo& ni = node_info[v_globalindex];
192 ni.addTouchingLevel(e_level);
199 const int f_level = fit->outside()->level();
202 if(fit->conforming())
206 assert(e_level != f_level);
210 if(e_level < f_level)
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);
218 node_info[v_globalindex].addTouchingLevel(f_level);
234 const typename GridView::IndexSet& indexSet =
grid.leafGridView().indexSet();
235 std::vector<NodeState> is_hanging;
237 const Dune::ReferenceElement<double,dim> &
239 Dune::ReferenceElements<double,dim>::general(e.geometry().type());
242 const IndexType v_size = reference_element.size(
dim);
245 is_hanging.resize(v_size);
250 for(IndexType i=0; i<v_size; ++i){
251 const IndexType v_globalindex = indexSet.subIndex(e,i,
dim);
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;
262 const Point & local = reference_element.position(i,
dim);
263 const Point global = e.geometry().global(local);
265 std::cout <<
"Found hanging node with id " << v_globalindex <<
" at " << global << std::endl;
270 is_hanging[i].is_hanging =
false;
271 is_hanging[i].is_boundary = v_info.is_boundary;
281 std::cout <<
"Begin isolation of hanging nodes" << std::endl;
283 const typename GridView::IndexSet& indexSet =
grid.leafGridView().indexSet();
285 size_t iterations(0);
291 size_t refinements(0);
296 Iterator it = gv.template begin<0>();
297 Iterator eit = gv.template end<0>();
302 const Dune::ReferenceElement<double,dim> &
304 Dune::ReferenceElements<double,dim>::general(it->geometry().type());
309 const unsigned short level = it->level();
314 const IndexType v_size = reference_element.size(
dim);
319 for(IndexType i=0; i<v_size; ++i){
321 const IndexType v_globalindex = indexSet.subIndex(*it,i,
dim);
323 const NodeInfo & v_info = node_info[v_globalindex];
327 const unsigned short level_diff = v_info.maximum_level - level;
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;
349 if( it->geometry().type().isSimplex() ){
360 unsigned int intersection_index = 0;
364 bool bJumpOut =
false;
367 for(;fit!=efit;++fit,++intersection_index){
370 if(!(*fit).boundary()) {
372 const int e_level = fit->inside()->level();
373 const int f_level = fit->outside()->level();
375 if( f_level > e_level ){
381 const int eLocalIndex = fit->indexInInside();
388 int nEdgeCenters = 0;
405 std::vector<Dune::FieldVector<ctype,dim> >
406 edgecenter( nEdgeCenters, Dune::FieldVector<ctype,dim>(0) );
410 for(
int counter=0; counter<nEdgeCenters; ++counter){
412 int cornerIndex1 = counter %
dim;
413 int cornerIndex2 = (counter+1) %
dim;
415 const int e_v_index_1 =
416 reference_element.subEntity(eLocalIndex,1,cornerIndex1,
dim);
418 const int e_v_index_2 =
419 reference_element.subEntity(eLocalIndex,1,cornerIndex2,
dim);
422 it->template subEntity<dim>(e_v_index_1);
425 it->template subEntity<dim>(e_v_index_2);
427 edgecenter[counter] += vertex1->geometry().center();
428 edgecenter[counter] += vertex2->geometry().center();
429 edgecenter[counter] /=
ctype( 2 );
436 const Dune::ReferenceElement<double,dim> &
437 nb_reference_element =
438 Dune::ReferenceElements<double,dim>::general( fit->outside()->geometry().type() );
441 const IndexType nb_v_size = nb_reference_element.size(
dim);
444 for(IndexType i=0; i<nb_v_size; ++i){
447 fit->outside()->template subEntity<dim>(i);
449 bool doExtraCheck =
false;
451 Dune::FieldVector<ctype,dim> center_diff ( edgecenter[counter] );
453 center_diff -= nb_vertex->geometry().center();
457 if( center_diff.two_norm2() < 5
e-12 ){
467 const IndexType nb_v_globalindex = indexSet.index(*nb_vertex);
469 const NodeInfo & nb_v_info = node_info[nb_v_globalindex];
471 const unsigned short level_diff = nb_v_info.maximum_level - level;
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;
494 if( bJumpOut )
break;
496 if( bJumpOut )
break;
502 if( bJumpOut )
break;
514 std::cout <<
"Re-adapt for isolation of hanging nodes..." << std::endl;
523 std::cout <<
"In iteration " << iterations <<
" there were "
524 << refinements <<
" grid cells to be refined additionally to isolate hanging nodes." << std::endl;
GridView::template Codim< 0 >::Entity Cell
Definition: constraints/hangingnodemanager.hh:85
Dune::FieldVector< ctype, dim > Point
Definition: constraints/hangingnodemanager.hh:91
bool isBoundary() const
Definition: constraints/hangingnodemanager.hh:73
Grid & grid
Definition: constraints/hangingnodemanager.hh:97
NodeState()
Definition: constraints/hangingnodemanager.hh:76
void analyzeView()
Definition: constraints/hangingnodemanager.hh:103
Definition: constraints/hangingnodemanager.hh:17
Definition: constraints/hangingnodemanager.hh:64
Dune::MultipleCodimMultipleGeomTypeMapper< GridView, MCMGElementLayout > CellMapper
Definition: constraints/hangingnodemanager.hh:95
Definition: constraints/hangingnodemanager.hh:83
const std::vector< NodeState > hangingNodes(const Cell &e) const
Definition: constraints/hangingnodemanager.hh:232
const BoundaryFunction & boundaryFunction
Definition: constraints/hangingnodemanager.hh:98
GridView::template Codim< dim >::EntityPointer VertexEntityPointer
Definition: constraints/hangingnodemanager.hh:87
void adaptToIsolatedHangingNodes()
Definition: constraints/hangingnodemanager.hh:278
GridView::Grid::ctype ctype
Definition: constraints/hangingnodemanager.hh:90
CellMapper cell_mapper
Definition: constraints/hangingnodemanager.hh:99
GridView::IntersectionIterator IntersectionIterator
Definition: constraints/hangingnodemanager.hh:89
GridView::template Codim< 0 >::Iterator Iterator
Definition: constraints/hangingnodemanager.hh:88
Dune::FieldVector< ctype, dim-1 > FacePoint
Definition: constraints/hangingnodemanager.hh:92
bool isHanging() const
Definition: constraints/hangingnodemanager.hh:74
Grid::LeafGridView GridView
Definition: constraints/hangingnodemanager.hh:82
const E & e
Definition: interpolate.hh:172
HangingNodeManager(Grid &_grid, const BoundaryFunction &_boundaryFunction)
Definition: constraints/hangingnodemanager.hh:226
GridView::template Codim< 0 >::EntityPointer CellEntityPointer
Definition: constraints/hangingnodemanager.hh:84
Wrap intersection.
Definition: geometrywrapper.hh:56