dune-grid  2.2.0
structuredgridfactory.hh
Go to the documentation of this file.
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