1 #ifndef DUNE_FACEUTILITY_IMP_HH 2 #define DUNE_FACEUTILITY_IMP_HH 3 4 namespace Dune 5 { 6 7 8 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 9 inline ALU3dGridFaceInfo< dim, dimworld, type, Comm >:: ALU3dGridFaceInfo(const bool levelIntersection)10 ALU3dGridFaceInfo( const bool levelIntersection ) : 11 face_(0), 12 innerElement_(0), 13 outerElement_(0), 14 innerFaceNumber_(-1), 15 outerFaceNumber_(-1), 16 innerTwist_(-665), 17 outerTwist_(-665), 18 segmentId_( -1 ), 19 bndId_( -1 ), 20 bndType_( noBoundary ), 21 conformanceState_(UNDEFINED), 22 conformingRefinement_( false ), 23 ghostCellsEnabled_( false ), 24 levelIntersection_( levelIntersection ) 25 { 26 } 27 28 // points face from inner element away? 29 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 30 inline void 31 ALU3dGridFaceInfo< dim, dimworld, type, Comm >:: setFlags(const bool conformingRefinement,const bool ghostCellsEnabled)32 setFlags(const bool conformingRefinement, 33 const bool ghostCellsEnabled ) 34 { 35 conformingRefinement_ = conformingRefinement; 36 ghostCellsEnabled_ = ghostCellsEnabled; 37 } 38 39 // points face from inner element away? 40 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 41 inline void 42 ALU3dGridFaceInfo< dim, dimworld, type, Comm >:: updateFaceInfo(const GEOFaceType & face,int innerLevel,int innerTwist)43 updateFaceInfo(const GEOFaceType& face, 44 int innerLevel, 45 int innerTwist) 46 { 47 face_ = &face; 48 49 innerElement_ = 0; 50 outerElement_ = 0; 51 innerFaceNumber_ = -1; 52 outerFaceNumber_ = -1; 53 bndType_ = noBoundary; 54 segmentId_ = -1; 55 bndId_ = 0; // inner face 56 57 // points face from inner element away? 58 if (innerTwist < 0) 59 { 60 innerElement_ = face.nb.rear().first; 61 innerFaceNumber_ = face.nb.rear().second; 62 outerElement_ = face.nb.front().first; 63 outerFaceNumber_ = face.nb.front().second; 64 } 65 else 66 { 67 innerElement_ = face.nb.front().first; 68 innerFaceNumber_ = face.nb.front().second; 69 outerElement_ = face.nb.rear().first; 70 outerFaceNumber_ = face.nb.rear().second; 71 } // end if 72 73 // if not true we are accessing a fake bnd 74 alugrid_assert ( innerElement_->isRealObject() ); 75 // if not true we are accessing a fake bnd 76 alugrid_assert ( outerElement_->isRealObject() ); 77 78 // we only have to do this in parallel runs 79 if( parallel() && innerElement_->isboundary() ) 80 { 81 bndType_ = innerGhostBoundary; 82 alugrid_assert ( ! dynamic_cast< const GEOPeriodicType* > ( innerElement_ ) ); 83 } 84 85 if( parallel() && innerBoundary() ) 86 { 87 // check for ghosts 88 // this check is only need in the parallel case 89 const BNDFaceType * bnd = static_cast<const BNDFaceType *> (innerElement_); 90 91 if(bnd->bndtype() == ALU3DSPACE ProcessorBoundary_t) 92 { 93 // if nonconformity occurs then go up one level 94 if( bnd->level () != bnd->ghostLevel() && !conformingRefinement_) 95 { 96 bnd = static_cast<const BNDFaceType *>(bnd->up()); 97 alugrid_assert ( bnd ); 98 innerElement_ = static_cast<const HasFaceType*> (bnd); 99 } 100 101 // get ghost and internal number 102 GhostPairType p = bnd->getGhost(); 103 104 // get face number 105 innerFaceNumber_ = p.second; 106 107 // this doesn't count as outer boundary 108 const GEOElementType* ghost = static_cast<const GEOElementType*> (p.first); 109 alugrid_assert (ghost); 110 111 innerTwist_ = ghost->twist(innerFaceNumber_); 112 } 113 else 114 { 115 innerTwist_ = innerFace().twist(innerALUFaceIndex()); 116 } 117 } 118 else 119 { 120 // set inner twist 121 alugrid_assert (innerTwist == innerEntity().twist(innerFaceNumber_)); 122 innerTwist_ = innerTwist; 123 } 124 125 //in the case of a levelIntersectionIterator and conforming elements 126 //we assume the macro grid view. So we go up to level 0 127 //after that we have to get new twist and facenumbers 128 if(levelIntersection_ && conformingRefinement_ && ! (innerElement_->isboundary() ) ) 129 { 130 const GEOElementType * inner = static_cast<const GEOElementType *> (innerElement_); 131 while( inner -> up () ) inner = static_cast<const GEOElementType *> ( inner ->up() ); 132 innerElement_ = static_cast<const HasFaceType *> (inner); 133 innerTwist_ = innerEntity().twist(innerFaceNumber_); 134 } 135 136 if( outerElement_->isboundary() ) 137 { 138 alugrid_assert ( ! innerBoundary() ); 139 // set to default boundary (with domain boundary) 140 bndType_ = domainBoundary ; 141 142 // check for ghosts 143 // this check is only need in the parallel case 144 // if this cast fails we have a periodic element 145 const BNDFaceType * bnd = dynamic_cast<const BNDFaceType *> (outerElement_); 146 const bool periodicBnd = ( bnd == 0 ) ; 147 148 if( periodicBnd ) // the periodic case 149 { 150 bndType_ = periodicBoundary ; 151 alugrid_assert ( dynamic_cast< const GEOPeriodicType* > ( outerElement_ ) ); 152 const GEOPeriodicType* periodicClosure = static_cast< const GEOPeriodicType* > ( outerElement_ ) ; 153 154 // previously, the segmentId( 1 - outerFaceNumber_ ) was used, why? 155 // compute segment already since it's complicated to obtain 156 segmentId_ = periodicClosure->segmentId( outerFaceNumber_ ); 157 bndId_ = periodicClosure->bndtype( outerFaceNumber_ ); 158 159 const GEOFaceType* face = ImplTraits::getFace( *periodicClosure, 1 - outerFaceNumber_ ); 160 alugrid_assert ( (face->nb.front().first == periodicClosure) || (face->nb.rear().first == periodicClosure) ); 161 if( face->nb.rear().first == periodicClosure ) 162 { 163 alugrid_assert ( dynamic_cast< const GEOPeriodicType * >( face->nb.rear().first ) ); 164 outerElement_ = face->nb.front().first ; 165 outerFaceNumber_ = face->nb.front().second ; 166 } 167 else 168 { 169 alugrid_assert ( dynamic_cast< const GEOPeriodicType * >( face->nb.front().first ) ); 170 outerElement_ = face->nb.rear().first ; 171 outerFaceNumber_ = face->nb.rear().second ; 172 } 173 174 alugrid_assert ( outerElement_->isRealObject() ); 175 if( outerElement_->isboundary() ) 176 { 177 alugrid_assert ( dynamic_cast< const BNDFaceType * >( outerElement_ ) ); 178 bnd = static_cast< const BNDFaceType * >( outerElement_ ); 179 } 180 else 181 outerTwist_ = outerEntity().twist( outerFaceNumber_ ); 182 } 183 if ( bnd ) // the boundary case 184 { 185 alugrid_assert ( bnd ); 186 187 // if this cast is valid we have either 188 // a boundary or a ghost element 189 // the ghost element case 190 if( parallel() && bnd->bndtype() == ALU3DSPACE ProcessorBoundary_t) 191 { 192 // if nonconformity occurs then go up one level 193 if( bnd->level () != bnd->ghostLevel() && !conformingRefinement_) 194 { 195 bnd = static_cast<const BNDFaceType *>(bnd->up()); 196 alugrid_assert ( bnd ); 197 outerElement_ = static_cast<const HasFaceType*> (bnd); 198 } 199 200 // set boundary type to ghost boundary 201 bndType_ = outerGhostBoundary ; 202 203 if(conformingRefinement_) 204 outerTwist_ = boundaryFace().twist(outerALUFaceIndex()); 205 206 // access ghost only when ghost cells are enabled 207 if( ghostCellsEnabled_ ) 208 { 209 // get ghost and internal number 210 GhostPairType p = bnd->getGhost(); 211 outerFaceNumber_ = p.second; 212 213 const GEOElementType* ghost = static_cast<const GEOElementType*> (p.first); 214 alugrid_assert ( ghost ); 215 outerTwist_ = ghost->twist(outerFaceNumber_); 216 } 217 } 218 else // the normal boundary case 219 { 220 // get outer twist 221 outerTwist_ = boundaryFace().twist(outerALUFaceIndex()); 222 // compute segment index when needed 223 segmentId_ = boundaryFace().segmentId(); 224 bndId_ = boundaryFace().bndtype(); 225 } 226 } 227 } // if outerElement_->isboundary 228 else 229 { 230 // get outer twist 231 outerTwist_ = outerEntity().twist(outerALUFaceIndex()); 232 } 233 234 //in the case of a levelIntersectionIterator and conforming elements 235 //we assume the macro grid view. So we go up to level 0 236 //after that we have to get new twist and facenumbers 237 if(levelIntersection_ && conformingRefinement_ && ! (outerElement_->isboundary() ) ) 238 { 239 const GEOElementType * outer = static_cast<const GEOElementType *> (outerElement_); 240 while( outer -> up () ) outer = static_cast<const GEOElementType *> ( outer ->up() ); 241 outerElement_ = static_cast<const HasFaceType *> (outer); 242 outerTwist_ = outerEntity().twist(outerFaceNumber_); 243 } 244 245 // make sure we got boundary id correctly 246 alugrid_assert ( bndType_ == periodicBoundary || bndType_ == domainBoundary ? bndId_ > 0 : bndId_ == 0 ); 247 248 //make sure twists are set 249 alugrid_assert( innerTwist_ != -665); 250 alugrid_assert( outerTwist_ != -665); 251 252 // set conformance information 253 conformanceState_ = getConformanceState(innerLevel); 254 } 255 256 // points face from inner element away? 257 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 258 inline ALU3dGridFaceInfo< dim, dimworld, type, Comm >:: ALU3dGridFaceInfo(const GEOFaceType & face,int innerTwist)259 ALU3dGridFaceInfo(const GEOFaceType& face, 260 int innerTwist) 261 { 262 updateFaceInfo(face,innerTwist); 263 } 264 265 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > ~ALU3dGridFaceInfo()266 inline ALU3dGridFaceInfo< dim, dimworld, type, Comm >::~ALU3dGridFaceInfo() {} 267 268 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 269 ALU3dGridFaceInfo< dim, dimworld, type, Comm >:: ALU3dGridFaceInfo(const ALU3dGridFaceInfo & orig)270 ALU3dGridFaceInfo ( const ALU3dGridFaceInfo &orig ) 271 : face_(orig.face_), 272 innerElement_(orig.innerElement_), 273 outerElement_(orig.outerElement_), 274 innerFaceNumber_(orig.innerFaceNumber_), 275 outerFaceNumber_(orig.outerFaceNumber_), 276 innerTwist_(orig.innerTwist_), 277 outerTwist_(orig.outerTwist_), 278 segmentId_( orig.segmentId_ ), 279 bndId_( orig.bndId_ ), 280 bndType_( orig.bndType_ ), 281 conformanceState_(orig.conformanceState_), 282 conformingRefinement_( orig.conformingRefinement_ ), 283 ghostCellsEnabled_( orig.ghostCellsEnabled_ ), 284 levelIntersection_( orig.levelIntersection_ ) 285 {} 286 287 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > isElementLike() const288 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::isElementLike() const { 289 return bndType_ < domainBoundary; 290 } 291 292 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > innerBoundary() const293 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::innerBoundary() const { 294 return bndType_ == innerGhostBoundary; 295 } 296 297 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > outerBoundary() const298 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::outerBoundary() const { 299 return bndType_ == domainBoundary; 300 } 301 302 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > boundary() const303 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::boundary() const { 304 return outerBoundary() || (bndType_ == periodicBoundary); 305 } 306 307 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > neighbor() const308 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::neighbor() const 309 { 310 return isElementLike() || ( ghostBoundary() && ghostCellsEnabled_ ); 311 } 312 313 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > ghostBoundary() const314 inline bool ALU3dGridFaceInfo< dim, dimworld, type, Comm >::ghostBoundary () const 315 { 316 // when communicator is No_Comm there is no ghost boundary 317 return parallel() ? ( bndType_ == outerGhostBoundary ) : false ; 318 } 319 320 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 321 inline const typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::GEOFaceType& face() const322 ALU3dGridFaceInfo< dim, dimworld, type, Comm >::face() const 323 { 324 return *face_; 325 } 326 327 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 328 inline const typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::GEOElementType& innerEntity() const329 ALU3dGridFaceInfo< dim, dimworld, type, Comm >::innerEntity() const 330 { 331 alugrid_assert ( ! innerElement_->isboundary() ); 332 return static_cast<const GEOElementType&>(*innerElement_); 333 } 334 335 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 336 inline const typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::GEOElementType& outerEntity() const337 ALU3dGridFaceInfo< dim, dimworld, type, Comm >::outerEntity() const 338 { 339 alugrid_assert ( isElementLike() ); 340 return static_cast<const GEOElementType&>(*outerElement_); 341 } 342 343 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 344 inline const typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::BNDFaceType& innerFace() const345 ALU3dGridFaceInfo< dim, dimworld, type, Comm >::innerFace() const 346 { 347 alugrid_assert ( innerElement_->isboundary() ); 348 return static_cast<const BNDFaceType&>(*innerElement_); 349 } 350 351 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 352 inline const typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::BNDFaceType& boundaryFace() const353 ALU3dGridFaceInfo< dim, dimworld, type, Comm >::boundaryFace() const { 354 alugrid_assert ( ! isElementLike() ); 355 return static_cast<const BNDFaceType&>(*outerElement_); 356 } 357 358 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > outsideLevel() const359 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::outsideLevel() const 360 { 361 alugrid_assert ( outerElement_ ); 362 alugrid_assert ( !isElementLike() || outerEntity().level() == outerElement_->nbLevel() ); 363 alugrid_assert ( isElementLike() || boundaryFace().level() == outerElement_->nbLevel() ); 364 return outerElement_->nbLevel(); 365 } 366 367 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > segmentId() const368 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::segmentId() const 369 { 370 alugrid_assert ( segmentId_ >= 0 ); 371 return segmentId_; 372 } 373 374 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > boundaryId() const375 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::boundaryId() const 376 { 377 return bndId_; 378 } 379 380 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > innerTwist() const381 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::innerTwist() const 382 { 383 // don't check ghost boundaries here 384 alugrid_assert ( ( ! innerBoundary() ) ? 385 innerEntity().twist(innerALUFaceIndex()) == innerTwist_ : true ); 386 return innerTwist_; 387 } 388 389 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > duneTwist(const int faceIdx,const int aluTwist) const390 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::duneTwist(const int faceIdx, const int aluTwist) const 391 { 392 typedef ElementTopologyMapping<type> ElementTopo; 393 typedef FaceTopologyMapping<type> FaceTopo; 394 395 const int mappedZero = 396 FaceTopo::twist(ElementTopo::dune2aluFaceVertex( faceIdx, 0), aluTwist); 397 398 const int twist = 399 (ElementTopo::faceOrientation( faceIdx ) * sign(aluTwist) < 0 ? 400 mappedZero : -mappedZero-1); 401 // see topology.* files for aluTwistMap 402 if( dim == 2 ) 403 { 404 // in 2d twists are either 0 or 1, but because 405 // of the underlying 3d alu grid they could be different 406 // therefore we adjust to the right range 407 const int duneTwst = FaceTopo :: aluTwistMap( twist ); 408 return (duneTwst == 0) ? 0 : 1; 409 } 410 else 411 { 412 return FaceTopo :: aluTwistMap( twist ); 413 } 414 } 415 416 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > outerTwist() const417 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::outerTwist() const 418 { 419 // don't check ghost boundaries here 420 //alugrid_assert ( (outerBoundary_) ? 421 // (outerTwist_ == boundaryFace().twist(0)) : 422 // (! ghostBoundary_) ? 423 // (outerTwist_ == outerEntity().twist(outerALUFaceIndex())) : true 424 // ); 425 return outerTwist_; 426 } 427 428 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > innerALUFaceIndex() const429 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::innerALUFaceIndex() const { 430 return innerFaceNumber_; 431 } 432 433 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > outerALUFaceIndex() const434 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::outerALUFaceIndex() const { 435 return outerFaceNumber_; 436 } 437 438 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 439 typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::ConformanceState conformanceState() const440 inline ALU3dGridFaceInfo< dim, dimworld, type, Comm >::conformanceState() const 441 { 442 alugrid_assert ( conformanceState_ != UNDEFINED ); 443 return conformanceState_; 444 } 445 446 // calculate conformance state 447 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 448 typename ALU3dGridFaceInfo< dim, dimworld, type, Comm >::ConformanceState getConformanceState(const int innerLevel) const449 inline ALU3dGridFaceInfo< dim, dimworld, type, Comm >::getConformanceState(const int innerLevel) const 450 { 451 ConformanceState result = CONFORMING; 452 453 // in case of non-conforming refinement check level difference 454 if( ! conformingRefinement_ ) 455 { 456 // A boundary is always unrefined 457 int levelDifference = 0 ; 458 if ( isElementLike() ) 459 levelDifference = innerLevel - outerEntity().level(); 460 else 461 levelDifference = innerLevel - boundaryFace().level(); 462 463 if (levelDifference < 0) { 464 result = REFINED_OUTER; 465 } 466 else if (levelDifference > 0) { 467 result = REFINED_INNER; 468 } 469 } 470 471 return result; 472 } 473 474 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 475 inline ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >:: ALU3dGridGeometricFaceInfoBase(const ConnectorType & connector)476 ALU3dGridGeometricFaceInfoBase(const ConnectorType& connector) : 477 connector_(connector), 478 coordsSelfLocal_(-1.0), 479 coordsNeighborLocal_(-1.0), 480 generatedGlobal_(false), 481 generatedLocal_(false) 482 { 483 const auto& refFace = (type == tetra) ? 484 Dune::ReferenceElements< alu3d_ctype, 2 >::simplex() : 485 Dune::ReferenceElements< alu3d_ctype, 2 >::cube() ; 486 487 // do the mappings 488 const int numCorners = refFace.size( 2 ); 489 alugrid_assert( numCorners == int(childLocal_.size()) ); 490 for( int i = 0; i < numCorners; ++i ) 491 { 492 childLocal_[ i ] = refFace.position( i, 2 ); 493 } 494 } 495 496 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 497 inline void 498 ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >:: resetFaceGeom()499 resetFaceGeom() 500 { 501 generatedGlobal_ = false; 502 generatedLocal_ = false; 503 } 504 505 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 506 inline const typename ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::CoordinateType& intersectionSelfLocal() const507 ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::intersectionSelfLocal() const { 508 generateLocalGeometries(); 509 alugrid_assert (generatedLocal_); 510 return coordsSelfLocal_; 511 } 512 513 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 514 inline const typename ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::CoordinateType& intersectionNeighborLocal() const515 ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::intersectionNeighborLocal() const { 516 alugrid_assert (!connector_.outerBoundary()); 517 generateLocalGeometries(); 518 alugrid_assert (generatedLocal_); 519 return coordsNeighborLocal_; 520 } 521 522 523 //sepcialisation for tetra and hexa 524 template< int dim, int dimworld, class Comm > 525 inline ALU3dGridGeometricFaceInfoTetra< dim, dimworld, Comm >:: ALU3dGridGeometricFaceInfoTetra(const ConnectorType & connector)526 ALU3dGridGeometricFaceInfoTetra(const ConnectorType& connector) 527 : Base( connector ), normalUp2Date_( false ) 528 {} 529 530 template< int dim, int dimworld, class Comm > 531 inline void ALU3dGridGeometricFaceInfoTetra< dim, dimworld, Comm >:: resetFaceGeom()532 resetFaceGeom() 533 { 534 Base::resetFaceGeom(); 535 normalUp2Date_ = false; 536 } 537 538 template< int dim, int dimworld, class Comm > 539 inline ALU3dGridGeometricFaceInfoTetra< dim, dimworld, Comm >:: ALU3dGridGeometricFaceInfoTetra(const ALU3dGridGeometricFaceInfoTetra & orig)540 ALU3dGridGeometricFaceInfoTetra(const ALU3dGridGeometricFaceInfoTetra& orig) 541 : Base( orig ), normalUp2Date_( orig.normalUp2Date_ ) 542 {} 543 544 template< int dim, int dimworld, class Comm > 545 template <class GeometryImp> 546 inline void 547 ALU3dGridGeometricFaceInfoTetra< dim, dimworld, Comm >:: buildGlobalGeom(GeometryImp & geo) const548 buildGlobalGeom(GeometryImp& geo) const 549 { 550 if (! this->generatedGlobal_) 551 { 552 // calculate the normal 553 const GEOFaceType & face = this->connector_.face(); 554 555 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() , 556 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() , 557 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() ); 558 559 this->generatedGlobal_ = true ; 560 } 561 } 562 563 template< int dim, int dimworld, class Comm > 564 inline FieldVector<alu3d_ctype, 3> & 565 ALU3dGridGeometricFaceInfoTetra< dim, dimworld, Comm >:: outerNormal(const FieldVector<alu3d_ctype,2> & local) const566 outerNormal(const FieldVector<alu3d_ctype, 2>& local) const 567 { 568 // if geomInfo was not reseted then normal is still correct 569 if(!normalUp2Date_) 570 { 571 // calculate the normal 572 const GEOFaceType & face = this->connector_.face(); 573 const alu3d_ctype (&_p0)[3] = face.myvertex(0)->Point(); 574 const alu3d_ctype (&_p1)[3] = face.myvertex(1)->Point(); 575 const alu3d_ctype (&_p2)[3] = face.myvertex(2)->Point(); 576 577 // change sign if face normal points into inner element 578 // factor is 1.0 to get integration outer normal and not volume outer normal 579 const double factor = (this->connector_.innerTwist() < 0) ? 1.0 : -1.0; 580 581 // see mapp_tetra_3d.h for this piece of code 582 outerNormal_[0] = factor * ((_p1[1]-_p0[1]) *(_p2[2]-_p1[2]) - (_p2[1]-_p1[1]) *(_p1[2]-_p0[2])); 583 outerNormal_[1] = factor * ((_p1[2]-_p0[2]) *(_p2[0]-_p1[0]) - (_p2[2]-_p1[2]) *(_p1[0]-_p0[0])); 584 outerNormal_[2] = factor * ((_p1[0]-_p0[0]) *(_p2[1]-_p1[1]) - (_p2[0]-_p1[0]) *(_p1[1]-_p0[1])); 585 586 normalUp2Date_ = true; 587 } // end if mapp ... 588 589 return outerNormal_; 590 } 591 592 //-sepcialisation for and hexa 593 template< int dim, int dimworld, class Comm > 594 inline ALU3dGridGeometricFaceInfoHexa< dim, dimworld, Comm >:: ALU3dGridGeometricFaceInfoHexa(const ConnectorType & connector)595 ALU3dGridGeometricFaceInfoHexa(const ConnectorType& connector) 596 : Base( connector ) 597 , mappingGlobal_() 598 , mappingGlobalUp2Date_(false) 599 {} 600 601 template< int dim, int dimworld, class Comm > 602 inline void ALU3dGridGeometricFaceInfoHexa< dim, dimworld, Comm >:: resetFaceGeom()603 resetFaceGeom() 604 { 605 Base::resetFaceGeom(); 606 mappingGlobalUp2Date_ = false; 607 } 608 609 template< int dim, int dimworld, class Comm > 610 inline ALU3dGridGeometricFaceInfoHexa< dim, dimworld, Comm >:: ALU3dGridGeometricFaceInfoHexa(const ALU3dGridGeometricFaceInfoHexa & orig)611 ALU3dGridGeometricFaceInfoHexa(const ALU3dGridGeometricFaceInfoHexa& orig) 612 : Base( orig ) 613 , mappingGlobal_(orig.mappingGlobal_) 614 , mappingGlobalUp2Date_(orig.mappingGlobalUp2Date_) 615 {} 616 617 template< int dim, int dimworld, class Comm > 618 template <class GeometryImp> 619 inline void 620 ALU3dGridGeometricFaceInfoHexa< dim, dimworld, Comm >:: buildGlobalGeom(GeometryImp & geo) const621 buildGlobalGeom(GeometryImp& geo) const 622 { 623 if (! this->generatedGlobal_) 624 { 625 // calculate the normal 626 const GEOFaceType & face = this->connector_.face(); 627 628 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() , 629 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() , 630 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() , 631 face.myvertex(FaceTopo::dune2aluVertex(3))->Point() ); 632 this->generatedGlobal_ = true ; 633 } 634 } 635 636 template< int dim, int dimworld, class Comm > 637 inline FieldVector<alu3d_ctype, 3> & 638 ALU3dGridGeometricFaceInfoHexa< dim, dimworld, Comm >:: outerNormal(const FieldVector<alu3d_ctype,2> & local) const639 outerNormal(const FieldVector<alu3d_ctype, 2>& local) const 640 { 641 // if mapping calculated and affine, nothing more to do 642 if ( mappingGlobal_.affine () && mappingGlobalUp2Date_ ) 643 return outerNormal_ ; 644 645 // update surface mapping 646 if(! mappingGlobalUp2Date_ ) 647 { 648 const GEOFaceType & face = connector_.face(); 649 // update mapping to actual face 650 mappingGlobal_.buildMapping( 651 face.myvertex( FaceTopo::dune2aluVertex(0) )->Point(), 652 face.myvertex( FaceTopo::dune2aluVertex(1) )->Point(), 653 face.myvertex( FaceTopo::dune2aluVertex(2) )->Point(), 654 face.myvertex( FaceTopo::dune2aluVertex(3) )->Point() 655 ); 656 mappingGlobalUp2Date_ = true; 657 } 658 659 // calculate the normal 660 // has to be calculated every time normal called, because 661 // depends on local 662 if (connector_.innerTwist() < 0) 663 mappingGlobal_.negativeNormal(local,outerNormal_); 664 else 665 mappingGlobal_.normal(local,outerNormal_); 666 667 // end if 668 return outerNormal_; 669 } 670 671 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 672 inline void ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >:: generateLocalGeometries() const673 generateLocalGeometries() const 674 { 675 if (!generatedLocal_) 676 { 677 // Get the coordinates of the face in the reference element of the 678 // adjoining inner and outer elements and initialise the respective 679 // geometries 680 switch (connector_.conformanceState()) 681 { 682 case (ConnectorType::CONFORMING) : 683 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_); 684 // generate outer local geometry only when not at boundary 685 // * in the parallel case, this needs to be altered for the ghost cells 686 if (!connector_.outerBoundary()) { 687 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_); 688 } // end if 689 break; 690 case (ConnectorType::REFINED_INNER) : 691 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_); 692 referenceElementCoordinatesUnrefined(OUTER, coordsNeighborLocal_); 693 break; 694 case (ConnectorType::REFINED_OUTER) : 695 referenceElementCoordinatesUnrefined(INNER, coordsSelfLocal_); 696 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_); 697 break; 698 default : 699 std::cerr << "ERROR: Wrong conformanceState in generateLocalGeometries! in: " << __FILE__ << " line: " << __LINE__<< std::endl; 700 alugrid_assert (false); 701 exit(1); 702 } // end switch 703 704 generatedLocal_ = true; 705 } // end if 706 } 707 708 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 709 inline int ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >:: globalVertexIndex(const int duneFaceIndex,const int aluFaceTwist,const int duneFaceVertexIndex) const710 globalVertexIndex(const int duneFaceIndex, 711 const int aluFaceTwist, 712 const int duneFaceVertexIndex) const 713 { 714 const int localALUIndex = 715 FaceTopo::dune2aluVertex(duneFaceVertexIndex, 716 aluFaceTwist); 717 718 // get local ALU vertex number on the element's face 719 const int localDuneIndex = ElementTopo:: 720 alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex), 721 localALUIndex); 722 723 return getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 3); 724 } 725 726 727 template< int dim, int dimworld, ALU3dGridElementType type, class Comm > 728 inline void ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >:: referenceElementCoordinatesRefined(SideIdentifier side,CoordinateType & result) const729 referenceElementCoordinatesRefined(SideIdentifier side, 730 CoordinateType& result) const 731 { 732 // this is a dune face index 733 const int faceIndex = 734 (side == INNER ? 735 ElementTopo::alu2duneFace(connector_.innerALUFaceIndex()) : 736 ElementTopo::alu2duneFace(connector_.outerALUFaceIndex())); 737 const int faceTwist = 738 (side == INNER ? 739 connector_.innerTwist() : 740 connector_.outerTwist()); 741 742 const ReferenceElementType& refElem = getReferenceElement(); 743 744 for (int i = 0; i < numVerticesPerFace; ++i) 745 { 746 int duneVertexIndex = globalVertexIndex(faceIndex, faceTwist, i); 747 result[i] = refElem.position(duneVertexIndex, 3); 748 } 749 } 750 751 template< int dimworld, ALU3dGridElementType type, class Comm > 752 inline ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: ALU3dGridGeometricFaceInfoBase(const ConnectorType & connector)753 ALU3dGridGeometricFaceInfoBase(const ConnectorType& connector) : 754 connector_(connector), 755 coordsSelfLocal_(-1.0), 756 coordsNeighborLocal_(-1.0), 757 generatedGlobal_(false), 758 generatedLocal_(false) 759 {} 760 761 template< int dimworld, ALU3dGridElementType type, class Comm > 762 inline void 763 ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: resetFaceGeom()764 resetFaceGeom() 765 { 766 generatedGlobal_ = false; 767 generatedLocal_ = false; 768 } 769 770 template< int dimworld, ALU3dGridElementType type, class Comm > 771 inline ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: ALU3dGridGeometricFaceInfoBase(const ALU3dGridGeometricFaceInfoBase & orig)772 ALU3dGridGeometricFaceInfoBase ( const ALU3dGridGeometricFaceInfoBase &orig ) 773 : connector_(orig.connector_), 774 coordsSelfLocal_(orig.coordsSelfLocal_), 775 coordsNeighborLocal_(orig.coordsNeighborLocal_), 776 generatedGlobal_(orig.generatedGlobal_), 777 generatedLocal_(orig.generatedLocal_) 778 {} 779 780 template< int dimworld, ALU3dGridElementType type, class Comm > 781 inline const typename ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >::LocalCoordinateType& intersectionSelfLocal() const782 ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >::intersectionSelfLocal() const { 783 generateLocalGeometries(); 784 alugrid_assert (generatedLocal_); 785 return coordsSelfLocal_; 786 } 787 788 template< int dimworld, ALU3dGridElementType type, class Comm > 789 inline const typename ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >::LocalCoordinateType& intersectionNeighborLocal() const790 ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >::intersectionNeighborLocal() const { 791 alugrid_assert (!connector_.outerBoundary()); 792 generateLocalGeometries(); 793 alugrid_assert (generatedLocal_); 794 return coordsNeighborLocal_; 795 } 796 797 798 //sepcialisation for tetra and hexa 799 template< int dimworld, class Comm > 800 inline ALU3dGridGeometricFaceInfoTetra< 2, dimworld, Comm >:: ALU3dGridGeometricFaceInfoTetra(const ConnectorType & connector)801 ALU3dGridGeometricFaceInfoTetra(const ConnectorType& connector) 802 : Base( connector ), normalUp2Date_( false ) 803 {} 804 805 template< int dimworld, class Comm > 806 inline void ALU3dGridGeometricFaceInfoTetra< 2, dimworld, Comm >:: resetFaceGeom()807 resetFaceGeom() 808 { 809 Base::resetFaceGeom(); 810 normalUp2Date_ = false; 811 } 812 813 template< int dimworld, class Comm > 814 inline ALU3dGridGeometricFaceInfoTetra< 2, dimworld, Comm >:: ALU3dGridGeometricFaceInfoTetra(const ALU3dGridGeometricFaceInfoTetra & orig)815 ALU3dGridGeometricFaceInfoTetra(const ALU3dGridGeometricFaceInfoTetra& orig) 816 : Base( orig ), normalUp2Date_( orig.normalUp2Date_ ) 817 {} 818 819 template< int dimworld, class Comm > 820 template <class GeometryImp> 821 inline void 822 ALU3dGridGeometricFaceInfoTetra< 2, dimworld, Comm >:: buildGlobalGeom(GeometryImp & geo) const823 buildGlobalGeom(GeometryImp& geo) const 824 { 825 //could be wrong in twist sense 826 if (! this->generatedGlobal_) 827 { 828 // calculate the normal 829 const GEOFaceType & face = this->connector_.face(); 830 831 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(1))->Point() , 832 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() ); 833 834 this->generatedGlobal_ = true ; 835 } 836 } 837 838 template< int dimworld, class Comm > 839 inline FieldVector<alu3d_ctype, dimworld> & 840 ALU3dGridGeometricFaceInfoTetra< 2, dimworld, Comm >:: outerNormal(const FieldVector<alu3d_ctype,1> & local) const841 outerNormal(const FieldVector<alu3d_ctype, 1>& local) const 842 { 843 // if geomInfo was not reseted then normal is still correct 844 if(!normalUp2Date_) 845 { 846 // calculate the normal 847 const GEOFaceType & face = this->connector_.face(); 848 const alu3d_ctype (&_p1)[3] = face.myvertex(1)->Point(); 849 const alu3d_ctype (&_p2)[3] = face.myvertex(2)->Point(); 850 851 // change sign if face normal points into inner element 852 // factor is 1.0 to get integration outer normal and not volume outer normal 853 const double factor = (this->connector_.innerTwist() < 0) ? 1.0 : -1.0; 854 855 856 if(dimworld == 2) 857 { 858 // we want the outer normal orhtogonal to the intersection and with length of the intersection 859 outerNormal_[0] = factor * (_p2[1]-_p1[1]); 860 outerNormal_[1] = factor * (_p1[0]-_p2[0]); 861 } 862 //implemented in iterator_imp.cc 863 //else if(dimworld == 3) 864 865 normalUp2Date_ = true; 866 } // end if mapp ... 867 868 return outerNormal_; 869 } 870 871 //-sepcialisation for and hexa 872 template< int dimworld, class Comm > 873 inline ALU3dGridGeometricFaceInfoHexa< 2, dimworld, Comm >:: ALU3dGridGeometricFaceInfoHexa(const ConnectorType & connector)874 ALU3dGridGeometricFaceInfoHexa(const ConnectorType& connector) 875 : Base( connector ) 876 , normalUp2Date_(false) 877 {} 878 879 template< int dimworld, class Comm > 880 inline void ALU3dGridGeometricFaceInfoHexa< 2, dimworld, Comm >:: resetFaceGeom()881 resetFaceGeom() 882 { 883 Base::resetFaceGeom(); 884 normalUp2Date_ = false; 885 } 886 887 template< int dimworld, class Comm > 888 inline ALU3dGridGeometricFaceInfoHexa< 2, dimworld, Comm >:: ALU3dGridGeometricFaceInfoHexa(const ALU3dGridGeometricFaceInfoHexa & orig)889 ALU3dGridGeometricFaceInfoHexa(const ALU3dGridGeometricFaceInfoHexa& orig) 890 : Base( orig ) 891 , normalUp2Date_(orig.normalUp2Date_) 892 {} 893 894 template< int dimworld, class Comm > 895 template <class GeometryImp> 896 inline void 897 ALU3dGridGeometricFaceInfoHexa< 2, dimworld, Comm >:: buildGlobalGeom(GeometryImp & geo) const898 buildGlobalGeom(GeometryImp& geo) const 899 { 900 //could be wrong in twist sense 901 if (! this->generatedGlobal_) 902 { 903 // calculate the normal 904 const GEOFaceType & face = this->connector_.face(); 905 906 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() , 907 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() ); 908 this->generatedGlobal_ = true ; 909 } 910 } 911 912 template< int dimworld, class Comm > 913 inline FieldVector<alu3d_ctype, dimworld> & 914 ALU3dGridGeometricFaceInfoHexa< 2, dimworld, Comm >:: outerNormal(const FieldVector<alu3d_ctype,1> & local) const915 outerNormal(const FieldVector<alu3d_ctype, 1>& local) const 916 { 917 // if geomInfo was not reseted then normal is still correct 918 if(!normalUp2Date_) 919 { 920 // calculate the normal 921 const GEOFaceType & face = this->connector_.face(); 922 const alu3d_ctype (&_p0)[3] = face.myvertex(0)->Point(); 923 const alu3d_ctype (&_p3)[3] = face.myvertex(3)->Point(); 924 925 // change sign if face normal points into inner element 926 // factor is 1.0 to get integration outer normal and not volume outer normal 927 const double factor = (this->connector_.innerTwist() < 0) ? -1.0 : 1.0; 928 929 if(dimworld == 2) 930 { 931 // we want the length of the intersection and orthogonal to it 932 outerNormal_[0] = factor * (_p0[1] - _p3[1]); 933 outerNormal_[1] = factor * (_p3[0] - _p0[0]); 934 } 935 //implemented in iterator_imp.cc 936 //else if(dimworld == 3) 937 938 normalUp2Date_ = true; 939 } // end if mapp ... 940 941 return outerNormal_; 942 } 943 944 template< int dimworld, ALU3dGridElementType type, class Comm > 945 inline void ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: generateLocalGeometries() const946 generateLocalGeometries() const 947 { 948 if (!generatedLocal_) 949 { 950 // Get the coordinates of the face in the reference element of the 951 // adjoining inner and outer elements and initialise the respective 952 // geometries 953 switch (connector_.conformanceState()) 954 { 955 case (ConnectorType::CONFORMING) : 956 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_); 957 // generate outer local geometry only when not at boundary 958 // * in the parallel case, this needs to be altered for the ghost cells 959 if (!connector_.outerBoundary() && !(connector_.conformingRefinement() && connector_.ghostBoundary())) { 960 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_); 961 } // end if 962 break; 963 case (ConnectorType::REFINED_INNER) : 964 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_); 965 referenceElementCoordinatesUnrefined(OUTER, coordsNeighborLocal_); 966 break; 967 case (ConnectorType::REFINED_OUTER) : 968 referenceElementCoordinatesUnrefined(INNER, coordsSelfLocal_); 969 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_); 970 break; 971 default : 972 std::cerr << "ERROR: Wrong conformanceState in generateLocalGeometries! in: " << __FILE__ << " line: " << __LINE__<< std::endl; 973 alugrid_assert (false); 974 exit(1); 975 } // end switch 976 977 generatedLocal_ = true; 978 } // end if 979 } 980 981 template< int dimworld, ALU3dGridElementType type, class Comm > 982 inline int ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: globalVertexIndex(const int duneFaceIndex,const int aluFaceTwist,const int duneFaceVertexIndex) const983 globalVertexIndex(const int duneFaceIndex, 984 const int aluFaceTwist, 985 const int duneFaceVertexIndex) const 986 { 987 //we want vertices 1,2 of the real 3d DUNE face for tetras and 0,1 for hexas 988 const int localALUIndex = 989 FaceTopo::dune2aluVertex(type == tetra ? duneFaceVertexIndex + 1 : duneFaceVertexIndex, 990 aluFaceTwist); 991 992 // get local DUNE vertex number on the element's face - for tetra map 1,2 of real 3d face back to 0,1 by subtracting 1 993 const int localDuneIndex = (type == tetra) ? 994 ElementTopo::alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex), localALUIndex) - 1 995 : 996 ElementTopo::alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex), localALUIndex) 997 ; 998 999 /* std::cout << "duneFaceIndex: " << duneFaceIndex << std::endl; 1000 std::cout << "aluFaceTwist: " << aluFaceTwist << std::endl; 1001 std::cout << "duneFaceVertexIndex: " << duneFaceVertexIndex << std::endl; 1002 std::cout << "localALUIndex: " << localALUIndex << std::endl; 1003 std::cout << "localDuneIndex: " << localDuneIndex << std::endl; 1004 std ::cout << "ReferenceElementindex: " << getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 2) << std::endl << std::endl; */ 1005 assert( localDuneIndex == 0 || localDuneIndex == 1 ); 1006 return getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 2); 1007 } 1008 1009 1010 template< int dimworld, ALU3dGridElementType type, class Comm > 1011 inline void ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >:: referenceElementCoordinatesRefined(SideIdentifier side,LocalCoordinateType & result) const1012 referenceElementCoordinatesRefined(SideIdentifier side, 1013 LocalCoordinateType& result) const 1014 { 1015 // this is a dune face index 1016 const int faceIndex = 1017 (side == INNER ? 1018 ElementTopo::alu2duneFace(connector_.innerALUFaceIndex()) : 1019 ElementTopo::alu2duneFace(connector_.outerALUFaceIndex())); 1020 const int faceTwist = 1021 (side == INNER ? 1022 connector_.innerTwist() : 1023 connector_.outerTwist()); 1024 1025 const ReferenceElementType& refElem = getReferenceElement(); 1026 1027 for (int i = 0; i < numVerticesPerFace; ++i) 1028 { 1029 int duneVertexIndex = globalVertexIndex(faceIndex, faceTwist, i); 1030 result[i] = refElem.position(duneVertexIndex, 2); 1031 } 1032 } 1033 } //end namespace Dune 1034 #endif 1035