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;
81 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
84 typename T::RowType contribution;
88 assert(nodeState.size() == 8);
90 const SizeType i = 4*j;
94 const unsigned int fi[16] = {0,1,2,3, 1,0,3,2, 2,0,3,1, 3,1,2,0};
97 if(nodeState[m[j]].isHanging()){
101 if(nodeState[m[fi[i+1]]].isHanging() && nodeState[m[fi[i+2]]].isHanging())
103 GeometryIndexAccessor::store(dof_index.entityIndex(),
104 GeometryTypes::vertex,
105 global_vertex_idx[fi[i+3]]);
107 contribution[dof_index] = 0.25;
109 GeometryIndexAccessor::store(dof_index.entityIndex(),
110 GeometryTypes::vertex,
111 global_vertex_idx[j]);
113 trafo[dof_index] = contribution;
116 else if(!nodeState[m[fi[i+1]]].isHanging())
118 GeometryIndexAccessor::store(dof_index.entityIndex(),
119 GeometryTypes::vertex,
120 global_vertex_idx[fi[i+1]]);
122 contribution[dof_index] = 0.5;
124 GeometryIndexAccessor::store(dof_index.entityIndex(),
125 GeometryTypes::vertex,
126 global_vertex_idx[j]);
128 trafo[dof_index] = contribution;
131 else if(!nodeState[m[fi[i+2]]].isHanging())
133 GeometryIndexAccessor::store(dof_index.entityIndex(),
134 GeometryTypes::vertex,
135 global_vertex_idx[fi[i+2]]);
137 contribution[dof_index] = 0.5;
139 GeometryIndexAccessor::store(dof_index.entityIndex(),
140 GeometryTypes::vertex,
141 global_vertex_idx[j]);
143 trafo[dof_index] = contribution;
147 }
else if(dimension == 2){
149 assert(nodeState.size() == 4);
153 if(nodeState[m[j]].isHanging()){
155 const SizeType n_j = 1-j;
157 assert( !nodeState[m[n_j]].isHanging() );
159 GeometryIndexAccessor::store(dof_index.entityIndex(),
160 GeometryTypes::vertex,
161 global_vertex_idx[n_j]);
163 contribution[dof_index] = 0.5;
165 GeometryIndexAccessor::store(dof_index.entityIndex(),
166 GeometryTypes::vertex,
167 global_vertex_idx[j]);
169 trafo[dof_index] = contribution;
184 template<
typename IG,
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,
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;
201 auto e =
ig.inside();
202 auto f = !
ig.boundary() ?
ig.outside() :
ig.inside();
204 const std::size_t dimension = Intersection::mydimension;
206 auto refelement_e = referenceElement(
e.geometry());
207 auto refelement_f = referenceElement(f.geometry());
211 if(e_has_hangingnodes && f_has_hangingnodes)
215 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
216 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
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;
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);
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);
238 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
240 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
243 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
246 typename T::RowType contribution;
250 assert(nodeState.size() == 4);
252 if(nodeState[m[j]].isHanging()){
253 for(
int k=1; k<=2; ++k ){
255 const SizeType n_j = (j+k)%3;
257 if( !nodeState[m[n_j]].isHanging() )
259 GeometryIndexAccessor::store(dof_index.entityIndex(),
260 GeometryTypes::vertex,
261 global_vertex_idx[n_j]);
263 contribution[dof_index] = 0.5;
265 GeometryIndexAccessor::store(dof_index.entityIndex(),
266 GeometryTypes::vertex,
267 global_vertex_idx[j]);
269 trafo[dof_index] = contribution;
273 }
else if(dimension == 2){
275 assert(nodeState.size() == 3);
277 if(nodeState[m[j]].isHanging()){
278 const SizeType n_j = 1-j;
279 assert( !nodeState[m[n_j]].isHanging() );
282 GeometryIndexAccessor::store(dof_index.entityIndex(),
283 GeometryTypes::vertex,
284 global_vertex_idx[n_j]);
286 contribution[dof_index] = 0.5;
288 GeometryIndexAccessor::store(dof_index.entityIndex(),
289 GeometryTypes::vertex,
290 global_vertex_idx[j]);
292 trafo[dof_index] = contribution;
309 template <
class Gr
id,
class HangingNodesConstra
intsAssemblerType,
class BoundaryFunction>
323 bool adaptToIsolatedHangingNodes,
324 const BoundaryFunction & _boundaryFunction )
325 : manager(grid, _boundaryFunction)
327 if(adaptToIsolatedHangingNodes)
343 template<
typename IG,
typename LFS,
typename T>
345 const LFS& lfs_e,
const LFS& lfs_f,
346 T& trafo_e, T& trafo_f)
const
348 auto e =
ig.inside();
349 auto f =
ig.outside();
351 auto refelem_e = referenceElement(
e.geometry());
352 auto refelem_f = referenceElement(f.geometry());
355 typedef typename std::vector<typename HangingNodeManager::NodeState> FlagVector;
357 const FlagVector isHangingNode_f(manager.
hangingNodes(f));
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());
367 const int faceindex_e =
ig.indexInInside();
368 const int faceindex_f =
ig.indexInOutside();
370 bool e_has_hangingnodes =
false;
372 for (
int j=0; j<refelem_e.size(faceindex_e,1,
dimension); j++){
374 if(isHangingNode_e[
index].isHanging())
376 e_has_hangingnodes =
true;
381 bool f_has_hangingnodes =
false;
383 for (
int j=0; j<refelem_f.size(faceindex_f,1,
dimension); j++){
385 if(isHangingNode_f[
index].isHanging())
387 f_has_hangingnodes =
true;
393 if(! e_has_hangingnodes && ! f_has_hangingnodes)
396 HangingNodesConstraintsAssemblerType::
397 assembleConstraints(isHangingNode_e, isHangingNode_f,
398 e_has_hangingnodes, f_has_hangingnodes,
410 #endif // DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH