DUNE PDELab (git)

hangingnode.hh
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
8#include<dune/geometry/referenceelements.hh>
10#include<dune/pdelab/common/geometrywrapper.hh>
11#include"conforming.hh"
12#include"hangingnodemanager.hh"
13
14namespace Dune {
15 namespace PDELab {
16
20
21 class HangingNodesConstraintsAssemblers
22 {
23 public:
24 class CubeGridQ1Assembler
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(),
105 global_vertex_idx[fi[i+3]]);
106
107 contribution[dof_index] = 0.25;
108
109 GeometryIndexAccessor::store(dof_index.entityIndex(),
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(),
120 global_vertex_idx[fi[i+1]]);
121
122 contribution[dof_index] = 0.5;
123
124 GeometryIndexAccessor::store(dof_index.entityIndex(),
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(),
135 global_vertex_idx[fi[i+2]]);
136
137 contribution[dof_index] = 0.5;
138
139 GeometryIndexAccessor::store(dof_index.entityIndex(),
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(),
161 global_vertex_idx[n_j]);
162
163 contribution[dof_index] = 0.5;
164
165 GeometryIndexAccessor::store(dof_index.entityIndex(),
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
181 class SimplexGridP1Assembler
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(),
261 global_vertex_idx[n_j]);
262
263 contribution[dof_index] = 0.5;
264
265 GeometryIndexAccessor::store(dof_index.entityIndex(),
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(),
284 global_vertex_idx[n_j]);
285
286 contribution[dof_index] = 0.5;
287
288 GeometryIndexAccessor::store(dof_index.entityIndex(),
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:
313 typedef Dune::PDELab::HangingNodeManager<Grid,BoundaryFunction> HangingNodeManager;
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
Grid abstract base class.
Definition: grid.hh:375
static constexpr int dimension
The dimension of the grid.
Definition: grid.hh:387
static constexpr int mydimension
dimension of the intersection
Definition: intersection.hh:207
GridImp::template Codim< 0 >::Entity Entity
Type of entity that this Intersection belongs to.
Definition: intersection.hh:192
Dirichlet Constraints construction.
Definition: conforming.hh:38
Hanging Node constraints construction.
Definition: hangingnode.hh:311
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
A few common exception classes.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
Dune namespace.
Definition: alignedallocator.hh:13
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)