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