1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ALBERTA_ELEMENTINFO_HH
4 #define DUNE_ALBERTA_ELEMENTINFO_HH
5 
6 /** \file
7  *  \author Martin Nolte
8  *  \brief  provides a wrapper for ALBERTA's el_info structure
9  */
10 
11 #include <cassert>
12 #include <vector>
13 #include <utility>
14 
15 #include <dune/grid/albertagrid/geometrycache.hh>
16 #include <dune/grid/albertagrid/macroelement.hh>
17 
18 #if HAVE_ALBERTA
19 
20 namespace Dune
21 {
22 
23   namespace Alberta
24   {
25 
26     // External Forward Declarations
27     // -----------------------------
28 
29     template< int dim >
30     class MeshPointer;
31 
32     struct BasicNodeProjection;
33 
34 
35 
36     // ElementInfo
37     // -----------
38 
39     template< int dim >
40     class ElementInfo
41     {
42       struct Instance;
43       class Stack;
44 
45       template< int >
46       struct Library;
47 
48       typedef Instance *InstancePtr;
49 
50     public:
51       static const int dimension = dim;
52 
53       static const int numVertices = NumSubEntities< dimension, dimension >::value;
54       static const int numFaces = NumSubEntities< dimension, 1 >::value;
55 
56       typedef Alberta::MacroElement< dimension > MacroElement;
57       typedef Alberta::MeshPointer< dimension > MeshPointer;
58       typedef Alberta::FillFlags< dimension > FillFlags;
59 
60       static const int maxNeighbors = N_NEIGH_MAX;
61 
62       static const int maxLevelNeighbors = Library< dimWorld >::maxLevelNeighbors;
63 
64 #if !DUNE_ALBERTA_CACHE_COORDINATES
65       typedef GeometryCacheProxy< dim > GeometryCache;
66 #endif
67 
68       struct Seed;
69 
70     private:
71       explicit ElementInfo ( const InstancePtr &instance );
72 
73     public:
74       ElementInfo ();
75       ElementInfo ( const MeshPointer &mesh, const MacroElement &macroElement,
76                     typename FillFlags::Flags fillFlags = FillFlags::standard );
77       ElementInfo ( const MeshPointer &mesh, const Seed &seed,
78                     typename FillFlags::Flags fillFlags = FillFlags::standard );
79       ElementInfo ( const ElementInfo &other );
80       ElementInfo ( ElementInfo&& other );
81 
82       ~ElementInfo ();
83 
84       ElementInfo &operator= ( const ElementInfo &other );
85       ElementInfo &operator= ( ElementInfo &&other );
86 
operator bool() const87       explicit operator bool () const { return (instance_ != null()); }
88 
89       bool operator== ( const ElementInfo &other ) const;
90       bool operator!= ( const ElementInfo &other ) const;
91 
92       const MacroElement &macroElement () const;
93       ElementInfo father () const;
94       int indexInFather () const;
95       ElementInfo child ( int i ) const;
96       bool isLeaf () const;
97 
98       Seed seed () const;
99 
100       MeshPointer mesh () const;
101 
102       bool mightVanish () const;
103 
104       int level () const;
105       // see ALBERTA documentation for definition of element type
106       // values are 0, 1, 2
107       int type () const;
108 
109       int getMark () const;
110       void setMark ( int refCount ) const;
111 
112       bool hasLeafNeighbor ( const int face ) const;
113       ElementInfo leafNeighbor ( const int face ) const;
114 
115       /* obtain all level neighbors of a face
116        *
117        * param[in]  face            face for which the neighbors are desired
118        * param[out] neighbor        array storing the neighbors
119        * param[out] faceInNeighbor  array storing the faces in neighbor
120        *                            (-1, if this neighbor does not exist)
121        *
122        * returns (potential) number of neighbors (i.e., the number of valid
123        *         entries in the output arrays
124        */
125       int levelNeighbors ( const int face, ElementInfo (&neighbor)[ maxLevelNeighbors ], int (&faceInNeighbor)[ maxLevelNeighbors ] ) const;
126 
127       template< int codim >
128       int twist ( int subEntity ) const;
129       int twistInNeighbor ( int face ) const;
130       bool isBoundary ( int face ) const;
131       int boundaryId ( int face ) const;
132       AffineTransformation *transformation ( int face ) const;
133       BasicNodeProjection *boundaryProjection ( int face ) const;
134 
135       bool hasCoordinates () const;
136       const GlobalVector &coordinate ( int vertex ) const;
137 #if !DUNE_ALBERTA_CACHE_COORDINATES
geometryCache() const138       GeometryCache geometryCache () const
139       {
140         return GeometryCache( instance_->geometryCache, instance_->elInfo );
141       }
142 #endif
143 
144       template< class Functor >
145       void hierarchicTraverse ( Functor &functor ) const;
146 
147       template< class Functor >
148       void leafTraverse ( Functor &functor ) const;
149 
150       const Element *element () const;
151       const Element *neighbor ( int face ) const;
152       Element *el () const;
153       ALBERTA EL_INFO &elInfo () const;
154 
155       static ElementInfo
156       createFake ( const MeshPointer &mesh,
157                    const Element *element, int level, int type = 0 );
158       static ElementInfo createFake ( const ALBERTA EL_INFO &elInfo );
159 
160     private:
161       static bool isLeaf ( Element *element );
162       static bool mightVanish ( Element *element, int depth );
163 
164       static void fill ( Mesh *mesh, const ALBERTA MACRO_EL *mel, ALBERTA EL_INFO &elInfo );
165       static void fill ( int ichild, const ALBERTA EL_INFO &parentInfo, ALBERTA EL_INFO &elInfo );
166 
167       void addReference () const;
168       void removeReference () const;
169 
170       static InstancePtr null ();
171       static Stack &stack ();
172 
173       InstancePtr instance_;
174     };
175 
176 
177 
178     // ElementInfo::Instance
179     // ---------------------
180 
181     template< int dim >
182     struct ElementInfo< dim >::Instance
183     {
184       ALBERTA EL_INFO elInfo;
185       unsigned int refCount;
186 
parentDune::Alberta::ElementInfo::Instance187       InstancePtr &parent ()
188       {
189         return parent_;
190       }
191 
192     private:
193       InstancePtr parent_;
194 
195 #if !DUNE_ALBERTA_CACHE_COORDINATES
196     public:
197       Alberta::GeometryCache< dim > geometryCache;
198 #endif
199     };
200 
201 
202 
203     // ElementInfo::Stack
204     // ------------------
205 
206     template< int dim >
207     class ElementInfo< dim >::Stack
208     {
209       InstancePtr top_;
210       Instance null_;
211 
212     public:
213       Stack ();
214       ~Stack ();
215 
216       InstancePtr allocate ();
217       void release ( InstancePtr &p );
218       InstancePtr null ();
219     };
220 
221 
222 
223     // ElementInfo::Library
224     // --------------------
225 
226     template< int dim >
227     template< int >
228     struct ElementInfo< dim >::Library
229     {
230       typedef Alberta::ElementInfo< dim > ElementInfo;
231 
232       static const int maxLevelNeighbors = (1 << (dim-1));
233 
234       static int
235       leafNeighbor ( const ElementInfo &element, const int face, ElementInfo &neighbor );
236 
237       static int
238       levelNeighbors ( const ElementInfo &element, const int face,
239                        ElementInfo (&neighbor)[ maxLevelNeighbors ], int (&faceInNeighbor)[ maxLevelNeighbors ] );
240 
241     private:
242       static int
243       macroNeighbor ( const ElementInfo &element, const int face, ElementInfo &neighbor );
244     };
245 
246 
247 
248     // ElementInfo::Seed
249     // -----------------
250 
251     template< int dim >
252     struct ElementInfo< dim >::Seed
253     {
SeedDune::Alberta::ElementInfo::Seed254       Seed ()
255         : macroIndex_( -1 ), level_( 0 ), path_( 0 )
256       {}
257 
SeedDune::Alberta::ElementInfo::Seed258       Seed ( const int macroIndex, const int level, const unsigned long path )
259         : macroIndex_( macroIndex ), level_( level ), path_( path )
260       {}
261 
operator ==Dune::Alberta::ElementInfo::Seed262       bool operator== ( const Seed &other ) const
263       {
264         return (macroIndex() == other.macroIndex()) && (level() == other.level()) && (path() == other.path());
265       }
266 
operator <Dune::Alberta::ElementInfo::Seed267       bool operator< ( const Seed &other ) const
268       {
269         const bool ml = (macroIndex() < other.macroIndex());
270         const bool me = (macroIndex() == other.macroIndex());
271         const bool ll = (level() < other.level());
272         const bool le = (level() == other.level());
273         const bool pl = (path() < other.path());
274         return ml | (me & (ll | (le & pl)));
275       }
276 
operator !=Dune::Alberta::ElementInfo::Seed277       bool operator!= ( const Seed &other ) const { return !(*this == other); }
operator <=Dune::Alberta::ElementInfo::Seed278       bool operator<= ( const Seed &other ) const { return !(other < *this); }
operator >Dune::Alberta::ElementInfo::Seed279       bool operator> ( const Seed &other ) const { return (other < *this); }
operator >=Dune::Alberta::ElementInfo::Seed280       bool operator>= ( const Seed &other ) const { return !(*this < other); }
281 
isValidDune::Alberta::ElementInfo::Seed282       bool isValid ( ) const { return macroIndex_ != -1; }
283 
macroIndexDune::Alberta::ElementInfo::Seed284       int macroIndex () const { return macroIndex_; }
levelDune::Alberta::ElementInfo::Seed285       int level () const { return level_; }
pathDune::Alberta::ElementInfo::Seed286       unsigned long path () const { return path_; }
287 
288     private:
289       int macroIndex_;
290       int level_;
291       unsigned long path_;
292     };
293 
294 
295 
296     // Implementation of ElementInfo
297     // -----------------------------
298 
299     template< int dim >
ElementInfo(const InstancePtr & instance)300     inline ElementInfo< dim >::ElementInfo ( const InstancePtr &instance )
301       : instance_( instance )
302     {
303       addReference();
304     }
305 
306 
307     template< int dim >
ElementInfo()308     inline ElementInfo< dim >::ElementInfo ()
309       : instance_( null() )
310     {
311       addReference();
312     }
313 
314 
315     template< int dim >
316     inline ElementInfo< dim >
ElementInfo(const MeshPointer & mesh,const MacroElement & macroElement,typename FillFlags::Flags fillFlags)317     ::ElementInfo ( const MeshPointer &mesh, const MacroElement &macroElement,
318                     typename FillFlags::Flags fillFlags )
319     {
320       instance_ = stack().allocate();
321       instance_->parent() = null();
322       ++(instance_->parent()->refCount);
323 
324       addReference();
325 
326       elInfo().fill_flag = fillFlags;
327 
328       // Alberta fills opp_vertex only if there is a neighbor
329       for( int k = 0; k < maxNeighbors; ++k )
330         elInfo().opp_vertex[ k ] = -1;
331 
332       fill( mesh, &macroElement, elInfo() );
333     }
334 
335 
336     template< int dim >
337     inline ElementInfo< dim >
ElementInfo(const MeshPointer & mesh,const Seed & seed,typename FillFlags::Flags fillFlags)338     ::ElementInfo ( const MeshPointer &mesh, const Seed &seed,
339                     typename FillFlags::Flags fillFlags )
340     {
341       instance_ = stack().allocate();
342       instance_->parent() = null();
343       ++(instance_->parent()->refCount);
344 
345       addReference();
346 
347       // fill in macro element info
348       elInfo().fill_flag = fillFlags;
349 
350       // Alberta fills opp_vertex only if there is a neighbor
351       for( int k = 0; k < maxNeighbors; ++k )
352         elInfo().opp_vertex[ k ] = -1;
353 
354       fill( mesh, ((Mesh *)mesh)->macro_els + seed.macroIndex(), elInfo() );
355 
356       // traverse the seed's path
357       unsigned long path = seed.path();
358       for( int i = 0; i < seed.level(); ++i )
359       {
360         InstancePtr child = stack().allocate();
361         child->parent() = instance_;
362 
363         // Alberta fills opp_vertex only if there is a neighbor
364         for( int k = 0; k < maxNeighbors; ++k )
365           child->elInfo.opp_vertex[ k ] = -2;
366 
367         fill( path & 1, elInfo(), child->elInfo );
368 
369         instance_ = child;
370         addReference();
371 
372         path = path >> 1;
373       }
374 
375       assert( this->seed() == seed );
376     }
377 
378 
379     template< int dim >
ElementInfo(const ElementInfo & other)380     inline ElementInfo< dim >::ElementInfo ( const ElementInfo &other )
381       : instance_( other.instance_ )
382     {
383       addReference();
384     }
385 
386     template< int dim >
ElementInfo(ElementInfo && other)387     inline ElementInfo< dim >::ElementInfo ( ElementInfo &&other )
388       : instance_( NULL )
389     {
390       using std::swap;
391       swap( instance_, other.instance_ );
392     }
393 
394     template< int dim >
~ElementInfo()395     inline ElementInfo< dim >::~ElementInfo ()
396     {
397       removeReference();
398     }
399 
400 
401     template< int dim >
402     inline ElementInfo< dim > &
operator =(const ElementInfo<dim> & other)403     ElementInfo< dim >::operator= ( const ElementInfo< dim > &other )
404     {
405       other.addReference();
406       removeReference();
407       instance_ = other.instance_;
408       return *this;
409     }
410 
411     template< int dim >
412     inline ElementInfo< dim > &
operator =(ElementInfo<dim> && other)413     ElementInfo< dim >::operator= ( ElementInfo< dim > &&other )
414     {
415       using std::swap;
416       swap( instance_, other.instance_ );
417       return *this;
418     }
419 
420     template< int dim >
421     inline bool
operator ==(const ElementInfo<dim> & other) const422     ElementInfo< dim >::operator== ( const ElementInfo< dim > &other ) const
423     {
424       return (instance_->elInfo.el == other.instance_->elInfo.el);
425     }
426 
427 
428     template< int dim >
429     inline bool
operator !=(const ElementInfo<dim> & other) const430     ElementInfo< dim >::operator!= ( const ElementInfo< dim > &other ) const
431     {
432       return (instance_->elInfo.el != other.instance_->elInfo.el);
433     }
434 
435 
436     template< int dim >
437     inline const typename ElementInfo< dim >::MacroElement &
macroElement() const438     ElementInfo< dim >::macroElement () const
439     {
440       assert( !!(*this) );
441       assert( elInfo().macro_el != NULL );
442       return static_cast< const MacroElement & >( *(elInfo().macro_el) );
443     }
444 
445 
446     template< int dim >
father() const447     inline ElementInfo< dim > ElementInfo< dim >::father () const
448     {
449       assert( !!(*this) );
450       return ElementInfo< dim >( instance_->parent() );
451     }
452 
453 
454     template< int dim >
indexInFather() const455     inline int ElementInfo< dim >::indexInFather () const
456     {
457       const Element *element = elInfo().el;
458       const Element *father = elInfo().parent->el;
459       assert( father != NULL );
460 
461       const int index = (father->child[ 0 ] == element ? 0 : 1);
462       assert( father->child[ index ] == element );
463       return index;
464     }
465 
466 
467     template< int dim >
child(int i) const468     inline ElementInfo< dim > ElementInfo< dim >::child ( int i ) const
469     {
470       assert( !isLeaf() );
471 
472       InstancePtr child = stack().allocate();
473       child->parent() = instance_;
474       addReference();
475 
476       // Alberta fills opp_vertex only if there is a neighbor
477       for( int k = 0; k < maxNeighbors; ++k )
478         child->elInfo.opp_vertex[ k ] = -2;
479 
480       fill( i, elInfo(), child->elInfo );
481       return ElementInfo< dim >( child );
482     }
483 
484 
485     template< int dim >
isLeaf() const486     inline bool ElementInfo< dim >::isLeaf () const
487     {
488       assert( !(*this) == false );
489       return isLeaf( el() );
490     }
491 
492 
493     template< int dim >
seed() const494     inline typename ElementInfo< dim >::Seed ElementInfo< dim >::seed () const
495     {
496       assert( !!(*this) );
497 
498       int level = 0;
499       unsigned long path = 0;
500       for( InstancePtr p = instance_; p->parent() != null(); p = p->parent() )
501       {
502         const Element *element = p->elInfo.el;
503         const Element *father = p->parent()->elInfo.el;
504         const unsigned long child = static_cast< unsigned long >( father->child[ 1 ] == element );
505         path = (path << 1) | child;
506         ++level;
507       }
508 
509       if( level != elInfo().level )
510         DUNE_THROW( NotImplemented, "Seed for fake elements not implemented." );
511 
512       return Seed( macroElement().index, level, path );
513     }
514 
515 
516     template< int dim >
mesh() const517     inline typename ElementInfo< dim >::MeshPointer ElementInfo< dim >::mesh () const
518     {
519       return MeshPointer( elInfo().mesh );
520     }
521 
522 
523     template< int dim >
mightVanish() const524     inline bool ElementInfo< dim >::mightVanish () const
525     {
526       return mightVanish( el(), 0 );
527     }
528 
529 
530     template< int dim >
level() const531     inline int ElementInfo< dim >::level () const
532     {
533       return elInfo().level;
534     }
535 
536 
537     template< int dim >
type() const538     inline int ElementInfo< dim >::type () const
539     {
540       return 0;
541     }
542 
543 
544     template<>
type() const545     inline int ElementInfo< 3 >::type () const
546     {
547       return instance_->elInfo.el_type;
548     }
549 
550 
551     template< int dim >
getMark() const552     inline int ElementInfo< dim >::getMark () const
553     {
554       return el()->mark;
555     }
556 
557 
558     template< int dim >
setMark(int refCount) const559     inline void ElementInfo< dim >::setMark ( int refCount ) const
560     {
561       assert( isLeaf() );
562       assert( (refCount >= -128) && (refCount < 127) );
563       el()->mark = refCount;
564     }
565 
566 
567     template< int dim >
hasLeafNeighbor(const int face) const568     inline bool ElementInfo< dim >::hasLeafNeighbor ( const int face ) const
569     {
570       assert( !!(*this) );
571       assert( (face >= 0) && (face < maxNeighbors) );
572 
573       assert( (elInfo().fill_flag & FillFlags::boundaryId) != 0 );
574       const int macroFace = elInfo().macro_wall[ face ];
575       if( macroFace >= 0 )
576         return (macroElement().neighbor( macroFace ) != NULL);
577       else
578         return true;
579     }
580 
581 
582     template< int dim >
leafNeighbor(const int face) const583     inline ElementInfo< dim > ElementInfo< dim >::leafNeighbor ( const int face ) const
584     {
585       assert( (face >= 0) && (face < numFaces) );
586       ElementInfo neighbor;
587       Library< dimWorld >::leafNeighbor( *this, face, neighbor );
588       return neighbor;
589     }
590 
591 
592     template< int dim >
593     inline int ElementInfo< dim >
levelNeighbors(const int face,ElementInfo (& neighbor)[maxLevelNeighbors],int (& faceInNeighbor)[maxLevelNeighbors]) const594     ::levelNeighbors ( const int face, ElementInfo (&neighbor)[ maxLevelNeighbors ], int (&faceInNeighbor)[ maxLevelNeighbors ] ) const
595     {
596       assert( (face >= 0) && (face < numFaces) );
597       return Library< dimWorld >::levelNeighbors( *this, face, neighbor, faceInNeighbor );
598     }
599 
600 
601     template< int dim >
602     template< int codim >
twist(int subEntity) const603     inline int ElementInfo< dim >::twist ( int subEntity ) const
604     {
605       return Twist< dim, dim-codim >::twist( element(), subEntity );
606     }
607 
608 
609     template< int dim >
twistInNeighbor(const int face) const610     inline int ElementInfo< dim >::twistInNeighbor ( const int face ) const
611     {
612       assert( neighbor( face ) != NULL );
613       return Twist< dim, dim-1 >::twist( neighbor( face ), elInfo().opp_vertex[ face ] );
614     }
615 
616 
617     template< int dim >
isBoundary(int face) const618     inline bool ElementInfo< dim >::isBoundary ( int face ) const
619     {
620       assert( !!(*this) );
621       assert( (face >= 0) && (face < maxNeighbors) );
622 
623       assert( (elInfo().fill_flag & FillFlags::boundaryId) != 0 );
624       const int macroFace = elInfo().macro_wall[ face ];
625       if( macroFace >= 0 )
626         return macroElement().isBoundary( macroFace );
627       else
628         return false;
629     }
630 
631 
632     template< int dim >
boundaryId(int face) const633     inline int ElementInfo< dim >::boundaryId ( int face ) const
634     {
635       assert( !!(*this) );
636       assert( (face >= 0) && (face < N_WALLS_MAX) );
637 
638       assert( (elInfo().fill_flag & FillFlags::boundaryId) != 0 );
639       const int macroFace = elInfo().macro_wall[ face ];
640       const int id = macroElement().boundaryId( macroFace );
641       // this assertion is only allowed, if FILL_BOUND is set
642       // assert( id == elInfo().wall_bound[ face ] );
643       return id;
644     }
645 
646 
647     template< int dim >
648     inline AffineTransformation *
transformation(int face) const649     ElementInfo< dim >::transformation ( int face ) const
650     {
651       assert( !!(*this) );
652       assert( (face >= 0) && (face < N_WALLS_MAX) );
653 
654       assert( (elInfo().fill_flag & FillFlags::boundaryId) != 0 );
655       const int macroFace = elInfo().macro_wall[ face ];
656       return (macroFace < 0 ? NULL : macroElement().wall_trafo[ macroFace ]);
657     }
658 
659 
660     template< int dim >
661     inline BasicNodeProjection *
boundaryProjection(int face) const662     ElementInfo< dim >::boundaryProjection ( int face ) const
663     {
664       assert( !!(*this) );
665       assert( (face >= 0) && (face < N_WALLS_MAX) );
666 
667       assert( (elInfo().fill_flag & FillFlags::boundaryId) != 0 );
668       const int macroFace = elInfo().macro_wall[ face ];
669       if( macroFace >= 0 )
670         return static_cast< BasicNodeProjection * >( macroElement().projection[ macroFace+1 ] );
671       else
672         return 0;
673     }
674 
675 
676     template< int dim >
hasCoordinates() const677     inline bool ElementInfo< dim >::hasCoordinates () const
678     {
679       return ((elInfo().fill_flag & FillFlags::coords) != 0);
680     }
681 
682     template< int dim >
coordinate(int vertex) const683     inline const GlobalVector &ElementInfo< dim >::coordinate ( int vertex ) const
684     {
685       assert( hasCoordinates() );
686       assert( (vertex >= 0) && (vertex < numVertices) );
687       return elInfo().coord[ vertex ];
688     }
689 
690 
691     template< int dim >
692     template< class Functor >
hierarchicTraverse(Functor & functor) const693     inline void ElementInfo< dim >::hierarchicTraverse ( Functor &functor ) const
694     {
695       functor( *this );
696       if( !isLeaf() )
697       {
698         child( 0 ).hierarchicTraverse( functor );
699         child( 1 ).hierarchicTraverse( functor );
700       }
701     }
702 
703 
704     template< int dim >
705     template< class Functor >
leafTraverse(Functor & functor) const706     inline void ElementInfo< dim >::leafTraverse ( Functor &functor ) const
707     {
708       if( !isLeaf() )
709       {
710         child( 0 ).leafTraverse( functor );
711         child( 1 ).leafTraverse( functor );
712       }
713       else
714         functor( *this );
715     }
716 
717 
718     template< int dim >
element() const719     inline const Element *ElementInfo< dim >::element () const
720     {
721       return elInfo().el;
722     }
723 
724 
725     template< int dim >
neighbor(int face) const726     inline const Element *ElementInfo< dim >::neighbor ( int face ) const
727     {
728       assert( (face >= 0) && (face < numFaces) );
729       assert( (elInfo().fill_flag & FillFlags::neighbor) != 0 );
730       return elInfo().neigh[ face ];
731     }
732 
733 
734     template< int dim >
el() const735     inline Element *ElementInfo< dim >::el () const
736     {
737       return elInfo().el;
738     }
739 
740 
741     template< int dim >
elInfo() const742     inline ALBERTA EL_INFO &ElementInfo< dim >::elInfo () const
743     {
744       return (instance_->elInfo);
745     }
746 
747 
748     template< int dim >
749     inline ElementInfo< dim >
createFake(const MeshPointer & mesh,const Element * element,int level,int type)750     ElementInfo< dim >::createFake ( const MeshPointer &mesh,
751                                      const Element *element, int level, int type )
752     {
753       InstancePtr instance = stack().allocate();
754       instance->parent() = null();
755       ++(instance->parent()->refCount);
756 
757       instance->elInfo.mesh = mesh;
758       instance->elInfo.macro_el = NULL;
759       instance->elInfo.el = const_cast< Element * >( element );
760       instance->elInfo.parent = NULL;
761       instance->elInfo.fill_flag = FillFlags::nothing;
762       instance->elInfo.level = level;
763       instance->elInfo.el_type = type;
764 
765       return ElementInfo< dim >( instance );
766     }
767 
768 
769     template< int dim >
770     inline ElementInfo< dim >
createFake(const ALBERTA EL_INFO & elInfo)771     ElementInfo< dim >::createFake ( const ALBERTA EL_INFO &elInfo )
772     {
773       InstancePtr instance = stack().allocate();
774       instance->parent() = null();
775       ++(instance->parent()->refCount);
776 
777       instance->elInfo = elInfo;
778       return ElementInfo< dim >( instance );
779     }
780 
781 
782     template< int dim >
isLeaf(Element * element)783     inline bool ElementInfo< dim >::isLeaf ( Element *element )
784     {
785       return IS_LEAF_EL( element );
786     }
787 
788 
789     template< int dim >
mightVanish(Alberta::Element * element,int depth)790     inline bool ElementInfo< dim >::mightVanish ( Alberta::Element *element, int depth )
791     {
792       if( isLeaf( element ) )
793         return (element->mark < depth);
794       else
795         return (mightVanish( element->child[ 0 ], depth-1 ) && mightVanish( element->child[ 1 ], depth-1 ));
796     }
797 
798 
799     template< int dim >
800     inline void ElementInfo< dim >
fill(Mesh * mesh,const ALBERTA MACRO_EL * mel,ALBERTA EL_INFO & elInfo)801     ::fill ( Mesh *mesh, const ALBERTA MACRO_EL *mel, ALBERTA EL_INFO &elInfo )
802     {
803       ALBERTA fill_macro_info( mesh, mel, &elInfo );
804     }
805 
806     template< int dim >
807     inline void ElementInfo< dim >
fill(int ichild,const ALBERTA EL_INFO & parentInfo,ALBERTA EL_INFO & elInfo)808     ::fill ( int ichild, const ALBERTA EL_INFO &parentInfo, ALBERTA EL_INFO &elInfo )
809     {
810       ALBERTA fill_elinfo( ichild, FILL_ANY, &parentInfo, &elInfo );
811     }
812 
813 
814     template< int dim >
addReference() const815     inline void ElementInfo< dim >::addReference () const
816     {
817       ++(instance_->refCount);
818     }
819 
820 
821     template< int dim >
removeReference() const822     inline void ElementInfo< dim >::removeReference () const
823     {
824       // short-circuit for rvalues that have been drained as argument to a move operation
825       if ( !instance_ )
826         return;
827       // this loop breaks when instance becomes null()
828       for( InstancePtr instance = instance_; --(instance->refCount) == 0; )
829       {
830         const InstancePtr parent = instance->parent();
831         stack().release( instance );
832         instance = parent;
833       }
834     }
835 
836 
837     template< int dim >
838     inline typename ElementInfo< dim >::InstancePtr
null()839     ElementInfo< dim >::null ()
840     {
841       return stack().null();
842     }
843 
844 
845     template< int dim >
846     inline typename ElementInfo< dim >::Stack &
stack()847     ElementInfo< dim >::stack ()
848     {
849       static Stack s;
850       return s;
851     }
852 
853 
854 
855     // Implementation of ElementInfo::Stack
856     // ------------------------------------
857 
858     template< int dim >
Stack()859     inline ElementInfo< dim >::Stack::Stack ()
860       : top_( 0 )
861     {
862       null_.elInfo.el = NULL;
863       null_.refCount = 1;
864       null_.parent() = 0;
865     }
866 
867 
868     template< int dim >
~Stack()869     inline ElementInfo< dim >::Stack::~Stack ()
870     {
871       while( top_ != 0 )
872       {
873         InstancePtr p = top_;
874         top_ = p->parent();
875         delete p;
876       }
877     }
878 
879 
880     template< int dim >
881     inline typename ElementInfo< dim >::InstancePtr
allocate()882     ElementInfo< dim >::Stack::allocate ()
883     {
884       InstancePtr p = top_;
885       if( p != 0 )
886         top_ = p->parent();
887       else
888         p = new Instance;
889       p->refCount = 0;
890       return p;
891     }
892 
893 
894     template< int dim >
release(InstancePtr & p)895     inline void ElementInfo< dim >::Stack::release ( InstancePtr &p )
896     {
897       assert( (p != null()) && (p->refCount == 0) );
898       p->parent() = top_;
899       top_ = p;
900     }
901 
902 
903     template< int dim >
904     inline typename ElementInfo< dim >::InstancePtr
null()905     ElementInfo< dim >::Stack::null ()
906     {
907       return &null_;
908     }
909 
910   } // namespace Alberta
911 
912 } // namespace Dune
913 
914 #endif // #if HAVE_ALBERTA
915 
916 #endif // #ifndef DUNE_ALBERTA_ELEMENTINFO_HH
917