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