dune-alugrid 3.0.0
faceutility_imp.cc
Go to the documentation of this file.
1#ifndef DUNE_FACEUTILITY_IMP_HH
2#define DUNE_FACEUTILITY_IMP_HH
3
4namespace Dune
5{
6
7
8 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
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 segmentIndex_( -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
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
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 segmentIndex_ = -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 }
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_);
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());
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());
117 }
118 else
120 // set inner twist
121 alugrid_assert (innerTwist == innerEntity().twist(innerFaceNumber_));
122 innerTwist_ = innerTwist;
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() ) )
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 segmentIndex( 1 - outerFaceNumber_ ) was used, why?
155 segmentIndex_ = periodicClosure->segmentIndex( outerFaceNumber_ );
156 bndId_ = periodicClosure->bndtype( outerFaceNumber_ );
157
158 const GEOFaceType* face = ImplTraits::getFace( *periodicClosure, 1 - outerFaceNumber_ );
159 alugrid_assert ( (face->nb.front().first == periodicClosure) || (face->nb.rear().first == periodicClosure) );
160 if( face->nb.rear().first == periodicClosure )
161 {
162 alugrid_assert ( dynamic_cast< const GEOPeriodicType * >( face->nb.rear().first ) );
163 outerElement_ = face->nb.front().first ;
164 outerFaceNumber_ = face->nb.front().second ;
165 }
166 else
167 {
168 alugrid_assert ( dynamic_cast< const GEOPeriodicType * >( face->nb.front().first ) );
169 outerElement_ = face->nb.rear().first ;
170 outerFaceNumber_ = face->nb.rear().second ;
171 }
172
173 alugrid_assert ( outerElement_->isRealObject() );
174 if( outerElement_->isboundary() )
175 {
176 alugrid_assert ( dynamic_cast< const BNDFaceType * >( outerElement_ ) );
177 bnd = static_cast< const BNDFaceType * >( outerElement_ );
178 }
179 else
180 outerTwist_ = outerEntity().twist( outerFaceNumber_ );
181 }
182 if ( bnd ) // the boundary case
183 {
184 alugrid_assert ( bnd );
185
186 // if this cast is valid we have either
187 // a boundary or a ghost element
188 // the ghost element case
189 if( parallel() && bnd->bndtype() == ALU3DSPACE ProcessorBoundary_t)
190 {
191 // if nonconformity occurs then go up one level
192 if( bnd->level () != bnd->ghostLevel() && !conformingRefinement_)
193 {
194 bnd = static_cast<const BNDFaceType *>(bnd->up());
195 alugrid_assert ( bnd );
196 outerElement_ = static_cast<const HasFaceType*> (bnd);
197 }
198
199 // set boundary type to ghost boundary
200 bndType_ = outerGhostBoundary ;
201
202 if(conformingRefinement_)
203 outerTwist_ = boundaryFace().twist(outerALUFaceIndex());
204
205 // access ghost only when ghost cells are enabled
206 if( ghostCellsEnabled_ )
207 {
208 // get ghost and internal number
209 GhostPairType p = bnd->getGhost();
210 outerFaceNumber_ = p.second;
211
212 const GEOElementType* ghost = static_cast<const GEOElementType*> (p.first);
213 alugrid_assert ( ghost );
214 outerTwist_ = ghost->twist(outerFaceNumber_);
215 }
216 }
217 else // the normal boundary case
218 {
219 // get outer twist
220 outerTwist_ = boundaryFace().twist(outerALUFaceIndex());
221 // store segment index
222 segmentIndex_ = boundaryFace().segmentIndex();
223 bndId_ = boundaryFace().bndtype();
224 }
225 }
226 } // if outerElement_->isboundary
227 else
228 {
229 // get outer twist
230 outerTwist_ = outerEntity().twist(outerALUFaceIndex());
231 }
232
233 //in the case of a levelIntersectionIterator and conforming elements
234 //we assume the macro grid view. So we go up to level 0
235 //after that we have to get new twist and facenumbers
236 if(levelIntersection_ && conformingRefinement_ && ! (outerElement_->isboundary() ) )
237 {
238 const GEOElementType * outer = static_cast<const GEOElementType *> (outerElement_);
239 while( outer -> up () ) outer = static_cast<const GEOElementType *> ( outer ->up() );
240 outerElement_ = static_cast<const HasFaceType *> (outer);
241 outerTwist_ = outerEntity().twist(outerFaceNumber_);
242 }
243
244 // make sure we got boundary id correctly
245 alugrid_assert ( bndType_ == periodicBoundary || bndType_ == domainBoundary ? bndId_ > 0 : bndId_ == 0 );
246
247 //make sure twists are set
248 alugrid_assert( innerTwist_ != -665);
249 alugrid_assert( outerTwist_ != -665);
250
251 // set conformance information
252 conformanceState_ = getConformanceState(innerLevel);
253 }
254
255 // points face from inner element away?
256 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
259 int innerTwist)
260 {
261 updateFaceInfo(face,innerTwist);
262 }
263
264 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
266
267 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
270 : face_(orig.face_),
271 innerElement_(orig.innerElement_),
272 outerElement_(orig.outerElement_),
273 innerFaceNumber_(orig.innerFaceNumber_),
274 outerFaceNumber_(orig.outerFaceNumber_),
275 innerTwist_(orig.innerTwist_),
276 outerTwist_(orig.outerTwist_),
277 segmentIndex_( orig.segmentIndex_ ),
278 bndId_( orig.bndId_ ),
279 bndType_( orig.bndType_ ),
280 conformanceState_(orig.conformanceState_),
281 conformingRefinement_( orig.conformingRefinement_ ),
282 ghostCellsEnabled_( orig.ghostCellsEnabled_ )
283 {}
284
285 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
287 return bndType_ < domainBoundary;
288 }
289
290 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
292 return bndType_ == innerGhostBoundary;
293 }
294
295 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
297 return bndType_ == domainBoundary;
298 }
299
300 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
302 return outerBoundary() || (bndType_ == periodicBoundary);
303 }
304
305 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
307 {
308 return isElementLike() || ( ghostBoundary() && ghostCellsEnabled_ );
309 }
310
311 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
313 {
314 // when communicator is No_Comm there is no ghost boundary
315 return parallel() ? ( bndType_ == outerGhostBoundary ) : false ;
316 }
317
318 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
321 {
322 return *face_;
323 }
324
325 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
328 {
329 alugrid_assert ( ! innerElement_->isboundary() );
330 return static_cast<const GEOElementType&>(*innerElement_);
331 }
332
333 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
336 {
337 alugrid_assert ( isElementLike() );
338 return static_cast<const GEOElementType&>(*outerElement_);
339 }
340
341 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
344 {
345 alugrid_assert ( innerElement_->isboundary() );
346 return static_cast<const BNDFaceType&>(*innerElement_);
347 }
348
349 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
352 alugrid_assert ( ! isElementLike() );
353 return static_cast<const BNDFaceType&>(*outerElement_);
354 }
355
356 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
358 {
359 alugrid_assert ( outerElement_ );
360 alugrid_assert ( !isElementLike() || outerEntity().level() == outerElement_->nbLevel() );
361 alugrid_assert ( isElementLike() || boundaryFace().level() == outerElement_->nbLevel() );
362 return outerElement_->nbLevel();
363 }
364
365 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
367 {
368 alugrid_assert ( segmentIndex_ >= 0 );
369 return segmentIndex_;
370 }
371
372 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
374 {
375 return bndId_;
376 }
377
378 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
380 {
381 // don't check ghost boundaries here
382 alugrid_assert ( ( ! innerBoundary() ) ?
383 innerEntity().twist(innerALUFaceIndex()) == innerTwist_ : true );
384 return innerTwist_;
385 }
386
387 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
388 inline int ALU3dGridFaceInfo< dim, dimworld, type, Comm >::duneTwist(const int faceIdx, const int aluTwist) const
389 {
390 typedef ElementTopologyMapping<type> ElementTopo;
391 typedef FaceTopologyMapping<type> FaceTopo;
392
393 const int mappedZero =
394 FaceTopo::twist(ElementTopo::dune2aluFaceVertex( faceIdx, 0), aluTwist);
395
396 const int twist =
397 (ElementTopo::faceOrientation( faceIdx ) * sign(aluTwist) < 0 ?
398 mappedZero : -mappedZero-1);
399 // see topology.* files for aluTwistMap
400 if( dim == 2 )
401 {
402 // in 2d twists are either 0 or 1, but because
403 // of the underlying 3d alu grid they could be different
404 // therefore we adjust to the right range
405 const int duneTwst = FaceTopo :: aluTwistMap( twist );
406 return (duneTwst == 0) ? 0 : 1;
407 }
408 else
409 {
410 return FaceTopo :: aluTwistMap( twist );
411 }
412 }
413
414 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
416 {
417 // don't check ghost boundaries here
418 //alugrid_assert ( (outerBoundary_) ?
419 // (outerTwist_ == boundaryFace().twist(0)) :
420 // (! ghostBoundary_) ?
421 // (outerTwist_ == outerEntity().twist(outerALUFaceIndex())) : true
422 // );
423 return outerTwist_;
424 }
425
426 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
428 return innerFaceNumber_;
429 }
430
431 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
433 return outerFaceNumber_;
434 }
435
436 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
439 {
440 alugrid_assert ( conformanceState_ != UNDEFINED );
441 return conformanceState_;
442 }
443
444 // calculate conformance state
445 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
448 {
449 ConformanceState result = CONFORMING;
450
451 // in case of non-conforming refinement check level difference
452 if( ! conformingRefinement_ )
453 {
454 // A boundary is always unrefined
455 int levelDifference = 0 ;
456 if ( isElementLike() )
457 levelDifference = innerLevel - outerEntity().level();
458 else
459 levelDifference = innerLevel - boundaryFace().level();
460
461 if (levelDifference < 0) {
462 result = REFINED_OUTER;
463 }
464 else if (levelDifference > 0) {
465 result = REFINED_INNER;
466 }
467 }
468
469 return result;
470 }
471
472 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
475 connector_(connector),
476 coordsSelfLocal_(-1.0),
477 coordsNeighborLocal_(-1.0),
478 generatedGlobal_(false),
479 generatedLocal_(false)
480 {}
481
482 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
483 inline void
486 {
487 generatedGlobal_ = false;
488 generatedLocal_ = false;
489 }
490
491 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
494 : connector_(orig.connector_),
495 coordsSelfLocal_(orig.coordsSelfLocal_),
496 coordsNeighborLocal_(orig.coordsNeighborLocal_),
497 generatedGlobal_(orig.generatedGlobal_),
498 generatedLocal_(orig.generatedLocal_)
499 {}
500
501 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
504 generateLocalGeometries();
505 alugrid_assert (generatedLocal_);
506 return coordsSelfLocal_;
507 }
508
509 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
512 alugrid_assert (!connector_.outerBoundary());
513 generateLocalGeometries();
514 alugrid_assert (generatedLocal_);
515 return coordsNeighborLocal_;
516 }
517
518
519 //sepcialisation for tetra and hexa
520 template< int dim, int dimworld, class Comm >
523 : Base( connector ), normalUp2Date_( false )
524 {}
525
526 template< int dim, int dimworld, class Comm >
529 {
530 Base::resetFaceGeom();
531 normalUp2Date_ = false;
532 }
533
534 template< int dim, int dimworld, class Comm >
537 : Base( orig ), normalUp2Date_( orig.normalUp2Date_ )
538 {}
539
540 template< int dim, int dimworld, class Comm >
541 template <class GeometryImp>
542 inline void
544 buildGlobalGeom(GeometryImp& geo) const
545 {
546 if (! this->generatedGlobal_)
547 {
548 // calculate the normal
549 const GEOFaceType & face = this->connector_.face();
550
551 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() ,
552 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() ,
553 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() );
554
555 this->generatedGlobal_ = true ;
556 }
557 }
558
559 template< int dim, int dimworld, class Comm >
560 inline FieldVector<alu3d_ctype, 3> &
562 outerNormal(const FieldVector<alu3d_ctype, 2>& local) const
563 {
564 // if geomInfo was not reseted then normal is still correct
565 if(!normalUp2Date_)
566 {
567 // calculate the normal
568 const GEOFaceType & face = this->connector_.face();
569 const alu3d_ctype (&_p0)[3] = face.myvertex(0)->Point();
570 const alu3d_ctype (&_p1)[3] = face.myvertex(1)->Point();
571 const alu3d_ctype (&_p2)[3] = face.myvertex(2)->Point();
572
573 // change sign if face normal points into inner element
574 // factor is 1.0 to get integration outer normal and not volume outer normal
575 const double factor = (this->connector_.innerTwist() < 0) ? 1.0 : -1.0;
576
577 // see mapp_tetra_3d.h for this piece of code
578 outerNormal_[0] = factor * ((_p1[1]-_p0[1]) *(_p2[2]-_p1[2]) - (_p2[1]-_p1[1]) *(_p1[2]-_p0[2]));
579 outerNormal_[1] = factor * ((_p1[2]-_p0[2]) *(_p2[0]-_p1[0]) - (_p2[2]-_p1[2]) *(_p1[0]-_p0[0]));
580 outerNormal_[2] = factor * ((_p1[0]-_p0[0]) *(_p2[1]-_p1[1]) - (_p2[0]-_p1[0]) *(_p1[1]-_p0[1]));
581
582 normalUp2Date_ = true;
583 } // end if mapp ...
584
585 return outerNormal_;
586 }
587
588 //-sepcialisation for and hexa
589 template< int dim, int dimworld, class Comm >
592 : Base( connector )
593 , mappingGlobal_()
594 , mappingGlobalUp2Date_(false)
595 {}
596
597 template< int dim, int dimworld, class Comm >
600 {
601 Base::resetFaceGeom();
602 mappingGlobalUp2Date_ = false;
603 }
604
605 template< int dim, int dimworld, class Comm >
608 : Base( orig )
609 , mappingGlobal_(orig.mappingGlobal_)
610 , mappingGlobalUp2Date_(orig.mappingGlobalUp2Date_)
611 {}
612
613 template< int dim, int dimworld, class Comm >
614 template <class GeometryImp>
615 inline void
617 buildGlobalGeom(GeometryImp& geo) const
618 {
619 if (! this->generatedGlobal_)
620 {
621 // calculate the normal
622 const GEOFaceType & face = this->connector_.face();
623
624 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() ,
625 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() ,
626 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() ,
627 face.myvertex(FaceTopo::dune2aluVertex(3))->Point() );
628 this->generatedGlobal_ = true ;
629 }
630 }
631
632 template< int dim, int dimworld, class Comm >
633 inline FieldVector<alu3d_ctype, 3> &
635 outerNormal(const FieldVector<alu3d_ctype, 2>& local) const
636 {
637 // if mapping calculated and affine, nothing more to do
638 if ( mappingGlobal_.affine () && mappingGlobalUp2Date_ )
639 return outerNormal_ ;
640
641 // update surface mapping
642 if(! mappingGlobalUp2Date_ )
643 {
644 const GEOFaceType & face = connector_.face();
645 // update mapping to actual face
646 mappingGlobal_.buildMapping(
647 face.myvertex( FaceTopo::dune2aluVertex(0) )->Point(),
648 face.myvertex( FaceTopo::dune2aluVertex(1) )->Point(),
649 face.myvertex( FaceTopo::dune2aluVertex(2) )->Point(),
650 face.myvertex( FaceTopo::dune2aluVertex(3) )->Point()
651 );
652 mappingGlobalUp2Date_ = true;
653 }
654
655 // calculate the normal
656 // has to be calculated every time normal called, because
657 // depends on local
658 if (connector_.innerTwist() < 0)
659 mappingGlobal_.negativeNormal(local,outerNormal_);
660 else
661 mappingGlobal_.normal(local,outerNormal_);
662
663 // end if
664 return outerNormal_;
665 }
666
667 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
670 {
671 if (!generatedLocal_)
672 {
673 // Get the coordinates of the face in the reference element of the
674 // adjoining inner and outer elements and initialise the respective
675 // geometries
676 switch (connector_.conformanceState())
677 {
678 case (ConnectorType::CONFORMING) :
679 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_);
680 // generate outer local geometry only when not at boundary
681 // * in the parallel case, this needs to be altered for the ghost cells
682 if (!connector_.outerBoundary()) {
683 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_);
684 } // end if
685 break;
686 case (ConnectorType::REFINED_INNER) :
687 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_);
688 referenceElementCoordinatesUnrefined(OUTER, coordsNeighborLocal_);
689 break;
690 case (ConnectorType::REFINED_OUTER) :
691 referenceElementCoordinatesUnrefined(INNER, coordsSelfLocal_);
692 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_);
693 break;
694 default :
695 std::cerr << "ERROR: Wrong conformanceState in generateLocalGeometries! in: " << __FILE__ << " line: " << __LINE__<< std::endl;
696 alugrid_assert (false);
697 exit(1);
698 } // end switch
699
700 generatedLocal_ = true;
701 } // end if
702 }
703
704 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
705 inline int ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::
706 globalVertexIndex(const int duneFaceIndex,
707 const int aluFaceTwist,
708 const int duneFaceVertexIndex) const
709 {
710 const int localALUIndex =
711 FaceTopo::dune2aluVertex(duneFaceVertexIndex,
712 aluFaceTwist);
713
714 // get local ALU vertex number on the element's face
715 const int localDuneIndex = ElementTopo::
716 alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex),
717 localALUIndex);
718
719 return getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 3);
720 }
721
722
723 template< int dim, int dimworld, ALU3dGridElementType type, class Comm >
724 inline void ALU3dGridGeometricFaceInfoBase< dim, dimworld, type, Comm >::
725 referenceElementCoordinatesRefined(SideIdentifier side,
726 CoordinateType& result) const
727 {
728 // this is a dune face index
729 const int faceIndex =
730 (side == INNER ?
731 ElementTopo::alu2duneFace(connector_.innerALUFaceIndex()) :
732 ElementTopo::alu2duneFace(connector_.outerALUFaceIndex()));
733 const int faceTwist =
734 (side == INNER ?
735 connector_.innerTwist() :
736 connector_.outerTwist());
737
738 const ReferenceElementType& refElem = getReferenceElement();
739
740 for (int i = 0; i < numVerticesPerFace; ++i)
741 {
742 int duneVertexIndex = globalVertexIndex(faceIndex, faceTwist, i);
743 result[i] = refElem.position(duneVertexIndex, 3);
744 }
745 }
746
747 template< int dimworld, ALU3dGridElementType type, class Comm >
750 connector_(connector),
751 coordsSelfLocal_(-1.0),
752 coordsNeighborLocal_(-1.0),
753 generatedGlobal_(false),
754 generatedLocal_(false)
755 {}
756
757 template< int dimworld, ALU3dGridElementType type, class Comm >
758 inline void
765
766 template< int dimworld, ALU3dGridElementType type, class Comm >
769 : connector_(orig.connector_),
770 coordsSelfLocal_(orig.coordsSelfLocal_),
771 coordsNeighborLocal_(orig.coordsNeighborLocal_),
772 generatedGlobal_(orig.generatedGlobal_),
773 generatedLocal_(orig.generatedLocal_)
774 {}
775
776 template< int dimworld, ALU3dGridElementType type, class Comm >
783
784 template< int dimworld, ALU3dGridElementType type, class Comm >
792
793
794 //sepcialisation for tetra and hexa
795 template< int dimworld, class Comm >
798 : Base( connector ), normalUp2Date_( false )
799 {}
800
801 template< int dimworld, class Comm >
804 {
806 normalUp2Date_ = false;
807 }
808
809 template< int dimworld, class Comm >
812 : Base( orig ), normalUp2Date_( orig.normalUp2Date_ )
813 {}
814
815 template< int dimworld, class Comm >
816 template <class GeometryImp>
817 inline void
819 buildGlobalGeom(GeometryImp& geo) const
820 {
821 //could be wrong in twist sense
822 if (! this->generatedGlobal_)
823 {
824 // calculate the normal
825 const GEOFaceType & face = this->connector_.face();
826
827 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(1))->Point() ,
828 face.myvertex(FaceTopo::dune2aluVertex(2))->Point() );
829
830 this->generatedGlobal_ = true ;
831 }
832 }
833
834 template< int dimworld, class Comm >
835 inline FieldVector<alu3d_ctype, dimworld> &
837 outerNormal(const FieldVector<alu3d_ctype, 1>& local) const
838 {
839 // if geomInfo was not reseted then normal is still correct
840 if(!normalUp2Date_)
841 {
842 // calculate the normal
843 const GEOFaceType & face = this->connector_.face();
844 const alu3d_ctype (&_p1)[3] = face.myvertex(1)->Point();
845 const alu3d_ctype (&_p2)[3] = face.myvertex(2)->Point();
846
847 // change sign if face normal points into inner element
848 // factor is 1.0 to get integration outer normal and not volume outer normal
849 const double factor = (this->connector_.innerTwist() < 0) ? 1.0 : -1.0;
850
851
852 if(dimworld == 2)
853 {
854 // we want the outer normal orhtogonal to the intersection and with length of the intersection
855 outerNormal_[0] = factor * (_p2[1]-_p1[1]);
856 outerNormal_[1] = factor * (_p1[0]-_p2[0]);
857 }
858 else if(dimworld == 3)
859 {
860 //we want the outer normal orthogonal to the intersection and to the normal of the inner element, with length of the intersection
861 const GEOElementType * innerElement (0);
862 if(this->connector_.innerBoundary())
863 innerElement = static_cast<const GEOElementType * > (this->connector_.innerFace().getGhost().first) ;
864 else
865 innerElement = static_cast<const GEOElementType * > (&(this->connector_.innerEntity() ));
866 //get the vertex that is opposed to the intersection
867 const alu3d_ctype (&_p3)[3] = innerElement->myvertex(this->connector_.innerALUFaceIndex())->Point();
868
869 FieldVector<alu3d_ctype, 3> normal;
870
871 // calculate normal of the element
872 normal[0] = (_p2[1] - _p3[1]) * (_p1[2] - _p2[2]) - (_p2[2] - _p3[2]) * (_p1[1] - _p2[1]) ;
873 normal[1] = (_p2[2] - _p3[2]) * (_p1[0] - _p2[0]) - (_p2[0] - _p3[0]) * (_p1[2] - _p2[2]) ;
874 normal[2] = (_p2[0] - _p3[0]) * (_p1[1] - _p2[1]) - (_p2[1] - _p3[1]) * (_p1[0] - _p2[0]) ;
875
876 //normalize
877 normal *= (1.0/normal.two_norm());
878
879 //calculate cross product of element normal and intersection
880 //we do not need to include the factor, this is encoded in the order of p1, p2
881 outerNormal_[0] = normal[2] * (_p1[1] - _p2[1]) - normal[1] * (_p1[2] - _p2[2]);
882 outerNormal_[1] = normal[0] * (_p1[2] - _p2[2]) - normal[2] * (_p1[0] - _p2[0]);
883 outerNormal_[2] = normal[1] * (_p1[0] - _p2[0]) - normal[0] * (_p1[1] - _p2[1]);
884 }
885
886 normalUp2Date_ = true;
887 } // end if mapp ...
888
889 return outerNormal_;
890 }
891
892 //-sepcialisation for and hexa
893 template< int dimworld, class Comm >
896 : Base( connector )
897 , normalUp2Date_(false)
898 {}
899
900 template< int dimworld, class Comm >
903 {
905 normalUp2Date_ = false;
906 }
907
908 template< int dimworld, class Comm >
911 : Base( orig )
912 , normalUp2Date_(orig.normalUp2Date_)
913 {}
914
915 template< int dimworld, class Comm >
916 template <class GeometryImp>
917 inline void
919 buildGlobalGeom(GeometryImp& geo) const
920 {
921 //could be wrong in twist sense
922 if (! this->generatedGlobal_)
923 {
924 // calculate the normal
925 const GEOFaceType & face = this->connector_.face();
926
927 geo.buildGeom( face.myvertex(FaceTopo::dune2aluVertex(0))->Point() ,
928 face.myvertex(FaceTopo::dune2aluVertex(1))->Point() );
929 this->generatedGlobal_ = true ;
930 }
931 }
932
933 template< int dimworld, class Comm >
934 inline FieldVector<alu3d_ctype, dimworld> &
936 outerNormal(const FieldVector<alu3d_ctype, 1>& local) const
937 {
938 // if geomInfo was not reseted then normal is still correct
939 if(!normalUp2Date_)
940 {
941 // calculate the normal
942 const GEOFaceType & face = this->connector_.face();
943 const alu3d_ctype (&_p0)[3] = face.myvertex(0)->Point();
944 const alu3d_ctype (&_p3)[3] = face.myvertex(3)->Point();
945
946 // change sign if face normal points into inner element
947 // factor is 1.0 to get integration outer normal and not volume outer normal
948 const double factor = (this->connector_.innerTwist() < 0) ? -1.0 : 1.0;
949
950 if(dimworld == 2)
951 {
952 // we want the length of the intersection and orthogonal to it
953 outerNormal_[0] = factor * (_p0[1] - _p3[1]);
954 outerNormal_[1] = factor * (_p3[0] - _p0[0]);
955 }
956 else if(dimworld == 3)
957 {
958 //Teh following is an adaption of
959 //
960 // const ReferenceElement< alu2d_ctype, dim > &refElement =
961 // ReferenceElements< alu2d_ctype, dim >::cube();
962 // typename LocalGeometry::GlobalCoordinate xInside = geometryInInside().global( local );
963 // typename LocalGeometry::GlobalCoordinate refNormal = refElement.integrationOuterNormal( indexInInside() );
964 // inside()->geometry().jacobianInverseTransposed( xInside ).mv( refNormal, outerNormal );
965 // outerNormal *= inside()->geometry().integrationElement( xInside );
966 //
967 // from the calculation of the former 2d code
968 // also it does not work for some reason - so this is implemented in iterator_imp.cc
969
970
971 /*
972 typedef typename ALU3dGrid<2, dimworld, hexa, Comm>::template Codim<1>::LocalGeometry LocalGeometry;
973 typedef Dune :: ALU3dGridGeometry< 2, dimworld, ALU3dGrid<2, dimworld, hexa, Comm> > GeometryImpl;
974 //generate the local geometries - if not already present
975 this->generateLocalGeometries();
976 //get the inside entity
977 const GEOElementType &inner = this->connector_.innerEntity();
978
979 //Get the 2d cube reference Element
980 const ReferenceElement< alu3d_ctype, 2 > &refElement =
981 ReferenceElements< alu3d_ctype, 2 >::cube();
982 //get xInside in the geometryInInside coordinates
983 typename LocalGeometry::GlobalCoordinate xInside = this->intersectionSelfLocal()[1];
984 xInside *= local[0];
985 xInside.axpy(1-local[0] , this->intersectionSelfLocal()[0]);
986
987 // get the normal of the reference element on the current face (DUNE index from the inside)
988 typename LocalGeometry::GlobalCoordinate refNormal = refElement.integrationOuterNormal( ElementTopologyMapping<hexa>::alu2duneFace(this->connector_.innerALUFaceIndex()) );
989
990 //construct the inner geometry
991 GeometryImpl geo ;
992 geo.buildGeom(inner.myvertex(0)->Point(), inner.myvertex(1)->Point(), inner.myvertex(3)->Point(), inner.myvertex(2)->Point());
993 //map the reference Normal back to the inner geometry
994 geo.jacobianInverseTransposed( xInside ).mv( refNormal, outerNormal_ );
995 outerNormal_ *= geo.integrationElement( xInside );
996
997 */
998 }
999
1000 normalUp2Date_ = true;
1001 } // end if mapp ...
1002
1003 return outerNormal_;
1004 }
1005
1006 template< int dimworld, ALU3dGridElementType type, class Comm >
1009 {
1010 if (!generatedLocal_)
1011 {
1012 // Get the coordinates of the face in the reference element of the
1013 // adjoining inner and outer elements and initialise the respective
1014 // geometries
1015 switch (connector_.conformanceState())
1016 {
1018 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_);
1019 // generate outer local geometry only when not at boundary
1020 // * in the parallel case, this needs to be altered for the ghost cells
1022 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_);
1023 } // end if
1024 break;
1026 referenceElementCoordinatesRefined(INNER, coordsSelfLocal_);
1027 referenceElementCoordinatesUnrefined(OUTER, coordsNeighborLocal_);
1028 break;
1030 referenceElementCoordinatesUnrefined(INNER, coordsSelfLocal_);
1031 referenceElementCoordinatesRefined(OUTER, coordsNeighborLocal_);
1032 break;
1033 default :
1034 std::cerr << "ERROR: Wrong conformanceState in generateLocalGeometries! in: " << __FILE__ << " line: " << __LINE__<< std::endl;
1035 alugrid_assert (false);
1036 exit(1);
1037 } // end switch
1038
1039 generatedLocal_ = true;
1040 } // end if
1041 }
1042
1043 template< int dimworld, ALU3dGridElementType type, class Comm >
1045 globalVertexIndex(const int duneFaceIndex,
1046 const int aluFaceTwist,
1047 const int duneFaceVertexIndex) const
1048 {
1049 //we want vertices 1,2 of the real 3d DUNE face for tetras and 0,1 for hexas
1050 const int localALUIndex =
1051 FaceTopo::dune2aluVertex(type == tetra ? duneFaceVertexIndex + 1 : duneFaceVertexIndex,
1052 aluFaceTwist);
1053
1054 // 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
1055 const int localDuneIndex = (type == tetra) ?
1056 ElementTopo::alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex), localALUIndex) - 1
1057 :
1058 ElementTopo::alu2duneFaceVertex(ElementTopo::dune2aluFace(duneFaceIndex), localALUIndex)
1059 ;
1060
1061 /* std::cout << "duneFaceIndex: " << duneFaceIndex << std::endl;
1062 std::cout << "aluFaceTwist: " << aluFaceTwist << std::endl;
1063 std::cout << "duneFaceVertexIndex: " << duneFaceVertexIndex << std::endl;
1064 std::cout << "localALUIndex: " << localALUIndex << std::endl;
1065 std::cout << "localDuneIndex: " << localDuneIndex << std::endl;
1066 std ::cout << "ReferenceElementindex: " << getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 2) << std::endl << std::endl; */
1067 assert( localDuneIndex == 0 || localDuneIndex == 1 );
1068 return getReferenceElement().subEntity(duneFaceIndex, 1, localDuneIndex, 2);
1069 }
1070
1071
1072 template< int dimworld, ALU3dGridElementType type, class Comm >
1073 inline void ALU3dGridGeometricFaceInfoBase< 2, dimworld, type, Comm >::
1074 referenceElementCoordinatesRefined(SideIdentifier side,
1075 LocalCoordinateType& result) const
1076 {
1077 // this is a dune face index
1078 const int faceIndex =
1079 (side == INNER ?
1080 ElementTopo::alu2duneFace(connector_.innerALUFaceIndex()) :
1081 ElementTopo::alu2duneFace(connector_.outerALUFaceIndex()));
1082 const int faceTwist =
1083 (side == INNER ?
1086
1087 const ReferenceElementType& refElem = getReferenceElement();
1088
1089 for (int i = 0; i < numVerticesPerFace; ++i)
1090 {
1091 int duneVertexIndex = globalVertexIndex(faceIndex, faceTwist, i);
1092 result[i] = refElem.position(duneVertexIndex, 2);
1093 }
1094 }
1095} //end namespace Dune
1096#endif
#define ALU3DSPACE
Definition alu3dinclude.hh:24
#define alugrid_assert(EX)
Definition alugrid_assert.hh:20
Definition alu3dinclude.hh:80
@ tetra
Definition topology.hh:12
double alu3d_ctype
Definition alu3dinclude.hh:85
Definition faceutility.hh:43
bool isElementLike() const
returns true if outerEntity casts into a helement
Definition faceutility_imp.cc:286
int outsideLevel() const
Definition faceutility_imp.cc:357
int outerTwist() const
Twist of the face seen from the outer element.
Definition faceutility_imp.cc:415
ImplTraits::GEOElementType GEOElementType
Definition faceutility.hh:51
ImplTraits::GhostPairType GhostPairType
Definition faceutility.hh:54
ConformanceState
Definition faceutility.hh:48
@ CONFORMING
Definition faceutility.hh:48
@ REFINED_INNER
Definition faceutility.hh:48
@ REFINED_OUTER
Definition faceutility.hh:48
~ALU3dGridFaceInfo()
Destructor.
Definition faceutility_imp.cc:265
void setFlags(const bool conformingRefinement, const bool ghostCellsEnabled)
reset flags
Definition faceutility_imp.cc:32
void updateFaceInfo(const GEOFaceType &face, int innerLevel, int innerTwist)
Definition faceutility_imp.cc:43
bool conformingRefinement() const
return true if conforming refinement is enabled
Definition faceutility.hh:138
const GEOElementType & innerEntity() const
Returns the inner element at that face.
Definition faceutility_imp.cc:327
int innerALUFaceIndex() const
Local number of the face in inner element (ALU3dGrid reference element)
Definition faceutility_imp.cc:427
const BNDFaceType & innerFace() const
Definition faceutility_imp.cc:343
bool ghostBoundary() const
Definition faceutility_imp.cc:312
ALU3dGridFaceInfo(const bool levelIntersection=false)
constructor creating empty face info
Definition faceutility_imp.cc:10
int duneTwist(const int faceIdx, const int aluTwist) const
Twist of the face seen from the inner element.
Definition faceutility_imp.cc:388
bool neighbor() const
returns true if outside is something meaningfull
Definition faceutility_imp.cc:306
bool innerBoundary() const
returns true if inside is a ghost entity
Definition faceutility_imp.cc:291
int boundaryId() const
return boundary id if intersection is with domain boundary
Definition faceutility_imp.cc:373
bool boundary() const
returns true if the face lies on the domain boundary
Definition faceutility_imp.cc:301
const GEOFaceType & face() const
Returns the ALU3dGrid face.
Definition faceutility_imp.cc:320
int innerTwist() const
Twist of the face seen from the inner element.
Definition faceutility_imp.cc:379
bool outerBoundary() const
Definition faceutility_imp.cc:296
int outerALUFaceIndex() const
Local number of the face in outer element (ALU3dGrid reference element)
Definition faceutility_imp.cc:432
const GEOElementType & outerEntity() const
Definition faceutility_imp.cc:335
ImplTraits::GEOPeriodicType GEOPeriodicType
Definition faceutility.hh:52
ImplTraits::GEOFaceType GEOFaceType
Definition faceutility.hh:50
ConformanceState conformanceState() const
Description of conformance on the face.
Definition faceutility_imp.cc:438
const BNDFaceType & boundaryFace() const
Definition faceutility_imp.cc:351
int segmentIndex() const
return boundary segment index if intersection is with domain boundary
Definition faceutility_imp.cc:366
ImplTraits::BNDFaceType BNDFaceType
Definition faceutility.hh:55
Definition faceutility.hh:235
FieldMatrix< alu3d_ctype, numVerticesPerFace, dimworld > CoordinateType
Definition faceutility.hh:261
bool generatedLocal_
Definition faceutility.hh:304
CoordinateType coordsNeighborLocal_
Definition faceutility.hh:301
const CoordinateType & intersectionNeighborLocal() const
Definition faceutility_imp.cc:511
@ numVerticesPerFace
Definition faceutility.hh:254
ReferenceElement< alu3d_ctype, 3 > ReferenceElementType
Definition faceutility.hh:249
static const ReferenceElementType & getReferenceElement()
Definition faceutility.hh:306
const ConnectorType & connector_
Definition faceutility.hh:298
ALU3dGridGeometricFaceInfoBase(const ConnectorType &)
Definition faceutility_imp.cc:474
ElementTopologyMapping< type > ElementTopo
Definition faceutility.hh:239
@ OUTER
Definition faceutility.hh:253
@ INNER
Definition faceutility.hh:253
void resetFaceGeom()
reset status of faceGeomInfo
Definition faceutility_imp.cc:485
CoordinateType coordsSelfLocal_
Definition faceutility.hh:300
const CoordinateType & intersectionSelfLocal() const
Definition faceutility_imp.cc:503
bool generatedGlobal_
Definition faceutility.hh:303
Definition faceutility.hh:326
ALU3dGridGeometricFaceInfoTetra(const ConnectorType &ctor)
Definition faceutility_imp.cc:522
void resetFaceGeom()
reset status of faceGeomInfo
Definition faceutility_imp.cc:528
void buildGlobalGeom(GeometryImp &geo) const
update global geometry
Definition faceutility_imp.cc:544
NormalType & outerNormal(const FieldVector< alu3d_ctype, 2 > &local) const
Definition faceutility_imp.cc:562
ALU3dGridFaceInfo< dim, dimworld, tetra, Comm >::GEOFaceType GEOFaceType
Definition faceutility.hh:333
Definition faceutility.hh:370
ALU3dGridFaceInfo< dim, dimworld, hexa, Comm >::GEOFaceType GEOFaceType
Definition faceutility.hh:377
void buildGlobalGeom(GeometryImp &geo) const
update global geometry
Definition faceutility_imp.cc:617
NormalType & outerNormal(const FieldVector< alu3d_ctype, 2 > &local) const
Definition faceutility_imp.cc:635
void resetFaceGeom()
reset status of faceGeomInfo
Definition faceutility_imp.cc:599
ALU3dGridGeometricFaceInfoHexa(const ConnectorType &)
Definition faceutility_imp.cc:591
ALU3dGridFaceInfo< 2, dimworld, tetra, Comm >::GEOElementType GEOElementType
Definition faceutility.hh:522
ALU3dGridFaceInfo< 2, dimworld, tetra, Comm >::GEOFaceType GEOFaceType
Definition faceutility.hh:521
ALU3dGridFaceInfo< 2, dimworld, hexa, Comm >::GEOFaceType GEOFaceType
Definition faceutility.hh:565
Definition topology.hh:40
Definition topology.hh:151