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_ONE_D_GRID_INTERSECTIONS_HH 4 #define DUNE_ONE_D_GRID_INTERSECTIONS_HH 5 6 /** \file 7 * \brief The OneDGridLevelIntersection and OneDGridLeafIntersection classes 8 */ 9 10 #include <dune/grid/onedgrid/onedgridentity.hh> 11 12 namespace Dune { 13 14 /** \brief Intersection between two neighboring elements on a level grid */ 15 template<class GridImp> 16 class OneDGridLevelIntersection 17 { 18 enum { dim=GridImp::dimension }; 19 enum { dimworld=GridImp::dimensionworld }; 20 21 // The corresponding iterator needs to access all members 22 friend class OneDGridLevelIntersectionIterator<GridImp>; 23 24 template<typename,typename> 25 friend class Intersection; 26 OneDGridLevelIntersection()27 OneDGridLevelIntersection() 28 : center_(nullptr) 29 , neighbor_(-1) // marker for invalid intersection 30 {} 31 32 //! Constructor for a given grid entity and a given neighbor OneDGridLevelIntersection(OneDEntityImp<1> * center,int nb)33 OneDGridLevelIntersection(OneDEntityImp<1>* center, int nb) 34 : center_(center), neighbor_(nb) 35 {} 36 37 /** \brief Constructor creating the 'one-after-last'-iterator */ OneDGridLevelIntersection(OneDEntityImp<1> * center)38 OneDGridLevelIntersection(OneDEntityImp<1>* center) 39 : center_(center), neighbor_(2) 40 {} 41 42 typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl; 43 typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl; 44 45 public: 46 47 typedef typename GridImp::template Codim<1>::Geometry Geometry; 48 typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry; 49 typedef typename GridImp::template Codim<0>::Entity Entity; 50 51 //! equality equals(const OneDGridLevelIntersection<GridImp> & other) const52 bool equals(const OneDGridLevelIntersection<GridImp>& other) const { 53 return (center_ == other.center_) && (neighbor_ == other.neighbor_); 54 } 55 target() const56 OneDEntityImp<1>* target() const { 57 const bool isValid = center_ && neighbor_>=0 && neighbor_<2; 58 59 if (!isValid) 60 return center_; 61 else if (neighbor_==0) 62 return center_->pred_; 63 else 64 return center_->succ_; 65 66 } 67 68 //! return true if intersection is with boundary. boundary() const69 bool boundary () const { 70 71 // Check whether we're on the left boundary 72 if (neighbor_==0) { 73 74 // If there's an element to the left we can't be on the boundary 75 if (center_->pred_) 76 return false; 77 78 const OneDEntityImp<1>* ancestor = center_; 79 80 while (ancestor->level_!=0) { 81 82 // Check if we're the left son of our father 83 if (ancestor != ancestor->father_->sons_[0]) 84 return false; 85 86 ancestor = ancestor->father_; 87 } 88 89 // We have reached level 0. If there is no element of the left 90 // we're truly on the boundary 91 return !ancestor->pred_; 92 } 93 94 // //////////////////////////////// 95 // Same for the right boundary 96 // //////////////////////////////// 97 // If there's an element to the right we can't be on the boundary 98 if (center_->succ_) 99 return false; 100 101 const OneDEntityImp<1>* ancestor = center_; 102 103 while (ancestor->level_!=0) { 104 105 // Check if we're the left son of our father 106 if (ancestor != ancestor->father_->sons_[1]) 107 return false; 108 109 ancestor = ancestor->father_; 110 } 111 112 // We have reached level 0. If there is no element of the left 113 // we're truly on the boundary 114 return !ancestor->succ_; 115 116 } 117 118 //! return true if across the edge a neighbor on this level exists neighbor() const119 bool neighbor () const { 120 assert(neighbor_ >= 0 && neighbor_ < 2); 121 122 return (neighbor_==0) 123 ? center_->pred_ && center_->pred_->vertex_[1] == center_->vertex_[0] 124 : center_->succ_ && center_->succ_->vertex_[0] == center_->vertex_[1]; 125 126 } 127 128 //! return true if intersection is conform. conforming() const129 bool conforming () const { 130 return true; 131 } 132 133 //! return the Entity on the inside of this intersection 134 //! (that is the Entity where we started this Iterator) inside() const135 Entity inside() const 136 { 137 return Entity(OneDGridEntity<0,dim,GridImp>(center_)); 138 } 139 140 //! return the Entity on the outside of this intersection 141 //! (that is the neighboring Entity) outside() const142 Entity outside() const 143 { 144 assert(neighbor()); 145 return Entity(OneDGridEntity<0,dim,GridImp>(target())); 146 } 147 148 //! return index of the boundary segment boundarySegmentIndex() const149 int boundarySegmentIndex () const { 150 // It is hardwired here that the domain is connected, i.e., the boundary consists of two points 151 return ((neighbor_==0 && center_->reversedBoundarySegmentNumbering_==false) 152 || (neighbor_==1 && center_->reversedBoundarySegmentNumbering_==true)) ? 0 : 1; 153 } 154 155 //! Here returned element is in LOCAL coordinates of the element 156 //! where iteration started. geometryInInside() const157 LocalGeometry geometryInInside () const 158 { 159 return LocalGeometry( LocalGeometryImpl( (indexInInside() == 0) ? 0 : 1 ) ); 160 } 161 162 //! intersection of codimension 1 of this neighbor with element where iteration started. 163 //! Here returned element is in LOCAL coordinates of neighbor geometryInOutside() const164 LocalGeometry geometryInOutside () const 165 { 166 return LocalGeometry( LocalGeometryImpl( (indexInInside() == 0) ? 1 : 0 ) ); 167 } 168 169 //! intersection of codimension 1 of this neighbor with element where iteration started. 170 //! Here returned element is in GLOBAL coordinates of the element where iteration started. geometry() const171 Geometry geometry () const 172 { 173 return Geometry( GeometryImpl( center_->vertex_[neighbor_]->pos_ ) ); 174 } 175 176 /** \brief obtain the type of reference element for this intersection */ type() const177 GeometryType type () const 178 { 179 return GeometryTypes::vertex; 180 } 181 182 //! local index of codim 1 entity in self where intersection is contained in indexInInside() const183 int indexInInside () const 184 { 185 return neighbor_; 186 } 187 188 //! local index of codim 1 entity in neighbor where intersection is contained indexInOutside() const189 int indexInOutside () const 190 { 191 // If numberInSelf is 0 then numberInNeighbor is 1 and vice versa 192 return 1-neighbor_; 193 } 194 195 //! return outer normal outerNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const196 FieldVector<typename GridImp::ctype, dimworld> outerNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local ) const { 197 return centerUnitOuterNormal(); 198 } 199 200 //! Return outer normal scaled with the integration element integrationOuterNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const201 FieldVector<typename GridImp::ctype, dimworld> integrationOuterNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local ) const { 202 return centerUnitOuterNormal(); 203 } 204 205 //! return unit outer normal unitOuterNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const206 FieldVector<typename GridImp::ctype, dimworld> unitOuterNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local ) const { 207 return centerUnitOuterNormal(); 208 } 209 210 //! return unit outer normal at center of intersection centerUnitOuterNormal() const211 FieldVector<typename GridImp::ctype, dimworld> centerUnitOuterNormal () const { 212 return FieldVector<typename GridImp::ctype, dimworld>(2 * neighbor_ - 1); 213 } 214 215 private: 216 //********************************************************** 217 // private methods 218 //********************************************************** 219 220 OneDEntityImp<1>* center_; 221 222 /** \brief Count on which neighbor we are lookin' at. Can be only 0 or 1. */ 223 int neighbor_; 224 225 }; 226 227 228 /** \brief Intersection between two neighboring elements on a leaf grid */ 229 template<class GridImp> 230 class OneDGridLeafIntersection 231 { 232 enum { dim=GridImp::dimension }; 233 enum { dimworld=GridImp::dimensionworld }; 234 235 // The corresponding iterator needs to access all members 236 friend class OneDGridLeafIntersectionIterator<GridImp>; 237 238 template<typename,typename> 239 friend class Intersection; 240 OneDGridLeafIntersection()241 OneDGridLeafIntersection() 242 : center_(nullptr) 243 , neighbor_(-1) // marker for invalid intersection 244 {} 245 246 //! Constructor for a given grid entity and a given neighbor OneDGridLeafIntersection(OneDEntityImp<1> * center,int nb)247 OneDGridLeafIntersection(OneDEntityImp<1>* center, int nb) 248 : center_(center), neighbor_(nb) 249 {} 250 251 /** \brief Constructor creating the 'one-after-last'-iterator */ OneDGridLeafIntersection(OneDEntityImp<1> * center)252 OneDGridLeafIntersection(OneDEntityImp<1>* center) 253 : center_(center), neighbor_(2) 254 {} 255 256 typedef typename GridImp::Traits::template Codim< 1 >::GeometryImpl GeometryImpl; 257 typedef typename GridImp::Traits::template Codim< 1 >::LocalGeometryImpl LocalGeometryImpl; 258 259 public: 260 261 typedef typename GridImp::template Codim<1>::Geometry Geometry; 262 typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry; 263 typedef typename GridImp::template Codim<0>::Entity Entity; 264 265 //! equality equals(const OneDGridLeafIntersection<GridImp> & other) const266 bool equals(const OneDGridLeafIntersection<GridImp>& other) const { 267 return (center_ == other.center_) && (neighbor_ == other.neighbor_); 268 } 269 target() const270 OneDEntityImp<1>* target() const { 271 const bool isValid = center_ && neighbor_>=0 && neighbor_<2; 272 273 if (!isValid) 274 return center_; 275 276 if (neighbor_==0) { 277 278 // Get left leaf neighbor 279 if (center_->pred_ && center_->pred_->vertex_[1] == center_->vertex_[0]) { 280 281 OneDEntityImp<1>* leftLeafNeighbor = center_->pred_; 282 while (!leftLeafNeighbor->isLeaf()) { 283 assert (leftLeafNeighbor->sons_[1] != NULL); 284 leftLeafNeighbor = leftLeafNeighbor->sons_[1]; 285 } 286 return leftLeafNeighbor; 287 288 } else { 289 290 OneDEntityImp<1>* ancestor = center_; 291 while (ancestor->father_) { 292 ancestor = ancestor->father_; 293 if (ancestor->pred_ && ancestor->pred_->vertex_[1] == ancestor->vertex_[0]) { 294 assert(ancestor->pred_->isLeaf()); 295 return ancestor->pred_; 296 } 297 } 298 299 DUNE_THROW(GridError, "Programming error, apparently we're on the left boundary, neighbor_==2 should not occur!"); 300 } 301 302 } else { 303 304 // Get right leaf neighbor 305 if (center_->succ_ && center_->succ_->vertex_[0] == center_->vertex_[1]) { 306 307 OneDEntityImp<1>* rightLeafNeighbor = center_->succ_; 308 while (!rightLeafNeighbor->isLeaf()) { 309 assert (rightLeafNeighbor->sons_[0] != NULL); 310 rightLeafNeighbor = rightLeafNeighbor->sons_[0]; 311 } 312 return rightLeafNeighbor; 313 314 } else { 315 316 OneDEntityImp<1>* ancestor = center_; 317 while (ancestor->father_) { 318 ancestor = ancestor->father_; 319 if (ancestor->succ_ && ancestor->succ_->vertex_[0] == ancestor->vertex_[1]) { 320 assert(ancestor->succ_->isLeaf()); 321 return ancestor->succ_; 322 } 323 } 324 325 DUNE_THROW(GridError, "Programming error, apparently we're on the right boundary, neighbor_==3 should not occur!"); 326 } 327 328 } 329 330 } 331 332 //! return true if intersection is with boundary. boundary() const333 bool boundary () const { 334 335 // Check whether we're on the left boundary 336 if (neighbor_==0) { 337 338 // If there's an element to the left we can't be on the boundary 339 if (center_->pred_) 340 return false; 341 342 const OneDEntityImp<1>* ancestor = center_; 343 344 while (ancestor->level_!=0) { 345 346 // Check if we're the left son of our father 347 if (ancestor != ancestor->father_->sons_[0]) 348 return false; 349 350 ancestor = ancestor->father_; 351 } 352 353 // We have reached level 0. If there is no element of the left 354 // we're truly on the boundary 355 return !ancestor->pred_; 356 } 357 358 // //////////////////////////////// 359 // Same for the right boundary 360 // //////////////////////////////// 361 // If there's an element to the right we can't be on the boundary 362 if (center_->succ_) 363 return false; 364 365 const OneDEntityImp<1>* ancestor = center_; 366 367 while (ancestor->level_!=0) { 368 369 // Check if we're the left son of our father 370 if (ancestor != ancestor->father_->sons_[1]) 371 return false; 372 373 ancestor = ancestor->father_; 374 } 375 376 // We have reached level 0. If there is no element of the left 377 // we're truly on the boundary 378 return !ancestor->succ_; 379 380 } 381 382 //! return true if across the edge an neighbor on this level exists neighbor() const383 bool neighbor () const { 384 return !boundary(); 385 } 386 387 //! return true if intersection is conform. conforming() const388 bool conforming () const { 389 return true; 390 } 391 392 //! return Entity on the inside of this intersection 393 //! (that is the Entity where we started this Iterator) inside() const394 Entity inside() const 395 { 396 return Entity(OneDGridEntity<0,dim,GridImp>(center_)); 397 } 398 399 //! return the Entity on the outside of this intersection 400 //! (that is the neighboring Entity) outside() const401 Entity outside() const 402 { 403 return Entity(OneDGridEntity<0,dim,GridImp>(target())); 404 } 405 406 //! return index of the boundary segment boundarySegmentIndex() const407 int boundarySegmentIndex () const { 408 // It is hardwired here that the domain is connected, i.e., the boundary consists of two points 409 return ((neighbor_==0 && center_->reversedBoundarySegmentNumbering_==false) 410 || (neighbor_==1 && center_->reversedBoundarySegmentNumbering_==true)) ? 0 : 1; 411 } 412 413 //! Here returned element is in LOCAL coordinates of the element 414 //! where iteration started. geometryInInside() const415 LocalGeometry geometryInInside () const 416 { 417 return LocalGeometry( LocalGeometryImpl( (indexInInside() == 0) ? 0 : 1 ) ); 418 } 419 420 //! intersection of codimension 1 of this neighbor with element where iteration started. 421 //! Here returned element is in LOCAL coordinates of neighbor geometryInOutside() const422 LocalGeometry geometryInOutside () const 423 { 424 return LocalGeometry( LocalGeometryImpl( (indexInInside() == 0) ? 1 : 0 ) ); 425 } 426 427 //! intersection of codimension 1 of this neighbor with element where iteration started. 428 //! Here returned element is in GLOBAL coordinates of the element where iteration started. geometry() const429 Geometry geometry () const 430 { 431 return Geometry( GeometryImpl( center_->vertex_[neighbor_%2]->pos_ ) ); 432 } 433 434 /** \brief obtain the type of reference element for this intersection */ type() const435 GeometryType type () const 436 { 437 return GeometryTypes::vertex; 438 } 439 440 //! local index of codim 1 entity in self where intersection is contained in indexInInside() const441 int indexInInside () const 442 { 443 return neighbor_ % 2; 444 } 445 446 //! local index of codim 1 entity in neighbor where intersection is contained indexInOutside() const447 int indexInOutside () const 448 { 449 // If numberInSelf is 0 then numberInNeighbor is 1 and vice versa 450 return 1-(neighbor_ % 2); 451 } 452 453 //! return outer normal outerNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const454 FieldVector<typename GridImp::ctype, dimworld> outerNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local) const { 455 return centerUnitOuterNormal(); 456 } 457 458 //! Return outer normal scaled with the integration element integrationOuterNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const459 FieldVector<typename GridImp::ctype, dimworld> integrationOuterNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local) const 460 { 461 return centerUnitOuterNormal(); 462 } 463 464 //! return unit outer normal unitOuterNormal(const FieldVector<typename GridImp::ctype,dim-1> & local) const465 FieldVector<typename GridImp::ctype, dimworld> unitOuterNormal ([[maybe_unused]] const FieldVector<typename GridImp::ctype, dim-1>& local) const { 466 return centerUnitOuterNormal(); 467 } 468 469 //! return unit outer normal at center of intersection centerUnitOuterNormal() const470 FieldVector<typename GridImp::ctype, dimworld> centerUnitOuterNormal () const { 471 return FieldVector<typename GridImp::ctype, dimworld>(2 * neighbor_ - 1); 472 } 473 474 private: 475 //********************************************************** 476 // private methods 477 //********************************************************** 478 479 OneDEntityImp<1>* center_; 480 481 /** \brief Count on which neighbor we are lookin' at 482 483 0,1 are the level neighbors, 2 and 3 are the leaf neighbors, 484 if they differ from the level neighbors. */ 485 int neighbor_; 486 487 }; 488 489 } // namespace Dune 490 491 #endif 492