2 #ifndef DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH
3 #define DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/geometry/referenceelements.hh>
9 #include<dune/geometry/type.hh>
27 template<
typename IG,
typename LFS,
typename T,
typename FlagVector>
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,
34 typedef IG Intersection;
35 typedef typename Intersection::Entity Cell;
36 typedef typename LFS::Traits::SizeType SizeType;
38 typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
40 auto f = !
ig.boundary() ?
ig.outside() :
ig.inside();
42 const std::size_t dimension = Intersection::mydimension;
44 auto refelement_e = referenceElement(
e.geometry());
45 auto refelement_f = referenceElement(f.geometry());
49 if(e_has_hangingnodes && f_has_hangingnodes)
53 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
54 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
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;
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);
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);
76 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
78 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
79 const GeometryType vertex_gt(0);
82 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
85 typename T::RowType contribution;
89 assert(nodeState.size() == 8);
91 const SizeType i = 4*j;
95 const unsigned int fi[16] = {0,1,2,3, 1,0,3,2, 2,0,3,1, 3,1,2,0};
98 if(nodeState[m[j]].isHanging()){
102 if(nodeState[m[fi[i+1]]].isHanging() && nodeState[m[fi[i+2]]].isHanging())
104 GeometryIndexAccessor::store(dof_index.entityIndex(),
106 global_vertex_idx[fi[i+3]]);
108 contribution[dof_index] = 0.25;
110 GeometryIndexAccessor::store(dof_index.entityIndex(),
112 global_vertex_idx[j]);
114 trafo[dof_index] = contribution;
117 else if(!nodeState[m[fi[i+1]]].isHanging())
119 GeometryIndexAccessor::store(dof_index.entityIndex(),
121 global_vertex_idx[fi[i+1]]);
123 contribution[dof_index] = 0.5;
125 GeometryIndexAccessor::store(dof_index.entityIndex(),
127 global_vertex_idx[j]);
129 trafo[dof_index] = contribution;
132 else if(!nodeState[m[fi[i+2]]].isHanging())
134 GeometryIndexAccessor::store(dof_index.entityIndex(),
136 global_vertex_idx[fi[i+2]]);
138 contribution[dof_index] = 0.5;
140 GeometryIndexAccessor::store(dof_index.entityIndex(),
142 global_vertex_idx[j]);
144 trafo[dof_index] = contribution;
148 }
else if(dimension == 2){
150 assert(nodeState.size() == 4);
154 if(nodeState[m[j]].isHanging()){
156 const SizeType n_j = 1-j;
158 assert( !nodeState[m[n_j]].isHanging() );
160 GeometryIndexAccessor::store(dof_index.entityIndex(),
162 global_vertex_idx[n_j]);
164 contribution[dof_index] = 0.5;
166 GeometryIndexAccessor::store(dof_index.entityIndex(),
168 global_vertex_idx[j]);
170 trafo[dof_index] = contribution;
185 template<
typename IG,
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,
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;
202 auto e =
ig.inside();
203 auto f = !
ig.boundary() ?
ig.outside() :
ig.inside();
205 const std::size_t dimension = Intersection::mydimension;
207 auto refelement_e = referenceElement(
e.geometry());
208 auto refelement_f = referenceElement(f.geometry());
212 if(e_has_hangingnodes && f_has_hangingnodes)
216 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
217 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
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;
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);
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);
239 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
241 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
242 const GeometryType vertex_gt(0);
245 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
248 typename T::RowType contribution;
252 assert(nodeState.size() == 4);
254 if(nodeState[m[j]].isHanging()){
255 for(
int k=1; k<=2; ++k ){
257 const SizeType n_j = (j+k)%3;
259 if( !nodeState[m[n_j]].isHanging() )
261 GeometryIndexAccessor::store(dof_index.entityIndex(),
263 global_vertex_idx[n_j]);
265 contribution[dof_index] = 0.5;
267 GeometryIndexAccessor::store(dof_index.entityIndex(),
269 global_vertex_idx[j]);
271 trafo[dof_index] = contribution;
275 }
else if(dimension == 2){
277 assert(nodeState.size() == 3);
279 if(nodeState[m[j]].isHanging()){
280 const SizeType n_j = 1-j;
281 assert( !nodeState[m[n_j]].isHanging() );
284 GeometryIndexAccessor::store(dof_index.entityIndex(),
286 global_vertex_idx[n_j]);
288 contribution[dof_index] = 0.5;
290 GeometryIndexAccessor::store(dof_index.entityIndex(),
292 global_vertex_idx[j]);
294 trafo[dof_index] = contribution;
311 template <
class Gr
id,
class HangingNodesConstra
intsAssemblerType,
class BoundaryFunction>
325 bool adaptToIsolatedHangingNodes,
326 const BoundaryFunction & _boundaryFunction )
327 : manager(grid, _boundaryFunction)
329 if(adaptToIsolatedHangingNodes)
345 template<
typename IG,
typename LFS,
typename T>
347 const LFS& lfs_e,
const LFS& lfs_f,
348 T& trafo_e, T& trafo_f)
const
350 auto e =
ig.inside();
351 auto f =
ig.outside();
353 auto refelem_e = referenceElement(
e.geometry());
354 auto refelem_f = referenceElement(f.geometry());
357 typedef typename std::vector<typename HangingNodeManager::NodeState> FlagVector;
359 const FlagVector isHangingNode_f(manager.
hangingNodes(f));
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());
369 const int faceindex_e =
ig.indexInInside();
370 const int faceindex_f =
ig.indexInOutside();
372 bool e_has_hangingnodes =
false;
374 for (
int j=0; j<refelem_e.size(faceindex_e,1,
dimension); j++){
376 if(isHangingNode_e[
index].isHanging())
378 e_has_hangingnodes =
true;
383 bool f_has_hangingnodes =
false;
385 for (
int j=0; j<refelem_f.size(faceindex_f,1,
dimension); j++){
387 if(isHangingNode_f[
index].isHanging())
389 f_has_hangingnodes =
true;
395 if(! e_has_hangingnodes && ! f_has_hangingnodes)
398 HangingNodesConstraintsAssemblerType::
399 assembleConstraints(isHangingNode_e, isHangingNode_f,
400 e_has_hangingnodes, f_has_hangingnodes,
412 #endif // DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH