dune-grid
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- 00002 // vi: set et ts=4 sw=4 sts=4: 00003 #ifndef DUNE_STRUCTURED_GRID_FACTORY_HH 00004 #define DUNE_STRUCTURED_GRID_FACTORY_HH 00005 00010 #include <algorithm> 00011 #include <cstddef> 00012 #include <cstdlib> 00013 00014 #include <dune/common/array.hh> 00015 #include <dune/common/classname.hh> 00016 #include <dune/common/exceptions.hh> 00017 #include <dune/common/fvector.hh> 00018 #include <dune/common/mpihelper.hh> 00019 #include <dune/common/shared_ptr.hh> 00020 00021 #include <dune/grid/common/gridfactory.hh> 00022 #include <dune/grid/yaspgrid.hh> 00023 #include <dune/grid/sgrid.hh> 00024 00025 namespace Dune { 00026 00029 template <class GridType> 00030 class StructuredGridFactory 00031 { 00032 typedef typename GridType::ctype ctype; 00033 00034 static const int dim = GridType::dimension; 00035 00036 static const int dimworld = GridType::dimensionworld; 00037 00040 class MultiIndex 00041 : public array<unsigned int,dim> 00042 { 00043 00044 // The range of each component 00045 array<unsigned int,dim> limits_; 00046 00047 public: 00049 MultiIndex(const array<unsigned int,dim>& limits) 00050 : limits_(limits) 00051 { 00052 std::fill(this->begin(), this->end(), 0); 00053 } 00054 00056 MultiIndex& operator++() { 00057 00058 for (int i=0; i<dim; i++) { 00059 00060 // Augment digit 00061 (*this)[i]++; 00062 00063 // If there is no carry-over we can stop here 00064 if ((*this)[i]<limits_[i]) 00065 break; 00066 00067 (*this)[i] = 0; 00068 00069 } 00070 return *this; 00071 } 00072 00074 size_t cycle() const { 00075 size_t result = 1; 00076 for (int i=0; i<dim; i++) 00077 result *= limits_[i]; 00078 return result; 00079 } 00080 00081 }; 00082 00084 static void insertVertices(GridFactory<GridType>& factory, 00085 const FieldVector<ctype,dimworld>& lowerLeft, 00086 const FieldVector<ctype,dimworld>& upperRight, 00087 const array<unsigned int,dim>& vertices) 00088 { 00089 00090 MultiIndex index(vertices); 00091 00092 // Compute the total number of vertices to be created 00093 int numVertices = index.cycle(); 00094 00095 // Create vertices 00096 for (int i=0; i<numVertices; i++, ++index) { 00097 00098 // scale the multiindex to obtain a world position 00099 FieldVector<double,dimworld> pos(0); 00100 for (int j=0; j<dimworld; j++) 00101 pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1); 00102 00103 factory.insertVertex(pos); 00104 00105 } 00106 00107 } 00108 00109 // Compute the index offsets needed to move to the adjacent vertices 00110 // in the different coordinate directions 00111 static array<unsigned int, dim> computeUnitOffsets(const array<unsigned int,dim>& vertices) 00112 { 00113 array<unsigned int, dim> unitOffsets; 00114 if (dim>0) // paranoia 00115 unitOffsets[0] = 1; 00116 00117 for (int i=1; i<dim; i++) 00118 unitOffsets[i] = unitOffsets[i-1] * vertices[i-1]; 00119 00120 return unitOffsets; 00121 } 00122 00123 public: 00124 00130 static shared_ptr<GridType> createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00131 const FieldVector<ctype,dimworld>& upperRight, 00132 const array<unsigned int,dim>& elements) 00133 { 00134 // The grid factory 00135 GridFactory<GridType> factory; 00136 00137 if (MPIHelper::getCollectiveCommunication().rank() == 0) 00138 { 00139 // Insert uniformly spaced vertices 00140 array<unsigned int,dim> vertices = elements; 00141 for( size_t i = 0; i < vertices.size(); ++i ) 00142 vertices[i]++; 00143 00144 // Insert vertices for structured grid into the factory 00145 insertVertices(factory, lowerLeft, upperRight, vertices); 00146 00147 // Compute the index offsets needed to move to the adjacent 00148 // vertices in the different coordinate directions 00149 array<unsigned int, dim> unitOffsets = 00150 computeUnitOffsets(vertices); 00151 00152 // Compute an element template (the cube at (0,...,0). All 00153 // other cubes are constructed by moving this template around 00154 unsigned int nCorners = 1<<dim; 00155 00156 std::vector<unsigned int> cornersTemplate(nCorners,0); 00157 00158 for (size_t i=0; i<nCorners; i++) 00159 for (int j=0; j<dim; j++) 00160 if ( i & (1<<j) ) 00161 cornersTemplate[i] += unitOffsets[j]; 00162 00163 // Insert elements 00164 MultiIndex index(elements); 00165 00166 // Compute the total number of elementss to be created 00167 int numElements = index.cycle(); 00168 00169 for (int i=0; i<numElements; i++, ++index) { 00170 00171 // 'base' is the index of the lower left element corner 00172 unsigned int base = 0; 00173 for (int j=0; j<dim; j++) 00174 base += index[j] * unitOffsets[j]; 00175 00176 // insert new element 00177 std::vector<unsigned int> corners = cornersTemplate; 00178 for (size_t j=0; j<corners.size(); j++) 00179 corners[j] += base; 00180 00181 factory.insertElement 00182 (GeometryType(GeometryType::cube, dim), corners); 00183 00184 } 00185 00186 } // if(rank == 0) 00187 00188 // Create the grid and hand it to the calling method 00189 return shared_ptr<GridType>(factory.createGrid()); 00190 00191 } 00192 00199 static shared_ptr<GridType> createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00200 const FieldVector<ctype,dimworld>& upperRight, 00201 const array<unsigned int,dim>& elements) 00202 { 00203 // The grid factory 00204 GridFactory<GridType> factory; 00205 00206 if(MPIHelper::getCollectiveCommunication().rank() == 0) 00207 { 00208 // Insert uniformly spaced vertices 00209 array<unsigned int,dim> vertices = elements; 00210 for (std::size_t i=0; i<vertices.size(); i++) 00211 vertices[i]++; 00212 00213 insertVertices(factory, lowerLeft, upperRight, vertices); 00214 00215 // Compute the index offsets needed to move to the adjacent 00216 // vertices in the different coordinate directions 00217 array<unsigned int, dim> unitOffsets = 00218 computeUnitOffsets(vertices); 00219 00220 // Insert the elements 00221 std::vector<unsigned int> corners(dim+1); 00222 00223 // Loop over all "cubes", and split up each cube into dim! 00224 // (factorial) simplices 00225 MultiIndex elementsIndex(elements); 00226 size_t cycle = elementsIndex.cycle(); 00227 00228 for (size_t i=0; i<cycle; ++elementsIndex, i++) { 00229 00230 // 'base' is the index of the lower left element corner 00231 unsigned int base = 0; 00232 for (int j=0; j<dim; j++) 00233 base += elementsIndex[j] * unitOffsets[j]; 00234 00235 // each permutation of the unit vectors gives a simplex. 00236 std::vector<unsigned int> permutation(dim); 00237 for (int j=0; j<dim; j++) 00238 permutation[j] = j; 00239 00240 do { 00241 00242 // Make a simplex 00243 std::vector<unsigned int> corners(dim+1); 00244 corners[0] = base; 00245 00246 for (int j=0; j<dim; j++) 00247 corners[j+1] = 00248 corners[j] + unitOffsets[permutation[j]]; 00249 00250 factory.insertElement 00251 (GeometryType(GeometryType::simplex, dim), 00252 corners); 00253 00254 } while (std::next_permutation(permutation.begin(), 00255 permutation.end())); 00256 00257 } 00258 00259 } // if(rank == 0) 00260 00261 // Create the grid and hand it to the calling method 00262 return shared_ptr<GridType>(factory.createGrid()); 00263 } 00264 00265 }; 00266 00276 template<int dim> 00277 class StructuredGridFactory<YaspGrid<dim> > { 00278 typedef YaspGrid<dim> GridType; 00279 typedef typename GridType::ctype ctype; 00280 static const int dimworld = GridType::dimensionworld; 00281 00282 public: 00292 static shared_ptr<GridType> 00293 createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00294 const FieldVector<ctype,dimworld>& upperRight, 00295 const array<unsigned int,dim>& elements) 00296 { 00297 for(int d = 0; d < dimworld; ++d) 00298 if(std::abs(lowerLeft[d]) > std::abs(upperRight[d])*1e-10) 00299 DUNE_THROW(GridError, className<StructuredGridFactory>() 00300 << "::createCubeGrid(): The lower coordinates " 00301 "must be at the origin for YaspGrid."); 00302 00303 FieldVector<int, dim> elements_; 00304 std::copy(elements.begin(), elements.end(), elements_.begin()); 00305 00306 return shared_ptr<GridType> 00307 (new GridType(upperRight, elements_, 00308 FieldVector<bool,dim>(false), 0)); 00309 } 00310 00316 static shared_ptr<GridType> 00317 createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00318 const FieldVector<ctype,dimworld>& upperRight, 00319 const array<unsigned int,dim>& elements) 00320 { 00321 DUNE_THROW(GridError, className<StructuredGridFactory>() 00322 << "::createSimplexGrid(): Simplices are not supported " 00323 "by YaspGrid."); 00324 } 00325 00326 }; 00327 00334 template<int dim> 00335 class StructuredGridFactory<SGrid<dim, dim> > { 00336 typedef SGrid<dim, dim> GridType; 00337 typedef typename GridType::ctype ctype; 00338 static const int dimworld = GridType::dimensionworld; 00339 00340 public: 00347 static shared_ptr<GridType> 00348 createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00349 const FieldVector<ctype,dimworld>& upperRight, 00350 const array<unsigned int,dim>& elements) 00351 { 00352 FieldVector<int, dim> elements_; 00353 std::copy(elements.begin(), elements.end(), elements_.begin()); 00354 00355 return shared_ptr<GridType> 00356 (new GridType(elements_, lowerLeft, upperRight)); 00357 } 00358 00368 static shared_ptr<GridType> 00369 createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft, 00370 const FieldVector<ctype,dimworld>& upperRight, 00371 const array<unsigned int,dim>& elements) 00372 { 00373 DUNE_THROW(GridError, className<StructuredGridFactory>() 00374 << "::createSimplexGrid(): Simplices are not supported " 00375 "by SGrid."); 00376 } 00377 }; 00378 00379 } // namespace Dune 00380 00381 #endif