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