dune-pdelab  2.4-dev
hangingnode.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_HANGINGNODECONSTRAINTS_HH
3 #define DUNE_PDELAB_HANGINGNODECONSTRAINTS_HH
4 
5 #include <cstddef>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/geometry/referenceelements.hh>
9 #include<dune/geometry/type.hh>
11 #include"conforming.hh"
12 #include"hangingnodemanager.hh"
13 
14 namespace Dune {
15  namespace PDELab {
16 
20 
22  {
23  public:
25  {
26  public:
27  template<typename IG, typename LFS, typename T, typename FlagVector>
28  static void assembleConstraints(const FlagVector & nodeState_e, const FlagVector & nodeState_f,
29  const bool & e_has_hangingnodes, const bool & f_has_hangingnodes,
30  const LFS & lfs_e, const LFS & lfs_f,
31  T& trafo_e, T& trafo_f,
32  const IG& ig)
33  {
34  typedef IG Intersection;
35  typedef typename Intersection::EntityPointer CellEntityPointer;
36  typedef typename Intersection::Entity Cell;
37  typedef typename Intersection::Geometry FaceGeometry;
38  typedef typename FaceGeometry::ctype DT;
39  typedef typename LFS::Traits::SizeType SizeType;
40 
41  typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
42  const CellEntityPointer e = ig.inside();
43  const CellEntityPointer f = ! ig.boundary() ? ig.outside() : ig.inside();
44 
45  const std::size_t dimension = Intersection::dimension;
46 
47  typedef Dune::ReferenceElement<DT,dimension> GRE;
48  const GRE& refelement_e = Dune::ReferenceElements<DT,dimension>::general(e->type());
49  const GRE& refelement_f = Dune::ReferenceElements<DT,dimension>::general(f->type());
50 
51  // If both entities have hangingnodes, then the face is
52  // conforming and no constraints have to be applied.
53  if(e_has_hangingnodes && f_has_hangingnodes)
54  return;
55 
56  // Choose local function space etc for element with hanging nodes
57  const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
58  const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
59 
60  const Cell& cell = *(e_has_hangingnodes ? e : f);
61  const int faceindex = e_has_hangingnodes ? ig.indexInInside() : ig.indexInOutside();
62  const GRE & refelement = e_has_hangingnodes ? refelement_e : refelement_f;
63  const FlagVector & nodeState = e_has_hangingnodes ? nodeState_e : nodeState_f;
64  T & trafo = e_has_hangingnodes ? trafo_e : trafo_f;
65 
66  // A map mapping the local indices from the face to local
67  // indices of the cell
68  std::vector<int> m(refelement.size(faceindex,1,dimension));
69  for (int j=0; j<refelement.size(faceindex,1,dimension); j++)
70  m[j] = refelement.subEntity(faceindex,1,j,dimension);
71 
72  // A map mapping the local indices from the face to global gridview indices
73  std::vector<std::size_t> global_vertex_idx(refelement.size(faceindex,1,dimension));
74  for (int j=0; j<refelement.size(faceindex,1,dimension); ++j)
75  global_vertex_idx[j] = indexSet.subIndex(cell,refelement.subEntity(faceindex,1,j,dimension),dimension);
76 
77  // Create a DOFIndex that we will use to manually craft the correct dof indices for the constraints trafo
78  // We copy one of the indices from the LocalFunctionSpace; that way, we automatically get the correct
79  // TreeIndex into the DOFIndex and only have to fiddle with the EntityIndex.
80  typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
81 
82  typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
83  const GeometryType vertex_gt(0);
84 
85  // Find the corresponding entity in the reference element
86  for (int j=0; j<refelement.size(faceindex,1,dimension); j++){
87 
88  // The contribution factors of the base function bound to entity j
89  typename T::RowType contribution;
90 
91  if(dimension == 3){
92 
93  assert(nodeState.size() == 8);
94 
95  const SizeType i = 4*j;
96 
97  // Neigbor relations in local indices on a quadrilateral face:
98  // {Node, Direct Neighbor, Direct Neighbor, Diagonal Neighbor, Node, ...}
99  const unsigned int fi[16] = {0,1,2,3, 1,0,3,2, 2,0,3,1, 3,1,2,0};
100 
101  // Only hanging nodes have contribution to other nodes
102  if(nodeState[m[j]].isHanging()){
103 
104  // If both neighbors are hanging nodes, then this node
105  // is diagonal to the target of the contribution
106  if(nodeState[m[fi[i+1]]].isHanging() && nodeState[m[fi[i+2]]].isHanging())
107  {
108  GeometryIndexAccessor::store(dof_index.entityIndex(),
109  vertex_gt,
110  global_vertex_idx[fi[i+3]]);
111 
112  contribution[dof_index] = 0.25;
113 
114  GeometryIndexAccessor::store(dof_index.entityIndex(),
115  vertex_gt,
116  global_vertex_idx[j]);
117 
118  trafo[dof_index] = contribution;
119  }
120  // Direct neigbor
121  else if(!nodeState[m[fi[i+1]]].isHanging())
122  {
123  GeometryIndexAccessor::store(dof_index.entityIndex(),
124  vertex_gt,
125  global_vertex_idx[fi[i+1]]);
126 
127  contribution[dof_index] = 0.5;
128 
129  GeometryIndexAccessor::store(dof_index.entityIndex(),
130  vertex_gt,
131  global_vertex_idx[j]);
132 
133  trafo[dof_index] = contribution;
134  }
135  // Direct neigbor
136  else if(!nodeState[m[fi[i+2]]].isHanging())
137  {
138  GeometryIndexAccessor::store(dof_index.entityIndex(),
139  vertex_gt,
140  global_vertex_idx[fi[i+2]]);
141 
142  contribution[dof_index] = 0.5;
143 
144  GeometryIndexAccessor::store(dof_index.entityIndex(),
145  vertex_gt,
146  global_vertex_idx[j]);
147 
148  trafo[dof_index] = contribution;
149  }
150  }
151 
152  } else if(dimension == 2){
153 
154  assert(nodeState.size() == 4);
155 
156 
157  // Only hanging nodes have contribution to other nodes
158  if(nodeState[m[j]].isHanging()){
159 
160  const SizeType n_j = 1-j;
161 
162  assert( !nodeState[m[n_j]].isHanging() );
163 
164  GeometryIndexAccessor::store(dof_index.entityIndex(),
165  vertex_gt,
166  global_vertex_idx[n_j]);
167 
168  contribution[dof_index] = 0.5;
169 
170  GeometryIndexAccessor::store(dof_index.entityIndex(),
171  vertex_gt,
172  global_vertex_idx[j]);
173 
174  trafo[dof_index] = contribution;
175  }
176 
177  } // end if(dimension==3)
178 
179  } // j
180 
181  } // end of static void assembleConstraints
182 
183  }; // end of class CubeGridQ1Assembler
184 
185 
187  {
188  public:
189  template<typename IG,
190  typename LFS,
191  typename T,
192  typename FlagVector>
193  static void assembleConstraints( const FlagVector & nodeState_e,
194  const FlagVector & nodeState_f,
195  const bool & e_has_hangingnodes,
196  const bool & f_has_hangingnodes,
197  const LFS & lfs_e, const LFS & lfs_f,
198  T& trafo_e, T& trafo_f,
199  const IG& ig)
200  {
201  typedef IG Intersection;
202  typedef typename Intersection::EntityPointer CellEntityPointer;
203  typedef typename Intersection::Entity Cell;
204  typedef typename Intersection::Geometry FaceGeometry;
205  typedef typename FaceGeometry::ctype DT;
206  typedef typename LFS::Traits::SizeType SizeType;
207  typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
208 
209  const CellEntityPointer e = ig.inside();
210  const CellEntityPointer f = ! ig.boundary() ? ig.outside() : ig.inside();
211 
212  const std::size_t dimension = Intersection::dimension;
213 
214  typedef Dune::ReferenceElement<DT,dimension> GRE;
215  const GRE& refelement_e = Dune::ReferenceElements<DT,dimension>::general(e->type());
216  const GRE& refelement_f = Dune::ReferenceElements<DT,dimension>::general(f->type());
217 
218  // If both entities have hangingnodes, then the face is
219  // conforming and no constraints have to be applied.
220  if(e_has_hangingnodes && f_has_hangingnodes)
221  return;
222 
223  // Choose local function space etc for element with hanging nodes
224  const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
225  const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
226 
227  const Cell& cell = *(e_has_hangingnodes ? e : f);
228  const int faceindex = e_has_hangingnodes ? ig.indexInInside() : ig.indexInOutside();
229  const GRE & refelement = e_has_hangingnodes ? refelement_e : refelement_f;
230  const FlagVector & nodeState = e_has_hangingnodes ? nodeState_e : nodeState_f;
231  T & trafo = e_has_hangingnodes ? trafo_e : trafo_f;
232 
233  // A map mapping the local indices from the face to local
234  // indices of the cell
235  std::vector<int> m(refelement.size(faceindex,1,dimension));
236  for (int j=0; j<refelement.size(faceindex,1,dimension); j++)
237  m[j] = refelement.subEntity(faceindex,1,j,dimension);
238 
239  // A map mapping the local indices from the face to global gridview indices
240  std::vector<std::size_t> global_vertex_idx(refelement.size(faceindex,1,dimension));
241  for (int j=0; j<refelement.size(faceindex,1,dimension); ++j)
242  global_vertex_idx[j] = indexSet.subIndex(cell,refelement.subEntity(faceindex,1,j,dimension),dimension);
243 
244  // Create a DOFIndex that we will use to manually craft the correct dof indices for the constraints trafo
245  // We copy one of the indices from the LocalFunctionSpace; that way, we automatically get the correct
246  // TreeIndex into the DOFIndex and only have to fiddle with the EntityIndex.
247  typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
248 
249  typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
250  const GeometryType vertex_gt(0);
251 
252  // Find the corresponding entity in the reference element
253  for (int j=0; j<refelement.size(faceindex,1,dimension); j++){
254 
255  // The contribution factors of the base function bound to entity j
256  typename T::RowType contribution;
257 
258  if(dimension == 3){
259 
260  assert(nodeState.size() == 4);
261  // Only hanging nodes have contribution to other nodes
262  if(nodeState[m[j]].isHanging()){
263  for( int k=1; k<=2; ++k ){
264 
265  const SizeType n_j = (j+k)%3;
266 
267  if( !nodeState[m[n_j]].isHanging() )
268  {
269  GeometryIndexAccessor::store(dof_index.entityIndex(),
270  vertex_gt,
271  global_vertex_idx[n_j]);
272 
273  contribution[dof_index] = 0.5;
274 
275  GeometryIndexAccessor::store(dof_index.entityIndex(),
276  vertex_gt,
277  global_vertex_idx[j]);
278 
279  trafo[dof_index] = contribution;
280  }
281  }
282  }
283  } else if(dimension == 2){
284 
285  assert(nodeState.size() == 3);
286  // Only hanging nodes have contribution to other nodes
287  if(nodeState[m[j]].isHanging()){
288  const SizeType n_j = 1-j;
289  assert( !nodeState[m[n_j]].isHanging() );
290  // If both neighbors are hanging nodes, then this node
291  // is diagonal to the target of the contribution
292  GeometryIndexAccessor::store(dof_index.entityIndex(),
293  vertex_gt,
294  global_vertex_idx[n_j]);
295 
296  contribution[dof_index] = 0.5;
297 
298  GeometryIndexAccessor::store(dof_index.entityIndex(),
299  vertex_gt,
300  global_vertex_idx[j]);
301 
302  trafo[dof_index] = contribution;
303  }
304 
305 
306  } // end if(dimension==3)
307 
308  } // j
309 
310  } // end of static void assembleConstraints
311 
312  }; // end of class SimplexGridP1Assembler
313 
314  }; // end of class HangingNodesConstraintsAssemblers
315 
316 
318  // works in any dimension and on all element types
319  template <class Grid, class HangingNodesConstraintsAssemblerType, class BoundaryFunction>
321  {
322  private:
324  HangingNodeManager manager;
325 
326  public:
327  enum { doBoundary = true };
328  enum { doSkeleton = true };
329  enum { doVolume = false };
330  enum { dimension = Grid::dimension };
331 
333  bool adaptToIsolatedHangingNodes,
334  const BoundaryFunction & _boundaryFunction )
335  : manager(grid, _boundaryFunction)
336  {
337  if(adaptToIsolatedHangingNodes)
338  manager.adaptToIsolatedHangingNodes();
339  }
340 
341  void update( Grid & grid ){
342  manager.analyzeView();
343  manager.adaptToIsolatedHangingNodes();
344  }
345 
346 
348 
353  template<typename IG, typename LFS, typename T>
354  void skeleton (const IG& ig,
355  const LFS& lfs_e, const LFS& lfs_f,
356  T& trafo_e, T& trafo_f) const
357  {
358  typedef IG Intersection;
359  typedef typename Intersection::EntityPointer CellEntityPointer;
360  typedef typename Intersection::Geometry FaceGeometry;
361  typedef typename FaceGeometry::ctype DT;
362 
363  const CellEntityPointer e = ig.inside();
364  const CellEntityPointer f = ig.outside();
365 
366  const Dune::ReferenceElement<DT,dimension>& refelem_e
367  = Dune::ReferenceElements<DT,dimension>::general(e->type());
368  const Dune::ReferenceElement<DT,dimension>& refelem_f
369  = Dune::ReferenceElements<DT,dimension>::general(f->type());
370 
371  // the return values of the hanging node manager
372  typedef typename std::vector<typename HangingNodeManager::NodeState> FlagVector;
373  const FlagVector isHangingNode_e(manager.hangingNodes(*e));
374  const FlagVector isHangingNode_f(manager.hangingNodes(*f));
375 
376  // just to make sure that the hanging node manager is doing
377  // what is expected of him
378  assert(std::size_t(refelem_e.size(dimension))
379  == isHangingNode_e.size());
380  assert(std::size_t(refelem_f.size(dimension))
381  == isHangingNode_f.size());
382 
383  // the LOCAL indices of the faces in the reference element
384  const int faceindex_e = ig.indexInInside();
385  const int faceindex_f = ig.indexInOutside();
386 
387  bool e_has_hangingnodes = false;
388  {
389  for (int j=0; j<refelem_e.size(faceindex_e,1,dimension); j++){
390  const int index = refelem_e.subEntity(faceindex_e,1,j,dimension);
391  if(isHangingNode_e[index].isHanging())
392  {
393  e_has_hangingnodes = true;
394  break;
395  }
396  }
397  }
398  bool f_has_hangingnodes = false;
399  {
400  for (int j=0; j<refelem_f.size(faceindex_f,1,dimension); j++){
401  const int index = refelem_f.subEntity(faceindex_f,1,j,dimension);
402  if(isHangingNode_f[index].isHanging())
403  {
404  f_has_hangingnodes = true;
405  break;
406  }
407  }
408  }
409 
410  if(! e_has_hangingnodes && ! f_has_hangingnodes)
411  return;
412 
413  HangingNodesConstraintsAssemblerType::
414  assembleConstraints(isHangingNode_e, isHangingNode_f,
415  e_has_hangingnodes, f_has_hangingnodes,
416  lfs_e,lfs_f,
417  trafo_e, trafo_f,
418  ig);
419  } // skeleton
420 
421  }; // end of class HangingNodesDirichletConstraints
423 
424  }
425 }
426 
427 #endif
void update(Grid &grid)
Definition: hangingnode.hh:341
Hanging Node constraints construction.
Definition: hangingnode.hh:320
HangingNodesDirichletConstraints(Grid &grid, bool adaptToIsolatedHangingNodes, const BoundaryFunction &_boundaryFunction)
Definition: hangingnode.hh:332
void skeleton(const IG &ig, const LFS &lfs_e, const LFS &lfs_f, T &trafo_e, T &trafo_f) const
skeleton constraints
Definition: hangingnode.hh:354
const std::vector< NodeState > hangingNodes(const Cell &e) const
Definition: hangingnodemanager.hh:232
void analyzeView()
Definition: hangingnodemanager.hh:103
const E & e
Definition: interpolate.hh:172
static void assembleConstraints(const FlagVector &nodeState_e, const FlagVector &nodeState_f, const bool &e_has_hangingnodes, const bool &f_has_hangingnodes, const LFS &lfs_e, const LFS &lfs_f, T &trafo_e, T &trafo_f, const IG &ig)
Definition: hangingnode.hh:28
Dirichlet Constraints construction.
Definition: conforming.hh:36
const IG & ig
Definition: constraints.hh:147
Definition: adaptivity.hh:27
Definition: hangingnodemanager.hh:17
static void assembleConstraints(const FlagVector &nodeState_e, const FlagVector &nodeState_f, const bool &e_has_hangingnodes, const bool &f_has_hangingnodes, const LFS &lfs_e, const LFS &lfs_f, T &trafo_e, T &trafo_f, const IG &ig)
Definition: hangingnode.hh:193
void adaptToIsolatedHangingNodes()
Definition: hangingnodemanager.hh:278