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