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 Intersection::Geometry FaceGeometry;
37 typedef typename FaceGeometry::ctype DT;
38 typedef typename LFS::Traits::SizeType SizeType;
40 typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::IndexSet IndexSet;
42 auto f = ! ig.boundary() ? ig.outside() : ig.inside();
44 const std::size_t dimension = Intersection::dimension;
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());
52 if(e_has_hangingnodes && f_has_hangingnodes)
56 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
57 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
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;
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);
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);
79 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
81 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
82 const GeometryType vertex_gt(0);
85 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
88 typename T::RowType contribution;
92 assert(nodeState.size() == 8);
94 const SizeType i = 4*j;
98 const unsigned int fi[16] = {0,1,2,3, 1,0,3,2, 2,0,3,1, 3,1,2,0};
101 if(nodeState[m[j]].isHanging()){
105 if(nodeState[m[fi[i+1]]].isHanging() && nodeState[m[fi[i+2]]].isHanging())
107 GeometryIndexAccessor::store(dof_index.entityIndex(),
109 global_vertex_idx[fi[i+3]]);
111 contribution[dof_index] = 0.25;
113 GeometryIndexAccessor::store(dof_index.entityIndex(),
115 global_vertex_idx[j]);
117 trafo[dof_index] = contribution;
120 else if(!nodeState[m[fi[i+1]]].isHanging())
122 GeometryIndexAccessor::store(dof_index.entityIndex(),
124 global_vertex_idx[fi[i+1]]);
126 contribution[dof_index] = 0.5;
128 GeometryIndexAccessor::store(dof_index.entityIndex(),
130 global_vertex_idx[j]);
132 trafo[dof_index] = contribution;
135 else if(!nodeState[m[fi[i+2]]].isHanging())
137 GeometryIndexAccessor::store(dof_index.entityIndex(),
139 global_vertex_idx[fi[i+2]]);
141 contribution[dof_index] = 0.5;
143 GeometryIndexAccessor::store(dof_index.entityIndex(),
145 global_vertex_idx[j]);
147 trafo[dof_index] = contribution;
151 }
else if(dimension == 2){
153 assert(nodeState.size() == 4);
157 if(nodeState[m[j]].isHanging()){
159 const SizeType n_j = 1-j;
161 assert( !nodeState[m[n_j]].isHanging() );
163 GeometryIndexAccessor::store(dof_index.entityIndex(),
165 global_vertex_idx[n_j]);
167 contribution[dof_index] = 0.5;
169 GeometryIndexAccessor::store(dof_index.entityIndex(),
171 global_vertex_idx[j]);
173 trafo[dof_index] = contribution;
188 template<
typename IG,
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,
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;
207 auto e = ig.inside();
208 auto f = ! ig.boundary() ? ig.outside() : ig.inside();
210 const std::size_t dimension = Intersection::dimension;
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());
218 if(e_has_hangingnodes && f_has_hangingnodes)
222 const LFS & lfs = e_has_hangingnodes ? lfs_e : lfs_f;
223 const IndexSet& indexSet = lfs.gridFunctionSpace().gridView().indexSet();
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;
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);
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);
245 typename LFS::Traits::DOFIndex dof_index(lfs.dofIndex(0));
247 typedef typename LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::GeometryIndex GeometryIndexAccessor;
248 const GeometryType vertex_gt(0);
251 for (
int j=0; j<refelement.size(faceindex,1,dimension); j++){
254 typename T::RowType contribution;
258 assert(nodeState.size() == 4);
260 if(nodeState[m[j]].isHanging()){
261 for(
int k=1; k<=2; ++k ){
263 const SizeType n_j = (j+k)%3;
265 if( !nodeState[m[n_j]].isHanging() )
267 GeometryIndexAccessor::store(dof_index.entityIndex(),
269 global_vertex_idx[n_j]);
271 contribution[dof_index] = 0.5;
273 GeometryIndexAccessor::store(dof_index.entityIndex(),
275 global_vertex_idx[j]);
277 trafo[dof_index] = contribution;
281 }
else if(dimension == 2){
283 assert(nodeState.size() == 3);
285 if(nodeState[m[j]].isHanging()){
286 const SizeType n_j = 1-j;
287 assert( !nodeState[m[n_j]].isHanging() );
290 GeometryIndexAccessor::store(dof_index.entityIndex(),
292 global_vertex_idx[n_j]);
294 contribution[dof_index] = 0.5;
296 GeometryIndexAccessor::store(dof_index.entityIndex(),
298 global_vertex_idx[j]);
300 trafo[dof_index] = contribution;
317 template <
class Gr
id,
class HangingNodesConstra
intsAssemblerType,
class BoundaryFunction>
322 HangingNodeManager manager;
325 enum { doBoundary =
true };
326 enum { doSkeleton =
true };
327 enum { doVolume =
false };
328 enum { dimension = Grid::dimension };
331 bool adaptToIsolatedHangingNodes,
332 const BoundaryFunction & _boundaryFunction )
333 : manager(grid, _boundaryFunction)
335 if(adaptToIsolatedHangingNodes)
351 template<
typename IG,
typename LFS,
typename T>
353 const LFS& lfs_e,
const LFS& lfs_f,
354 T& trafo_e, T& trafo_f)
const 356 typedef IG Intersection;
357 typedef typename Intersection::Geometry FaceGeometry;
358 typedef typename FaceGeometry::ctype DT;
360 auto e = ig.inside();
361 auto f = ig.outside();
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());
369 typedef typename std::vector<typename HangingNodeManager::NodeState> FlagVector;
371 const FlagVector isHangingNode_f(manager.
hangingNodes(f));
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());
381 const int faceindex_e = ig.indexInInside();
382 const int faceindex_f = ig.indexInOutside();
384 bool e_has_hangingnodes =
false;
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())
390 e_has_hangingnodes =
true;
395 bool f_has_hangingnodes =
false;
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())
401 f_has_hangingnodes =
true;
407 if(! e_has_hangingnodes && ! f_has_hangingnodes)
410 HangingNodesConstraintsAssemblerType::
411 assembleConstraints(isHangingNode_e, isHangingNode_f,
412 e_has_hangingnodes, f_has_hangingnodes,
424 #endif // DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH const IG & ig
Definition: constraints.hh:148
Definition: hangingnode.hh:185
void update(Grid &grid)
Definition: hangingnode.hh:339
Definition: hangingnodemanager.hh:17
Definition: hangingnode.hh:24
void adaptToIsolatedHangingNodes()
Definition: hangingnodemanager.hh:272
HangingNodesDirichletConstraints(Grid &grid, bool adaptToIsolatedHangingNodes, const BoundaryFunction &_boundaryFunction)
Definition: hangingnode.hh:330
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
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
void analyzeView()
Definition: hangingnodemanager.hh:101
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const Entity & e
Definition: localfunctionspace.hh:111
Hanging Node constraints construction.
Definition: hangingnode.hh:318
Definition: hangingnode.hh:21
const std::vector< NodeState > hangingNodes(const Cell &e) const
Definition: hangingnodemanager.hh:226
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
Dirichlet Constraints construction.
Definition: conforming.hh:36