2 #ifndef DUNE_PDELAB_HANGINGNODECONSTRAINTS_HH
3 #define DUNE_PDELAB_HANGINGNODECONSTRAINTS_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::EntityPointer CellEntityPointer;
36 typedef typename Intersection::Entity Cell;
37 typedef typename Intersection::Geometry FaceGeometry;
38 typedef typename FaceGeometry::ctype DT;
39 typedef typename LFS::Traits::SizeType SizeType;
41 typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
42 const CellEntityPointer
e = ig.inside();
43 const CellEntityPointer f = ! ig.boundary() ? ig.outside() : ig.inside();
45 const std::size_t dimension = Intersection::dimension;
47 typedef Dune::ReferenceElement<DT,dimension> GRE;
48 const GRE& refelement_e = Dune::ReferenceElements<DT,dimension>::general(e->type());
49 const GRE& refelement_f = Dune::ReferenceElements<DT,dimension>::general(f->type());
53 if(e_has_hangingnodes && f_has_hangingnodes)
57 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
58 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
60 const Cell& cell = *(e_has_hangingnodes ? e : f);
61 const int faceindex = e_has_hangingnodes ? ig.indexInInside() : ig.indexInOutside();
62 const GRE & refelement = e_has_hangingnodes ? refelement_e : refelement_f;
63 const FlagVector & nodeState = e_has_hangingnodes ? nodeState_e : nodeState_f;
64 T & trafo = e_has_hangingnodes ? trafo_e : trafo_f;
68 std::vector<int> m(refelement.size(faceindex,1,dimension));
69 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++)
70 m[j] = refelement.subEntity(faceindex,1,j,dimension);
73 std::vector<std::size_t> global_vertex_idx(refelement.size(faceindex,1,dimension));
74 for (
int j=0; j<refelement.size(faceindex,1,dimension); ++j)
75 global_vertex_idx[j] = indexSet.subIndex(cell,refelement.subEntity(faceindex,1,j,dimension),dimension);
80 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
82 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
83 const GeometryType vertex_gt(0);
86 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
89 typename T::RowType contribution;
93 assert(nodeState.size() == 8);
95 const SizeType i = 4*j;
99 const unsigned int fi[16] = {0,1,2,3, 1,0,3,2, 2,0,3,1, 3,1,2,0};
102 if(nodeState[m[j]].isHanging()){
106 if(nodeState[m[fi[i+1]]].isHanging() && nodeState[m[fi[i+2]]].isHanging())
108 GeometryIndexAccessor::store(dof_index.entityIndex(),
110 global_vertex_idx[fi[i+3]]);
112 contribution[dof_index] = 0.25;
114 GeometryIndexAccessor::store(dof_index.entityIndex(),
116 global_vertex_idx[j]);
118 trafo[dof_index] = contribution;
121 else if(!nodeState[m[fi[i+1]]].isHanging())
123 GeometryIndexAccessor::store(dof_index.entityIndex(),
125 global_vertex_idx[fi[i+1]]);
127 contribution[dof_index] = 0.5;
129 GeometryIndexAccessor::store(dof_index.entityIndex(),
131 global_vertex_idx[j]);
133 trafo[dof_index] = contribution;
136 else if(!nodeState[m[fi[i+2]]].isHanging())
138 GeometryIndexAccessor::store(dof_index.entityIndex(),
140 global_vertex_idx[fi[i+2]]);
142 contribution[dof_index] = 0.5;
144 GeometryIndexAccessor::store(dof_index.entityIndex(),
146 global_vertex_idx[j]);
148 trafo[dof_index] = contribution;
152 }
else if(dimension == 2){
154 assert(nodeState.size() == 4);
158 if(nodeState[m[j]].isHanging()){
160 const SizeType n_j = 1-j;
162 assert( !nodeState[m[n_j]].isHanging() );
164 GeometryIndexAccessor::store(dof_index.entityIndex(),
166 global_vertex_idx[n_j]);
168 contribution[dof_index] = 0.5;
170 GeometryIndexAccessor::store(dof_index.entityIndex(),
172 global_vertex_idx[j]);
174 trafo[dof_index] = contribution;
189 template<
typename IG,
194 const FlagVector & nodeState_f,
195 const bool & e_has_hangingnodes,
196 const bool & f_has_hangingnodes,
197 const LFS & lfs_e,
const LFS & lfs_f,
198 T& trafo_e, T& trafo_f,
201 typedef IG Intersection;
202 typedef typename Intersection::EntityPointer CellEntityPointer;
203 typedef typename Intersection::Entity Cell;
204 typedef typename Intersection::Geometry FaceGeometry;
205 typedef typename FaceGeometry::ctype DT;
206 typedef typename LFS::Traits::SizeType SizeType;
207 typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
209 const CellEntityPointer
e = ig.inside();
210 const CellEntityPointer f = ! ig.boundary() ? ig.outside() : ig.inside();
212 const std::size_t dimension = Intersection::dimension;
214 typedef Dune::ReferenceElement<DT,dimension> GRE;
215 const GRE& refelement_e = Dune::ReferenceElements<DT,dimension>::general(e->type());
216 const GRE& refelement_f = Dune::ReferenceElements<DT,dimension>::general(f->type());
220 if(e_has_hangingnodes && f_has_hangingnodes)
224 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
225 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
227 const Cell& cell = *(e_has_hangingnodes ? e : f);
228 const int faceindex = e_has_hangingnodes ? ig.indexInInside() : ig.indexInOutside();
229 const GRE & refelement = e_has_hangingnodes ? refelement_e : refelement_f;
230 const FlagVector & nodeState = e_has_hangingnodes ? nodeState_e : nodeState_f;
231 T & trafo = e_has_hangingnodes ? trafo_e : trafo_f;
235 std::vector<int> m(refelement.size(faceindex,1,dimension));
236 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++)
237 m[j] = refelement.subEntity(faceindex,1,j,dimension);
240 std::vector<std::size_t> global_vertex_idx(refelement.size(faceindex,1,dimension));
241 for (
int j=0; j<refelement.size(faceindex,1,dimension); ++j)
242 global_vertex_idx[j] = indexSet.subIndex(cell,refelement.subEntity(faceindex,1,j,dimension),dimension);
247 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
249 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
250 const GeometryType vertex_gt(0);
253 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
256 typename T::RowType contribution;
260 assert(nodeState.size() == 4);
262 if(nodeState[m[j]].isHanging()){
263 for(
int k=1; k<=2; ++k ){
265 const SizeType n_j = (j+k)%3;
267 if( !nodeState[m[n_j]].isHanging() )
269 GeometryIndexAccessor::store(dof_index.entityIndex(),
271 global_vertex_idx[n_j]);
273 contribution[dof_index] = 0.5;
275 GeometryIndexAccessor::store(dof_index.entityIndex(),
277 global_vertex_idx[j]);
279 trafo[dof_index] = contribution;
283 }
else if(dimension == 2){
285 assert(nodeState.size() == 3);
287 if(nodeState[m[j]].isHanging()){
288 const SizeType n_j = 1-j;
289 assert( !nodeState[m[n_j]].isHanging() );
292 GeometryIndexAccessor::store(dof_index.entityIndex(),
294 global_vertex_idx[n_j]);
296 contribution[dof_index] = 0.5;
298 GeometryIndexAccessor::store(dof_index.entityIndex(),
300 global_vertex_idx[j]);
302 trafo[dof_index] = contribution;
319 template <
class Gr
id,
class HangingNodesConstra
intsAssemblerType,
class BoundaryFunction>
324 HangingNodeManager manager;
333 bool adaptToIsolatedHangingNodes,
334 const BoundaryFunction & _boundaryFunction )
335 : manager(grid, _boundaryFunction)
337 if(adaptToIsolatedHangingNodes)
353 template<
typename IG,
typename LFS,
typename T>
355 const LFS& lfs_e,
const LFS& lfs_f,
356 T& trafo_e, T& trafo_f)
const
358 typedef IG Intersection;
359 typedef typename Intersection::EntityPointer CellEntityPointer;
360 typedef typename Intersection::Geometry FaceGeometry;
361 typedef typename FaceGeometry::ctype DT;
363 const CellEntityPointer
e = ig.inside();
364 const CellEntityPointer f = ig.outside();
366 const Dune::ReferenceElement<DT,dimension>& refelem_e
367 = Dune::ReferenceElements<DT,dimension>::general(e->type());
368 const Dune::ReferenceElement<DT,dimension>& refelem_f
369 = Dune::ReferenceElements<DT,dimension>::general(f->type());
372 typedef typename std::vector<typename HangingNodeManager::NodeState> FlagVector;
373 const FlagVector isHangingNode_e(manager.
hangingNodes(*e));
374 const FlagVector isHangingNode_f(manager.
hangingNodes(*f));
378 assert(std::size_t(refelem_e.size(
dimension))
379 == isHangingNode_e.size());
380 assert(std::size_t(refelem_f.size(
dimension))
381 == isHangingNode_f.size());
384 const int faceindex_e = ig.indexInInside();
385 const int faceindex_f = ig.indexInOutside();
387 bool e_has_hangingnodes =
false;
389 for (
int j=0; j<refelem_e.size(faceindex_e,1,
dimension); j++){
390 const int index = refelem_e.subEntity(faceindex_e,1,j,
dimension);
391 if(isHangingNode_e[index].isHanging())
393 e_has_hangingnodes =
true;
398 bool f_has_hangingnodes =
false;
400 for (
int j=0; j<refelem_f.size(faceindex_f,1,
dimension); j++){
401 const int index = refelem_f.subEntity(faceindex_f,1,j,
dimension);
402 if(isHangingNode_f[index].isHanging())
404 f_has_hangingnodes =
true;
410 if(! e_has_hangingnodes && ! f_has_hangingnodes)
413 HangingNodesConstraintsAssemblerType::
414 assembleConstraints(isHangingNode_e, isHangingNode_f,
415 e_has_hangingnodes, f_has_hangingnodes,
void update(Grid &grid)
Definition: hangingnode.hh:341
Definition: hangingnode.hh:330
Hanging Node constraints construction.
Definition: hangingnode.hh:320
Definition: hangingnode.hh:186
HangingNodesDirichletConstraints(Grid &grid, bool adaptToIsolatedHangingNodes, const BoundaryFunction &_boundaryFunction)
Definition: hangingnode.hh:332
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:354
const std::vector< NodeState > hangingNodes(const Cell &e) const
Definition: hangingnodemanager.hh:232
void analyzeView()
Definition: hangingnodemanager.hh:103
const E & e
Definition: interpolate.hh:172
Definition: hangingnode.hh:329
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
Definition: hangingnode.hh:328
Dirichlet Constraints construction.
Definition: conforming.hh:36
const IG & ig
Definition: constraints.hh:147
Definition: hangingnode.hh:21
Definition: adaptivity.hh:27
Definition: hangingnodemanager.hh:17
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:193
void adaptToIsolatedHangingNodes()
Definition: hangingnodemanager.hh:278
Definition: hangingnode.hh:327
Definition: hangingnode.hh:24