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 7 #include<dune/common/exceptions.hh> 8 #include<dune/geometry/referenceelements.hh> 9 #include<dune/geometry/type.hh> 10 #include<dune/pdelab/common/geometrywrapper.hh> 11 #include"conforming.hh" 12 #include"hangingnodemanager.hh" 13 14 namespace Dune { 15 namespace PDELab { 16 17 //! \addtogroup Constraints 18 //! \ingroup FiniteElementMap 19 //! \{ 20 21 class HangingNodesConstraintsAssemblers 22 { 23 public: 24 class CubeGridQ1Assembler 25 { 26 public: 27 template<typename IG, typename LFS, typename T, typename FlagVector> 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)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(), 104 GeometryTypes::vertex, 105 global_vertex_idx[fi[i+3]]); 106 107 contribution[dof_index] = 0.25; 108 109 GeometryIndexAccessor::store(dof_index.entityIndex(), 110 GeometryTypes::vertex, 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(), 119 GeometryTypes::vertex, 120 global_vertex_idx[fi[i+1]]); 121 122 contribution[dof_index] = 0.5; 123 124 GeometryIndexAccessor::store(dof_index.entityIndex(), 125 GeometryTypes::vertex, 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(), 134 GeometryTypes::vertex, 135 global_vertex_idx[fi[i+2]]); 136 137 contribution[dof_index] = 0.5; 138 139 GeometryIndexAccessor::store(dof_index.entityIndex(), 140 GeometryTypes::vertex, 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(), 160 GeometryTypes::vertex, 161 global_vertex_idx[n_j]); 162 163 contribution[dof_index] = 0.5; 164 165 GeometryIndexAccessor::store(dof_index.entityIndex(), 166 GeometryTypes::vertex, 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> 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)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(), 260 GeometryTypes::vertex, 261 global_vertex_idx[n_j]); 262 263 contribution[dof_index] = 0.5; 264 265 GeometryIndexAccessor::store(dof_index.entityIndex(), 266 GeometryTypes::vertex, 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(), 283 GeometryTypes::vertex, 284 global_vertex_idx[n_j]); 285 286 contribution[dof_index] = 0.5; 287 288 GeometryIndexAccessor::store(dof_index.entityIndex(), 289 GeometryTypes::vertex, 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 307 //! Hanging Node constraints construction 308 // works in any dimension and on all element types 309 template <class Grid, class HangingNodesConstraintsAssemblerType, class BoundaryFunction> 310 class HangingNodesDirichletConstraints : public ConformingDirichletConstraints 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 HangingNodesDirichletConstraints(Grid & grid,bool adaptToIsolatedHangingNodes,const BoundaryFunction & _boundaryFunction)322 HangingNodesDirichletConstraints( Grid & grid, 323 bool adaptToIsolatedHangingNodes, 324 const BoundaryFunction & _boundaryFunction ) 325 : manager(grid, _boundaryFunction) 326 { 327 if(adaptToIsolatedHangingNodes) 328 manager.adaptToIsolatedHangingNodes(); 329 } 330 update(Grid & grid)331 void update( Grid & grid ){ 332 manager.analyzeView(); 333 manager.adaptToIsolatedHangingNodes(); 334 } 335 336 337 //! skeleton constraints 338 /** 339 * \tparam I intersection geometry 340 * \tparam LFS local function space 341 * \tparam T TransformationType 342 */ 343 template<typename IG, typename LFS, typename T> skeleton(const IG & ig,const LFS & lfs_e,const LFS & lfs_f,T & trafo_e,T & trafo_f) const344 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 405 //! \} 406 407 } 408 } 409 410 #endif // DUNE_PDELAB_CONSTRAINTS_HANGINGNODE_HH 411