1 #ifndef DUNE_ALUGRID_FROMTOGRIDFACTORY_HH
2 #define DUNE_ALUGRID_FROMTOGRIDFACTORY_HH
3 
4 #include <map>
5 #include <vector>
6 
7 #include <dune/common/version.hh>
8 
9 #include <dune/grid/common/gridfactory.hh>
10 #include <dune/grid/common/exceptions.hh>
11 
12 #include <dune/alugrid/common/alugrid_assert.hh>
13 #include <dune/alugrid/common/declaration.hh>
14 
15 namespace Dune
16 {
17 
18   // External Forward Declarations
19   // -----------------------------
20 
21   template< class ToGrid >
22   class FromToGridFactory;
23 
24 
25   // FromToGridFactory for ALUGrid
26   // -----------------------------
27 
28   template< int dim, int dimworld, ALUGridElementType eltype, ALUGridRefinementType refineType, class Comm >
29   class FromToGridFactory< ALUGrid< dim, dimworld, eltype, refineType, Comm > >
30   {
31   public:
32     // type of grid that is converted to
33     typedef ALUGrid< dim, dimworld, eltype, refineType, Comm > Grid;
34   protected:
35     typedef FromToGridFactory< Grid > This;
36 
37     // typedef grid pointer type based on what the grid factory interface defines
38     typedef decltype(std::declval< Dune::GridFactoryInterface< Grid >* >()->createGrid())  GridPtrType;
39 
40     std::vector< unsigned int > ordering_ ;
41 
42   public:
43     template <class FromGrid, class Vector>
convert(const FromGrid & grid,Vector & cellData,std::vector<unsigned int> & ordering)44     GridPtrType convert( const FromGrid& grid, Vector& cellData, std::vector<unsigned int>& ordering )
45     {
46       int rank = 0;
47 #if HAVE_MPI
48       MPI_Comm_rank( MPI_COMM_WORLD, &rank );
49 #endif
50 
51       // create ALUGrid GridFactory
52       GridFactory< Grid > factory;
53 
54       // only insert elements on rank 0
55       if( rank == 0 )
56       {
57         typedef typename FromGrid  :: LeafGridView GridView ;
58         typedef typename GridView  :: IndexSet  IndexSet ;
59         typedef typename IndexSet  :: IndexType IndexType ;
60         typedef typename GridView  :: template Codim< 0 > :: Iterator ElementIterator ;
61         typedef typename ElementIterator::Entity  Entity ;
62         typedef typename GridView :: IntersectionIterator     IntersectionIterator ;
63         typedef typename IntersectionIterator :: Intersection Intersection ;
64 
65         GridView gridView = grid.leafGridView();
66         const IndexSet &indexSet = gridView.indexSet();
67 
68         // map global vertex ids to local ones
69         std::map< IndexType, unsigned int > vtxMap;
70 
71         const int numVertices = (1 << dim);
72         std::vector< unsigned int > vertices( numVertices );
73         typedef std::pair< Dune::GeometryType, std::vector< unsigned int > > ElementPair;
74         std::vector< ElementPair > elements;
75         if( ! ordering.empty() )
76           elements.resize( ordering.size() );
77 
78         int nextElementIndex = 0;
79         const ElementIterator end = gridView.template end< 0 >();
80         for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
81         {
82           const Entity &entity = *it;
83 
84           // insert vertices and element
85           const typename Entity::Geometry geo = entity.geometry();
86           alugrid_assert( numVertices == geo.corners() );
87           for( int i = 0; i < numVertices; ++i )
88           {
89             const IndexType vtxId = indexSet.subIndex( entity, i, dim );
90             std::pair< typename std::map< IndexType, unsigned int >::iterator, bool > result
91               = vtxMap.insert( std::make_pair( vtxId, vtxMap.size() ) );
92             if( result.second )
93               factory.insertVertex( geo.corner( i ), vtxId );
94             vertices[ i ] = result.first->second;
95           }
96           if( ordering.empty() )
97           {
98             factory.insertElement( entity.type(), vertices );
99           }
100           else
101           {
102             // store element applying the reordering
103             elements[ ordering[ nextElementIndex++ ] ] = ElementPair( entity.type(), vertices ) ;
104           }
105         }
106 
107         if( ! ordering.empty() )
108         {
109           // insert elements using reordered list
110           for( auto it = elements.begin(), end = elements.end(); it != end; ++it )
111           {
112             factory.insertElement( (*it).first, (*it).second );
113           }
114         }
115 
116         nextElementIndex = 0;
117         for( ElementIterator it = gridView.template begin< 0 >(); it != end; ++it )
118         {
119           const Entity &entity = *it;
120 
121           const int elementIndex = ordering.empty() ? nextElementIndex++ : ordering[ nextElementIndex++ ];
122           const IntersectionIterator iend = gridView.iend( entity );
123           for( IntersectionIterator iit = gridView.ibegin( entity ); iit != iend; ++iit )
124           {
125             const Intersection &intersection = *iit;
126             const int faceNumber = intersection.indexInInside();
127             // insert boundary face in case of domain boundary
128             if( intersection.boundary() )
129               factory.insertBoundary( elementIndex, faceNumber );
130 
131             // for parallel grids we can check if we are at a process border
132             if( intersection.neighbor() &&
133                 intersection.outside().partitionType() != InteriorEntity )
134             {
135               // insert process boundary if the neighboring element has a different rank
136               factory.insertProcessBorder( elementIndex, faceNumber );
137             }
138           }
139         }
140       }
141 
142       // create grid pointer (behaving like a shared_ptr)
143       GridPtrType newgrid = factory.createGrid( true, true, std::string("FromToGrid") );
144 
145       if( ! cellData.empty() )
146       {
147         Vector oldCellData( cellData );
148         auto macroView = newgrid->levelGridView( 0 );
149         size_t idx = 0;
150         for( auto it = macroView.template begin<0>(), end = macroView.template end<0>();
151              it != end; ++it, ++idx )
152         {
153           const int insertionIndex = ordering.empty() ?
154             factory.insertionIndex( *it ) : ordering[ factory.insertionIndex( *it ) ];                                                        ;
155           cellData[ idx ] = oldCellData[ insertionIndex ] ;
156         }
157       }
158 
159       // store the ordering from the factory, if it was not provided
160       if( ordering.empty() )
161         ordering = factory.ordering();
162 
163 #if HAVE_MPI
164       MPI_Barrier( MPI_COMM_WORLD );
165 #endif
166 
167       return newgrid;
168     }
169 
170     template <class FromGrid, class Vector>
convert(const FromGrid & fromGrid,Vector & cellData)171     GridPtrType convert( const FromGrid& fromGrid, Vector& cellData )
172     {
173       return convert( fromGrid, cellData, ordering_ );
174     }
175 
176     template <class FromGrid>
convert(const FromGrid & fromGrid)177     GridPtrType convert( const FromGrid& fromGrid )
178     {
179       std::vector<int> dummy(0);
180       return convert( fromGrid, dummy, ordering_ );
181     }
182   protected:
183     template <int codim, class Entity>
subEntities(const Entity & entity) const184     int subEntities ( const Entity& entity ) const
185     {
186       return entity.subEntities( codim );
187     }
188   };
189 
190 } // namespace Dune
191 
192 #endif // #ifndef DUNE_ALUGRID_STRUCTUREDGRIDFACTORY_HH
193