1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_UGGRID_INDEXSETS_HH 4 #define DUNE_UGGRID_INDEXSETS_HH 5 6 /** \file 7 \brief The index and id sets for the UGGrid class 8 */ 9 10 #include <vector> 11 #include <set> 12 13 #include <dune/common/std/type_traits.hh> 14 #include <dune/grid/common/grid.hh> 15 16 namespace Dune { 17 18 template<class GridImp> 19 class UGGridLevelIndexSet : public IndexSet<GridImp,UGGridLevelIndexSet<GridImp>, UG::UINT> 20 { 21 enum {dim = GridImp::dimension}; 22 23 public: 24 25 /** \brief Default constructor 26 27 Unfortunately we can't force the user to init grid_ and level_, because 28 UGGridLevelIndexSets are meant to be stored in an array. 29 30 \todo I want to make this constructor private, but I can't, because 31 it is called by UGGrid through a std::vector::resize() 32 */ UGGridLevelIndexSet()33 UGGridLevelIndexSet () 34 : level_(0), 35 numSimplices_(0), 36 numPyramids_(0), 37 numPrisms_(0), 38 numCubes_(0), 39 numVertices_(0), 40 numEdges_(0), 41 numTriFaces_(0), 42 numQuadFaces_(0) 43 {} 44 45 //! get index of an entity 46 template<int cd> index(const typename GridImp::Traits::template Codim<cd>::Entity & e) const47 unsigned int index (const typename GridImp::Traits::template Codim<cd>::Entity& e) const 48 { 49 return UG_NS<dim>::levelIndex(e.impl().getTarget()); 50 } 51 52 /** \brief Get index of subEntity of a codim cc entity 53 * \param codim Codimension WITH RESPECT TO THE GRID of the subentity whose index we want 54 */ 55 template<int cc> subIndex(const typename GridImp::Traits::template Codim<cc>::Entity & entity,int i,unsigned int codim) const56 unsigned int subIndex (const typename GridImp::Traits::template Codim<cc>::Entity& entity, 57 int i, 58 unsigned int codim) const 59 { 60 // The entity is a vertex, so each subentity must be a vertex too (anything else is not supported) 61 if (cc==dim) 62 { 63 assert(codim==dim); 64 return UG_NS<dim>::levelIndex(entity.impl().getTarget()); 65 } 66 67 // The following block is for 2d grids 68 if constexpr (dim==2) 69 { 70 // The entity is an element 71 if constexpr (cc==0) 72 { 73 // Element indices 74 if (codim==0) 75 return UG_NS<dim>::levelIndex(entity.impl().getTarget()); 76 77 // Edge indices 78 if (codim==1) 79 { 80 auto ref_el = referenceElement<double,dim>(entity.type()); 81 auto a = ref_el.subEntity(i,dim-1,0,dim); 82 auto b = ref_el.subEntity(i,dim-1,1,dim); 83 return UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(entity.impl().getTarget(), 84 UGGridRenumberer<dim>::verticesDUNEtoUG(a,entity.type())), 85 UG_NS<dim>::Corner(entity.impl().getTarget(), 86 UGGridRenumberer<dim>::verticesDUNEtoUG(b,entity.type())))); 87 } 88 89 // Vertex indices 90 if (codim==dim) 91 return UG_NS<dim>::levelIndex(UG_NS<dim>::Corner(entity.impl().getTarget(), 92 UGGridRenumberer<dim>::verticesDUNEtoUG(i,entity.type()))); 93 } 94 95 // The entity is an edge 96 if constexpr (cc==1) 97 return UG_NS<dim>::levelIndex(entity.impl().getTarget()->links[i].nbnode); 98 } 99 100 // The following block is for 3d grids 101 if constexpr (dim==3) 102 { 103 // The entity is an element 104 if constexpr (cc==0) 105 { 106 // Element indices 107 if (codim==0) 108 return UG_NS<dim>::levelIndex(entity.impl().getTarget()); 109 110 // Face indices 111 if (codim==1) 112 return UG_NS<dim>::levelIndex(UG_NS<dim>::SideVector(entity.impl().getTarget(), 113 UGGridRenumberer<dim>::facesDUNEtoUG(i,entity.type()))); 114 115 // Edge indices 116 if (codim==2) 117 { 118 auto ref_el = referenceElement<double,dim>(entity.type()); 119 auto a = ref_el.subEntity(i,dim-1,0,dim); 120 auto b = ref_el.subEntity(i,dim-1,1,dim); 121 return UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(entity.impl().getTarget(), 122 UGGridRenumberer<dim>::verticesDUNEtoUG(a,entity.type())), 123 UG_NS<dim>::Corner(entity.impl().getTarget(), 124 UGGridRenumberer<dim>::verticesDUNEtoUG(b,entity.type())))); 125 } 126 127 // Vertex indices 128 if (codim==3) 129 return UG_NS<dim>::levelIndex(UG_NS<dim>::Corner(entity.impl().getTarget(), 130 UGGridRenumberer<dim>::verticesDUNEtoUG(i,entity.type()))); 131 } 132 133 // The entity is a face 134 if constexpr (cc==1) 135 { 136 DUNE_THROW(NotImplemented, "Subindices of an element face"); 137 } 138 139 // The entity is an edge 140 if constexpr (cc==2) 141 return UG_NS<dim>::levelIndex(entity.impl().getTarget()->links[i].nbnode); 142 } 143 144 // Should never happen 145 return std::numeric_limits<unsigned int>::max(); 146 } 147 148 //! get number of entities of given codim, type and on this level size(int codim) const149 std::size_t size (int codim) const { 150 if (codim==0) 151 return numSimplices_+numPyramids_+numPrisms_+numCubes_; 152 if (codim==dim) 153 return numVertices_; 154 if (codim==dim-1) 155 return numEdges_; 156 if (codim==1) 157 return numTriFaces_+numQuadFaces_; 158 DUNE_THROW(NotImplemented, "wrong codim!"); 159 } 160 161 //! get number of entities of given codim, type and on this level size(GeometryType type) const162 std::size_t size (GeometryType type) const 163 { 164 int codim = GridImp::dimension-type.dim(); 165 166 if (codim==0) { 167 if (type.isSimplex()) 168 return numSimplices_; 169 else if (type.isPyramid()) 170 return numPyramids_; 171 else if (type.isPrism()) 172 return numPrisms_; 173 else if (type.isCube()) 174 return numCubes_; 175 else 176 return 0; 177 178 } 179 180 if (codim==dim) { 181 return numVertices_; 182 } 183 if (codim==dim-1) { 184 return numEdges_; 185 } 186 if (codim==1) { 187 if (type.isSimplex()) 188 return numTriFaces_; 189 else if (type.isCube()) 190 return numQuadFaces_; 191 else 192 return 0; 193 } 194 195 DUNE_THROW(NotImplemented, "Wrong codim!"); 196 } 197 types(int codim) const198 std::vector< GeometryType > types ( int codim ) const { return myTypes_[ codim ]; } 199 200 /** \brief Deliver all geometry types used in this grid */ geomTypes(int codim) const201 const std::vector<GeometryType>& geomTypes (int codim) const 202 { 203 return myTypes_[codim]; 204 } 205 206 /** \brief Return true if e is contained in the index set. 207 208 \note If 'entity' is from another grid this method may still return 'true'. 209 This is acceptable by the Dune grid interface specification. 210 */ 211 template <class EntityType> contains(const EntityType & entity) const212 bool contains (const EntityType& entity) const 213 { 214 return entity.level() == level_; 215 } 216 217 /** \brief Update the level indices. This method is called after each grid change */ 218 void update(const GridImp& grid, int level, std::vector<unsigned int>* nodePermutation=0); 219 220 const GridImp* grid_; 221 int level_; 222 223 int numSimplices_; 224 int numPyramids_; 225 int numPrisms_; 226 int numCubes_; 227 int numVertices_; 228 int numEdges_; 229 int numTriFaces_; 230 int numQuadFaces_; 231 232 std::vector<GeometryType> myTypes_[dim+1]; 233 }; 234 235 template<class GridImp> 236 class UGGridLeafIndexSet : public IndexSet<GridImp,UGGridLeafIndexSet<GridImp>, UG::UINT> 237 { 238 public: 239 240 /* 241 We use the remove_const to extract the Type from the mutable class, 242 because the const class is not instantiated yet. 243 */ 244 enum {dim = std::remove_const<GridImp>::type::dimension}; 245 246 //! constructor stores reference to a grid and level UGGridLeafIndexSet(const GridImp & g)247 UGGridLeafIndexSet (const GridImp& g) 248 : grid_(g), coarsestLevelWithLeafElements_(0) 249 {} 250 251 //! get index of an entity 252 /* 253 We use the remove_const to extract the Type from the mutable class, 254 because the const class is not instantiated yet. 255 */ 256 template<int cd> index(const typename std::remove_const<GridImp>::type::Traits::template Codim<cd>::Entity & e) const257 int index (const typename std::remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 258 { 259 return UG_NS<dim>::leafIndex(e.impl().getTarget()); 260 } 261 262 //! get index of subEntity of a codim 0 entity 263 /* 264 We use the remove_const to extract the Type from the mutable class, 265 because the const class is not instantiated yet. 266 */ 267 template<int cc> subIndex(const typename std::remove_const<GridImp>::type::Traits::template Codim<cc>::Entity & entity,int i,unsigned int codim) const268 unsigned int subIndex (const typename std::remove_const<GridImp>::type::Traits::template Codim<cc>::Entity& entity, 269 int i, 270 unsigned int codim) const 271 { 272 // The entity is a vertex, so each subentity must be a vertex too (anything else is not supported) 273 if (cc==dim) 274 { 275 assert(codim==dim); 276 return UG_NS<dim>::leafIndex(entity.impl().getTarget()); 277 } 278 279 // The following block is for 2d grids 280 if constexpr (dim==2) 281 { 282 // The entity is an element 283 if constexpr (cc==0) 284 { 285 // Element indices 286 if (codim==0) 287 return UG_NS<dim>::leafIndex(entity.impl().getTarget()); 288 289 // Edge indices 290 if (codim==1) 291 { 292 auto ref_el = referenceElement<double,dim>(entity.type()); 293 auto a = ref_el.subEntity(i,dim-1,0,dim); 294 auto b = ref_el.subEntity(i,dim-1,1,dim); 295 return UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(entity.impl().getTarget(), 296 UGGridRenumberer<dim>::verticesDUNEtoUG(a,entity.type())), 297 UG_NS<dim>::Corner(entity.impl().getTarget(), 298 UGGridRenumberer<dim>::verticesDUNEtoUG(b,entity.type())))); 299 } 300 301 // Vertex indices 302 if (codim==dim) 303 return UG_NS<dim>::leafIndex(UG_NS<dim>::Corner(entity.impl().getTarget(), 304 UGGridRenumberer<dim>::verticesDUNEtoUG(i,entity.type()))); 305 } 306 307 // The entity is an edge 308 if constexpr (cc==1) 309 return UG_NS<dim>::leafIndex(entity.impl().getTarget()->links[i].nbnode); 310 } 311 312 // The following block is for 3d grids 313 if constexpr (dim==3) 314 { 315 // The entity is an element 316 if constexpr (cc==0) 317 { 318 // Element indices 319 if (codim==0) 320 return UG_NS<dim>::leafIndex(entity.impl().getTarget()); 321 322 // Face indices 323 if (codim==1) 324 return UG_NS<dim>::leafIndex(UG_NS<dim>::SideVector(entity.impl().getTarget(), 325 UGGridRenumberer<dim>::facesDUNEtoUG(i,entity.type()))); 326 327 // Edge indices 328 if (codim==2) 329 { 330 auto ref_el = referenceElement<double,dim>(entity.type()); 331 auto a = ref_el.subEntity(i,dim-1,0,dim); 332 auto b = ref_el.subEntity(i,dim-1,1,dim); 333 return UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(entity.impl().getTarget(), 334 UGGridRenumberer<dim>::verticesDUNEtoUG(a,entity.type())), 335 UG_NS<dim>::Corner(entity.impl().getTarget(), 336 UGGridRenumberer<dim>::verticesDUNEtoUG(b,entity.type())))); 337 } 338 339 // Vertex indices 340 if (codim==3) 341 return UG_NS<dim>::leafIndex(UG_NS<dim>::Corner(entity.impl().getTarget(), 342 UGGridRenumberer<dim>::verticesDUNEtoUG(i,entity.type()))); 343 } 344 345 // The entity is a face 346 if constexpr (cc==1) 347 { 348 DUNE_THROW(NotImplemented, "Subindices of an element face"); 349 } 350 351 // The entity is an edge 352 if constexpr (cc==2) 353 return UG_NS<dim>::leafIndex(entity.impl().getTarget()->links[i].nbnode); 354 } 355 356 return std::numeric_limits<unsigned int>::max(); 357 } 358 359 //! get number of entities of given codim and type size(GeometryType type) const360 std::size_t size (GeometryType type) const 361 { 362 if (type.dim()==GridImp::dimension) { 363 if (type.isSimplex()) 364 return numSimplices_; 365 else if (type.isPyramid()) 366 return numPyramids_; 367 else if (type.isPrism()) 368 return numPrisms_; 369 else if (type.isCube()) 370 return numCubes_; 371 else 372 return 0; 373 } 374 if (type.dim()==0) { 375 return numVertices_; 376 } 377 if (type.dim()==1) { 378 return numEdges_; 379 } 380 if (type.isTriangle()) 381 return numTriFaces_; 382 else if (type.isQuadrilateral()) 383 return numQuadFaces_; 384 385 return 0; 386 } 387 388 //! get number of entities of given codim size(int codim) const389 std::size_t size (int codim) const 390 { 391 int s=0; 392 const std::vector<GeometryType>& geomTs = geomTypes(codim); 393 for (unsigned int i=0; i<geomTs.size(); ++i) 394 s += size(geomTs[i]); 395 return s; 396 } 397 types(int codim) const398 std::vector< GeometryType > types ( int codim ) const { return myTypes_[ codim ]; } 399 400 /** deliver all geometry types used in this grid */ geomTypes(int codim) const401 const std::vector<GeometryType>& geomTypes (int codim) const 402 { 403 return myTypes_[codim]; 404 } 405 406 /** \brief Return true if e is contained in the index set. 407 408 \note If 'entity' is from another grid this method may still return 'true'. 409 This is acceptable by the Dune grid interface specification. 410 */ 411 template <class EntityType> contains(const EntityType & entity) const412 bool contains (const EntityType& entity) const 413 { 414 return UG_NS<dim>::isLeaf(entity.impl().getTarget()); 415 } 416 417 418 /** \brief Update the leaf indices. This method is called after each grid change. */ 419 void update(std::vector<unsigned int>* nodePermutation=0); 420 421 const GridImp& grid_; 422 423 /** \brief The lowest level that contains leaf elements 424 425 This corresponds to UG's fullRefineLevel, which is, unfortunately only 426 computed if you use some nontrivial UG algebra. Thus we compute it 427 ourselves, and use it to speed up the leaf iterators. 428 */ 429 unsigned int coarsestLevelWithLeafElements_; 430 431 432 433 int numSimplices_; 434 int numPyramids_; 435 int numPrisms_; 436 int numCubes_; 437 int numVertices_; 438 int numEdges_; 439 int numTriFaces_; 440 int numQuadFaces_; 441 442 std::vector<GeometryType> myTypes_[dim+1]; 443 }; 444 445 446 /** \brief Implementation class for the UGGrid Id sets 447 448 The UGGridGlobalIdSet and the UGGridLocalIdSet are identical. This 449 class implements them both at once. 450 */ 451 template <class GridImp> 452 class UGGridIdSet : public IdSet<GridImp,UGGridIdSet<GridImp>,typename UG_NS<std::remove_const<GridImp>::type::dimension>::UG_ID_TYPE> 453 { 454 enum {dim = std::remove_const<GridImp>::type::dimension}; 455 456 typedef typename std::pair<const typename UG_NS<dim>::Element*, int> Face; 457 458 /** \brief Look for copy of a face on the next-lower grid level. 459 460 \todo This method is not implemented very efficiently 461 */ getFatherFace(const Face & face)462 static Face getFatherFace(const Face& face) { 463 464 // set up result object 465 Face resultFace; 466 resultFace.first = UG_NS<dim>::EFather(face.first); 467 468 // If there is no father element then we know there is no father face 469 /** \bug This is not true when doing vertical load balancing. */ 470 if (resultFace.first == nullptr) 471 return resultFace; 472 473 // Get all corners of the face 474 std::set<const typename UG_NS<dim>::Vertex*> corners; 475 476 for (int i=0; i<UG_NS<dim>::Corners_Of_Side(face.first, face.second); i++) 477 corners.insert(UG_NS<dim>::Corner(face.first, UG_NS<dim>::Corner_Of_Side(face.first, face.second, i))->myvertex); 478 479 // Loop over all faces of the father element and look for a face that has the same vertices 480 for (int i=0; i<UG_NS<dim>::Sides_Of_Elem(resultFace.first); i++) { 481 482 // Copy father face into set 483 std::set<const typename UG_NS<dim>::Vertex*> fatherFaceCorners; 484 485 for (int j=0; j<UG_NS<dim>::Corners_Of_Side(resultFace.first, i); j++) 486 fatherFaceCorners.insert(UG_NS<dim>::Corner(resultFace.first, UG_NS<dim>::Corner_Of_Side(resultFace.first, i, j))->myvertex); 487 488 // Do the vertex sets match? 489 if (corners==fatherFaceCorners) { 490 resultFace.second = i; 491 return resultFace; 492 } 493 494 } 495 496 // No father face found 497 resultFace.first = nullptr; 498 return resultFace; 499 } 500 501 public: 502 //! constructor stores reference to a grid UGGridIdSet(const GridImp & g)503 UGGridIdSet (const GridImp& g) : grid_(g) {} 504 505 /** \brief Get id of an entity 506 507 We use the remove_const to extract the Type from the mutable class, 508 because the const class is not instantiated yet. 509 510 \bug Since copies of different entities on different levels are supposed to have the 511 same id, we look for the ancestor on the coarsest level that is still a copy of 512 the entity we are interested in. However, the current implementation only searches 513 on one processor, while with UG's vertical load balancing the ancestors of an entity 514 may be distributed across different processors. This will lead to very-difficult-to-fix 515 bugs. Unfortunately, the proper fix for this is not easy, either. 516 */ 517 template<int cd> id(const typename std::remove_const<GridImp>::type::Traits::template Codim<cd>::Entity & e) const518 typename UG_NS<dim>::UG_ID_TYPE id (const typename std::remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const 519 { 520 if (cd==0) { 521 // If we're asked for the id of an element, and that element is a copy of its father, then 522 // we return the id of the lowest ancestor that the element is a copy from. That way copies 523 // of elements have the same id 524 const typename UG_NS<dim>::Element* ancestor = (typename UG_NS<dim>::Element* const)(e.impl().getTarget()); 525 /** \todo We should really be using an isCopy() method rather than hasCopy() */ 526 while (UG_NS<dim>::EFather(ancestor) && UG_NS<dim>::hasCopy(UG_NS<dim>::EFather(ancestor))) 527 ancestor = UG_NS<dim>::EFather(ancestor); 528 529 return UG_NS<dim>::id(ancestor); 530 } 531 532 if (dim-cd==1) { 533 534 const typename UG_NS<dim>::Edge* edge = (typename UG_NS<dim>::Edge* const)(e.impl().getTarget()); 535 536 // If this edge is the copy of an edge on a lower level we return the id of that lower 537 // edge, because Dune wants entities which are copies of each other to have the same id. 538 // BUG: in the parallel setting, we only search on our own processor, but the lowest 539 // copy may actually be on a different processor! 540 const typename UG_NS<dim>::Edge* fatherEdge; 541 fatherEdge = GetFatherEdge(edge); 542 543 while (fatherEdge // fatherEdge exists 544 // ... and it must be a true copy father 545 && ( (fatherEdge->links[0].nbnode->myvertex == edge->links[0].nbnode->myvertex 546 && fatherEdge->links[1].nbnode->myvertex == edge->links[1].nbnode->myvertex) 547 || 548 (fatherEdge->links[0].nbnode->myvertex == edge->links[1].nbnode->myvertex 549 && fatherEdge->links[1].nbnode->myvertex == edge->links[0].nbnode->myvertex) ) ) { 550 edge = fatherEdge; 551 fatherEdge = GetFatherEdge(edge); 552 } 553 554 return UG_NS<dim>::id(edge); 555 } 556 557 558 if (cd == dim) { 559 typename UG_NS<dim>::Node *node = 560 reinterpret_cast<typename UG_NS<dim>::Node *>(e.impl().getTarget()); 561 562 return UG_NS<dim>::id(node); 563 } 564 565 // The entity must be a facet in a 3d grid 566 assert(cd==1 && dim==3); 567 568 typename UG_NS<dim>::Vector *facet = 569 reinterpret_cast<typename UG_NS<dim>::Vector *>(e.impl().getTarget()); 570 571 // If 'facet' is the copy of a facet on a lower level we return the id of that lower-level 572 // facet, because Dune wants entities that are copies of each other to have the same id. 573 // BUG: in a distributed setting, we only search on our own process, but the lowest 574 // copy may actually be on a different process! 575 576 // Going up and down the refinement hierarchy can only be done through an 577 // element that has 'facet' as a face. 578 // TODO: Simplify this if ever SideVector objects get hierarchical information 579 typename UG_NS<dim>::Element* supportingElement = (typename UG_NS<dim>::Element*)facet->object; 580 int thisSideVectorUGNumbering = -1; // Dummy value 581 for (int side=0; side<UG_NS<dim>::Sides_Of_Elem(supportingElement); side++) 582 { 583 auto sideVector = UG_NS<dim>::SideVector(supportingElement,side); 584 if (sideVector==facet) 585 thisSideVectorUGNumbering = side; 586 } 587 588 Face face(supportingElement, thisSideVectorUGNumbering); 589 590 // Go down in the grid hierarchy as long as there is a copy father-facet 591 Face fatherFace; 592 fatherFace = getFatherFace(face); 593 while (fatherFace.first) { 594 face = fatherFace; 595 fatherFace = getFatherFace(face); 596 } 597 598 // Actually return the facet id 599 return UG_NS<dim>::id(UG_NS<dim>::SideVector(face.first, face.second)); 600 } 601 602 /** \brief Get id of subentity 603 604 We use the remove_const to extract the Type from the mutable class, 605 because the const class is not instantiated yet. 606 607 \bug Since copies of different entities on different levels are supposed to have the 608 same id, we look for the ancestor on the coarsest level that is still a copy of 609 the entity we are interested in. However, the current implementation only searches 610 on one processor, while with UG's vertical load balancing the ancestors of an entity 611 may be distributed across different processors. This will lead to very-difficult-to-fix 612 bugs. Unfortunately, the proper fix for this is not easy, either. 613 */ subId(const typename std::remove_const<GridImp>::type::Traits::template Codim<0>::Entity & e,int i,unsigned int codim) const614 typename UG_NS<dim>::UG_ID_TYPE subId (const typename std::remove_const<GridImp>::type::Traits::template Codim<0>::Entity& e, 615 int i, 616 unsigned int codim) const 617 { 618 if (codim==0) 619 return id<0>(e); 620 621 const typename UG_NS<dim>::Element* target = e.impl().getTarget(); 622 GeometryType type = e.type(); 623 624 if (dim-codim==1) { 625 626 auto ref_el = referenceElement<double,dim>(type); 627 auto a = ref_el.subEntity(i,dim-1,0,dim); 628 auto b = ref_el.subEntity(i,dim-1,1,dim); 629 const typename UG_NS<dim>::Edge* edge = UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(a,type)), 630 UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(b,type))); 631 632 // If this edge is the copy of an edge on a lower level we return the id of that lower 633 // edge, because Dune wants entities which are copies of each other to have the same id. 634 // BUG: in the parallel setting, we only search on our own processor, but the lowest 635 // copy may actually be on a different processor! 636 const typename UG_NS<dim>::Edge* fatherEdge; 637 fatherEdge = GetFatherEdge(edge); 638 639 while (fatherEdge // fatherEdge exists 640 // ... and it must be a true copy father 641 && ( (fatherEdge->links[0].nbnode->myvertex == edge->links[0].nbnode->myvertex 642 && fatherEdge->links[1].nbnode->myvertex == edge->links[1].nbnode->myvertex) 643 || 644 (fatherEdge->links[0].nbnode->myvertex == edge->links[1].nbnode->myvertex 645 && fatherEdge->links[1].nbnode->myvertex == edge->links[0].nbnode->myvertex) ) ) { 646 edge = fatherEdge; 647 fatherEdge = GetFatherEdge(edge); 648 } 649 650 return UG_NS<dim>::id(edge); 651 } 652 653 if (codim==1) { // Faces 654 655 Face face(target, UGGridRenumberer<dim>::facesDUNEtoUG(i,type)); 656 657 // If this face is the copy of a face on a lower level we return the id of that lower 658 // face, because Dune wants entities which are copies of each other to have the same id. 659 // BUG: in the parallel setting, we only search on our own processor, but the lowest 660 // copy may actually be on a different processor! 661 Face fatherFace; 662 fatherFace = getFatherFace(face); 663 while (fatherFace.first) { 664 face = fatherFace; 665 fatherFace = getFatherFace(face); 666 } 667 668 return UG_NS<dim>::id(UG_NS<dim>::SideVector(face.first, face.second)); 669 } 670 671 if (codim==dim) { 672 return UG_NS<dim>::id(UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(i,type))); 673 } 674 675 DUNE_THROW(GridError, "UGGrid<" << dim << ">::subId isn't implemented for codim==" << codim ); 676 } 677 678 //private: 679 680 const GridImp& grid_; 681 }; 682 683 } // namespace Dune 684 685 #endif 686