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