1 #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH
2 #define DUNE_ALUGRIDGEOMETRYSTORAGE_HH
3 
4 #include <array>
5 #include <memory>
6 
7 #include <dune/common/exceptions.hh>
8 
9 #include <dune/grid/common/grid.hh>
10 #include <dune/grid/common/gridfactory.hh>
11 
12 #include <dune/alugrid/common/declaration.hh>
13 #include <dune/alugrid/3d/alu3dinclude.hh>
14 
15 namespace Dune
16 {
17   template< class Grid >
18   class ReferenceGridFactory;
19 
20 
21   template< class GridImp, class GeometryImpl, int nChild >
22   class ALULocalGeometryStorage
23   {
24     typedef ALULocalGeometryStorage< GridImp, GeometryImpl, nChild > ThisType;
25 
26     // array with pointers to the geometries
27     std::array< GeometryImpl, nChild > geoms_;
28 
29     // count local geometry creation
30     int count_;
31 
32     // true if geoms have been initialized
33     bool initialized_;
34 
35     // type of grid impl
36     typedef typename GridImp :: ctype ctype;
37     enum{ dimension       = GridImp :: dimension };
38     enum{ dimensionworld  = GridImp :: dimensionworld };
39 
40     template <int dummy, int dim, int dimworld, int >
41     struct CreateGeometries;
42 
43     template <int dummy, int dimworld>
44     struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE triangle >
45     {
46       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries47       static inline void createGeometries(Storage& storage,
48                                    const GeometryType& type,
49                                    const bool nonConform )
50       {
51         if( nonConform )
52         {
53           typedef ALUGrid< 2, dimworld, simplex, nonconforming, ALUGridNoComm > Grid;
54           storage.template createGeometries< Grid > (type);
55         }
56         else
57         {
58           typedef ALUGrid< 2, dimworld, simplex, conforming, ALUGridNoComm > Grid;
59           storage.template createGeometries< Grid > (type);
60         }
61       }
62     };
63 
64     template <int dummy, int dimworld>
65     struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE tetra >
66     {
67       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries68       static inline void createGeometries(Storage& storage,
69                                           const GeometryType& type,
70                                           const bool nonConform )
71       {
72         if( nonConform )
73         {
74           typedef ALUGrid< 2, dimworld, simplex, nonconforming, ALUGridNoComm > Grid;
75           storage.template createGeometries< Grid > (type);
76         }
77         else
78         {
79           typedef ALUGrid< 2, dimworld, simplex, conforming, ALUGridNoComm > Grid;
80           storage.template createGeometries< Grid > (type);
81         }
82       }
83     };
84 
85     template <int dummy>
86     struct CreateGeometries<dummy, 3, 3, ALU3DSPACE tetra >
87     {
88       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries89       static inline void createGeometries(Storage& storage,
90                                           const GeometryType& type,
91                                           const bool nonConform )
92       {
93         alugrid_assert ( nonConform ) ;
94         {
95           typedef ALUGrid< 3, 3, simplex, nonconforming, ALUGridNoComm > Grid;
96           storage.template createGeometries< Grid > (type);
97         }
98         /*
99          // TODO, implement this for refinement of all edges (conforming)
100         else
101         {
102           typedef ALUGrid< 3, 3, simplex, conforming, ALUGridNoComm > Grid;
103           storage.template createGeometries< Grid > (type);
104         }
105         */
106       }
107     };
108 
109     template <int dummy, int dimworld>
110     struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE quadrilateral >
111     {
112       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries113       static void createGeometries(Storage& storage,
114                                    const GeometryType& type,
115                                    const bool nonConform )
116       {
117         alugrid_assert ( nonConform ) ;
118         {
119           typedef ALUGrid< 2, dimworld, cube, nonconforming, ALUGridNoComm > Grid;
120           storage.template createGeometries< Grid > (type);
121         }
122       }
123     };
124 
125     template <int dummy, int dimworld>
126     struct CreateGeometries<dummy, 2, dimworld, ALU3DSPACE hexa >
127     {
128       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries129       static void createGeometries(Storage& storage,
130                                    const GeometryType& type,
131                                    const bool nonConform )
132       {
133         alugrid_assert ( nonConform ) ;
134         {
135           typedef ALUGrid< 2, dimworld, cube, nonconforming, ALUGridNoComm > Grid;
136           storage.template createGeometries< Grid > (type);
137         }
138       }
139     };
140 
141     template <int dummy>
142     struct CreateGeometries<dummy, 3, 3, ALU3DSPACE hexa >
143     {
144       template <class Storage>
createGeometriesDune::ALULocalGeometryStorage::CreateGeometries145       static void createGeometries(Storage& storage,
146                                    const GeometryType& type,
147                                    const bool nonConform )
148       {
149         alugrid_assert ( nonConform );
150         {
151           typedef ALUGrid< 3, 3, cube, nonconforming, ALUGridNoComm > Grid;
152           storage.template createGeometries< Grid > (type);
153         }
154       }
155     };
156 
157   public:
158     // create empty storage
ALULocalGeometryStorage(const GeometryType type,const bool nonConform)159     inline ALULocalGeometryStorage ( const GeometryType type, const bool nonConform )
160     : count_( 0 ), initialized_( false )
161     {
162       // initialize geometries
163       initialize( type, nonConform );
164     }
165 
166     // create empty storage
ALULocalGeometryStorage()167     inline ALULocalGeometryStorage ()
168     : count_( 0 ), initialized_( false )
169     {
170     }
171 
172     // return reference to local geometry
operator [](int child) const173     inline const GeometryImpl& operator [] (int child) const
174     {
175       alugrid_assert ( geomCreated(child) );
176       // this method is not thread safe yet
177       return geoms_[child];
178     }
179 
180     //! access local geometry storage
storage(const GeometryType type,const bool nonConforming)181     static inline const ThisType& storage( const GeometryType type, const bool nonConforming )
182     {
183       if( type.isSimplex() )
184       {
185         // create static variable for this thread
186         static ThisType simplexGeoms( type, nonConforming );
187         return simplexGeoms ;
188       }
189       else
190       {
191         // should be a cube geometry a this point
192         alugrid_assert( type.isCube() );
193 
194         // create static variable
195         static ThisType cubeGeoms( type, nonConforming );
196         return cubeGeoms ;
197       }
198     }
199 
200   protected:
201     // check if geometry has been created
geomCreated(int child) const202     inline bool geomCreated(int child) const { return geoms_[child].valid(); }
203 
204     //! initialize local geometries
initialize(const GeometryType type,const bool nonConform)205     inline bool initialize( const GeometryType type, const bool nonConform )
206     {
207       if( ! initialized_ )
208       {
209         // first set flag, because this method might be called again during
210         // creation of local geometries and then result in an infinite loop
211         initialized_ = true ;
212 
213         // the idea is to create a grid containing the reference element,
214         // refine once and the store the father - child relations
215         CreateGeometries<0, dimension, dimensionworld, GridImp :: elementType >
216           ::createGeometries(*this, type, nonConform);
217         return true;
218       }
219       return false;
220     }
221 
222     template < class Grid >
createGeometries(const GeometryType & type)223     inline void createGeometries(const GeometryType& type)
224     {
225       static bool firstCall = true ;
226       if( firstCall )
227       {
228         firstCall = false ;
229 
230         // create factory for the reference element grid
231         ReferenceGridFactory< Grid > factory;
232 
233         const auto& refElem =
234           Dune::ReferenceElements< ctype, dimension >::general( type );
235 
236         // insert vertices
237         FieldVector<ctype, dimensionworld> pos( 0 );
238         const int vxSize = refElem.size(dimension);
239         for(int i=0; i<vxSize; ++i)
240         {
241           FieldVector<ctype, dimension> position = refElem.position(i, dimension );
242           // copy position
243           for(int d = 0; d<dimension; ++d )
244             pos[ d ] = position[ d ];
245 
246           factory.insertVertex( pos );
247         }
248 
249         std::vector< unsigned int > vertices( vxSize );
250         // create grid with reference element
251         for(size_t i=0; i<vertices.size(); ++i) vertices[ i ] = i;
252         factory.insertElement(type, vertices);
253 
254         std::unique_ptr< Grid > gridPtr( factory.createGrid() );
255         Grid& grid = *gridPtr;
256 
257         // refine once to get children in the reference element
258         const int level = 1;
259         grid.globalRefine( level );
260 
261         {
262           typedef typename Grid :: MacroGridView MacroGridView;
263           MacroGridView macroView = grid.template macroGridView< All_Partition > ();
264           typedef typename MacroGridView :: template Codim< 0 > :: Iterator Iterator;
265 
266           Iterator it = macroView.template begin<0> ();
267 
268           if( it == macroView.template end<0>() )
269             DUNE_THROW(InvalidStateException,"Empty Grid, should contain at least 1 element");
270 
271           typedef typename Iterator :: Entity EntityType;
272 
273           const EntityType& entity = *it;
274           const typename EntityType :: Geometry& geo = entity.geometry();
275           typedef typename EntityType :: HierarchicIterator HierarchicIteratorType;
276           const HierarchicIteratorType end = entity.hend( level );
277 
278           int childNum = 0;
279           for( HierarchicIteratorType child = entity.hbegin( level );
280                child != end; ++child, ++childNum )
281           {
282             create( geo, child->geometry(), childNum );
283           }
284         }
285       }
286     }
287 
288     // create local geometry
289     template< class Geometry >
create(const Geometry & father,const Geometry & son,const int child)290     inline void create ( const Geometry &father,
291                   const Geometry &son,
292                   const int child )
293     {
294       alugrid_assert ( !geomCreated( child ) );
295       alugrid_assert ( (child >= 0) && (child < nChild) );
296 
297       alugrid_assert ( (count_ < nChild) );
298       ++count_;
299 
300       geoms_[ child ].buildGeomInFather( father, son );
301     }
302 
303   };
304 
305 } // namespace Dune
306 
307 #endif // #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH
308