dune-grid  2.3.1
structuredgridfactory.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_STRUCTURED_GRID_FACTORY_HH
4 #define DUNE_STRUCTURED_GRID_FACTORY_HH
5 
10 #include <algorithm>
11 #include <cstddef>
12 #include <cstdlib>
13 
14 #include <dune/common/array.hh>
15 #include <dune/common/classname.hh>
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/parallel/mpihelper.hh>
19 #include <dune/common/shared_ptr.hh>
20 
22 #include <dune/grid/yaspgrid.hh>
23 #include <dune/grid/sgrid.hh>
24 
25 namespace Dune {
26 
29  template <class GridType>
31  {
32  typedef typename GridType::ctype ctype;
33 
34  static const int dim = GridType::dimension;
35 
36  static const int dimworld = GridType::dimensionworld;
37 
40  class MultiIndex
41  : public array<unsigned int,dim>
42  {
43 
44  // The range of each component
45  array<unsigned int,dim> limits_;
46 
47  public:
49  MultiIndex(const array<unsigned int,dim>& limits)
50  : limits_(limits)
51  {
52  std::fill(this->begin(), this->end(), 0);
53  }
54 
56  MultiIndex& operator++() {
57 
58  for (int i=0; i<dim; i++) {
59 
60  // Augment digit
61  (*this)[i]++;
62 
63  // If there is no carry-over we can stop here
64  if ((*this)[i]<limits_[i])
65  break;
66 
67  (*this)[i] = 0;
68 
69  }
70  return *this;
71  }
72 
74  size_t cycle() const {
75  size_t result = 1;
76  for (int i=0; i<dim; i++)
77  result *= limits_[i];
78  return result;
79  }
80 
81  };
82 
84  static void insertVertices(GridFactory<GridType>& factory,
85  const FieldVector<ctype,dimworld>& lowerLeft,
86  const FieldVector<ctype,dimworld>& upperRight,
87  const array<unsigned int,dim>& vertices)
88  {
89 
90  MultiIndex index(vertices);
91 
92  // Compute the total number of vertices to be created
93  int numVertices = index.cycle();
94 
95  // Create vertices
96  for (int i=0; i<numVertices; i++, ++index) {
97 
98  // scale the multiindex to obtain a world position
99  FieldVector<double,dimworld> pos(0);
100  for (int j=0; j<dim; j++)
101  pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
102  for (int j=dim; j<dimworld; j++)
103  pos[j] = lowerLeft[j];
104 
105  factory.insertVertex(pos);
106 
107  }
108 
109  }
110 
111  // Compute the index offsets needed to move to the adjacent vertices
112  // in the different coordinate directions
113  static array<unsigned int, dim> computeUnitOffsets(const array<unsigned int,dim>& vertices)
114  {
115  array<unsigned int, dim> unitOffsets;
116  if (dim>0) // paranoia
117  unitOffsets[0] = 1;
118 
119  for (int i=1; i<dim; i++)
120  unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
121 
122  return unitOffsets;
123  }
124 
125  public:
126 
136  static shared_ptr<GridType> createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
137  const FieldVector<ctype,dimworld>& upperRight,
138  const array<unsigned int,dim>& elements)
139  {
140  // The grid factory
141  GridFactory<GridType> factory;
142 
143  if (MPIHelper::getCollectiveCommunication().rank() == 0)
144  {
145  // Insert uniformly spaced vertices
146  array<unsigned int,dim> vertices = elements;
147  for( size_t i = 0; i < vertices.size(); ++i )
148  vertices[i]++;
149 
150  // Insert vertices for structured grid into the factory
151  insertVertices(factory, lowerLeft, upperRight, vertices);
152 
153  // Compute the index offsets needed to move to the adjacent
154  // vertices in the different coordinate directions
155  array<unsigned int, dim> unitOffsets =
156  computeUnitOffsets(vertices);
157 
158  // Compute an element template (the cube at (0,...,0). All
159  // other cubes are constructed by moving this template around
160  unsigned int nCorners = 1<<dim;
161 
162  std::vector<unsigned int> cornersTemplate(nCorners,0);
163 
164  for (size_t i=0; i<nCorners; i++)
165  for (int j=0; j<dim; j++)
166  if ( i & (1<<j) )
167  cornersTemplate[i] += unitOffsets[j];
168 
169  // Insert elements
170  MultiIndex index(elements);
171 
172  // Compute the total number of elementss to be created
173  int numElements = index.cycle();
174 
175  for (int i=0; i<numElements; i++, ++index) {
176 
177  // 'base' is the index of the lower left element corner
178  unsigned int base = 0;
179  for (int j=0; j<dim; j++)
180  base += index[j] * unitOffsets[j];
181 
182  // insert new element
183  std::vector<unsigned int> corners = cornersTemplate;
184  for (size_t j=0; j<corners.size(); j++)
185  corners[j] += base;
186 
187  factory.insertElement
188  (GeometryType(GeometryType::cube, dim), corners);
189 
190  }
191 
192  } // if(rank == 0)
193 
194  // Create the grid and hand it to the calling method
195  return shared_ptr<GridType>(factory.createGrid());
196 
197  }
198 
212  static shared_ptr<GridType> createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
213  const FieldVector<ctype,dimworld>& upperRight,
214  const array<unsigned int,dim>& elements)
215  {
216  // The grid factory
217  GridFactory<GridType> factory;
218 
219  if(MPIHelper::getCollectiveCommunication().rank() == 0)
220  {
221  // Insert uniformly spaced vertices
222  array<unsigned int,dim> vertices = elements;
223  for (std::size_t i=0; i<vertices.size(); i++)
224  vertices[i]++;
225 
226  insertVertices(factory, lowerLeft, upperRight, vertices);
227 
228  // Compute the index offsets needed to move to the adjacent
229  // vertices in the different coordinate directions
230  array<unsigned int, dim> unitOffsets =
231  computeUnitOffsets(vertices);
232 
233  // Insert the elements
234  std::vector<unsigned int> corners(dim+1);
235 
236  // Loop over all "cubes", and split up each cube into dim!
237  // (factorial) simplices
238  MultiIndex elementsIndex(elements);
239  size_t cycle = elementsIndex.cycle();
240 
241  for (size_t i=0; i<cycle; ++elementsIndex, i++) {
242 
243  // 'base' is the index of the lower left element corner
244  unsigned int base = 0;
245  for (int j=0; j<dim; j++)
246  base += elementsIndex[j] * unitOffsets[j];
247 
248  // each permutation of the unit vectors gives a simplex.
249  std::vector<unsigned int> permutation(dim);
250  for (int j=0; j<dim; j++)
251  permutation[j] = j;
252 
253  do {
254 
255  // Make a simplex
256  std::vector<unsigned int> corners(dim+1);
257  corners[0] = base;
258 
259  for (int j=0; j<dim; j++)
260  corners[j+1] =
261  corners[j] + unitOffsets[permutation[j]];
262 
263  factory.insertElement
265  corners);
266 
267  } while (std::next_permutation(permutation.begin(),
268  permutation.end()));
269 
270  }
271 
272  } // if(rank == 0)
273 
274  // Create the grid and hand it to the calling method
275  return shared_ptr<GridType>(factory.createGrid());
276  }
277 
278  };
279 
289  template<int dim>
291  typedef YaspGrid<dim> GridType;
292  typedef typename GridType::ctype ctype;
293  static const int dimworld = GridType::dimensionworld;
294 
295  public:
305  static shared_ptr<GridType>
306  createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
307  const FieldVector<ctype,dimworld>& upperRight,
308  const array<unsigned int,dim>& elements)
309  {
310  for(int d = 0; d < dimworld; ++d)
311  if(std::abs(lowerLeft[d]) > std::abs(upperRight[d])*1e-10)
312  DUNE_THROW(GridError, className<StructuredGridFactory>()
313  << "::createCubeGrid(): The lower coordinates "
314  "must be at the origin for YaspGrid.");
315 
316  Dune::array<int, dim> elements_;
317  std::copy(elements.begin(), elements.end(), elements_.begin());
318 
319  return shared_ptr<GridType>
320  (new GridType(upperRight, elements_,
321  std::bitset<dim>(), 0)); // default constructor of bitset sets to zero
322  }
323 
329  static shared_ptr<GridType>
330  createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
331  const FieldVector<ctype,dimworld>& upperRight,
332  const array<unsigned int,dim>& elements)
333  {
334  DUNE_THROW(GridError, className<StructuredGridFactory>()
335  << "::createSimplexGrid(): Simplices are not supported "
336  "by YaspGrid.");
337  }
338 
339  };
340 
347  template<int dim>
348  class StructuredGridFactory<SGrid<dim, dim> > {
349  typedef SGrid<dim, dim> GridType;
350  typedef typename GridType::ctype ctype;
351  static const int dimworld = GridType::dimensionworld;
352 
353  public:
360  static shared_ptr<GridType>
361  createCubeGrid(const FieldVector<ctype,dimworld>& lowerLeft,
362  const FieldVector<ctype,dimworld>& upperRight,
363  const array<unsigned int,dim>& elements)
364  {
365  FieldVector<int, dim> elements_;
366  std::copy(elements.begin(), elements.end(), elements_.begin());
367 
368  return shared_ptr<GridType>
369  (new GridType(elements_, lowerLeft, upperRight));
370  }
371 
381  static shared_ptr<GridType>
382  createSimplexGrid(const FieldVector<ctype,dimworld>& lowerLeft,
383  const FieldVector<ctype,dimworld>& upperRight,
384  const array<unsigned int,dim>& elements)
385  {
386  DUNE_THROW(GridError, className<StructuredGridFactory>()
387  << "::createSimplexGrid(): Simplices are not supported "
388  "by SGrid.");
389  }
390  };
391 
392 } // namespace Dune
393 
394 #endif
static shared_ptr< GridType > createCubeGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured cube grid.
Definition: structuredgridfactory.hh:306
static shared_ptr< GridType > createCubeGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured cube grid.
Definition: structuredgridfactory.hh:361
static shared_ptr< GridType > createSimplexGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured simplex grid.
Definition: structuredgridfactory.hh:212
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
[ provides Dune::Grid ]
Definition: yaspgrid.hh:62
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: common/gridfactory.hh:290
Provide a generic factory class for unstructured grids.
Definition: common/gridfactory.hh:263
Definition: alugrid/common/declaration.hh:18
Provide a generic factory class for unstructured grids.
[ provides Dune::Grid ]
Definition: sgrid.hh:48
virtual GridType * createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:316
Base class for exceptions in Dune grid modules.
Definition: exceptions.hh:16
Construct structured cube and simplex grids in unstructured grid managers.
Definition: structuredgridfactory.hh:30
_ctype ctype
define type used for coordinates in grid module
Definition: sgrid.hh:1284
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: common/gridfactory.hh:279
Definition: alugrid/common/declaration.hh:18
static shared_ptr< GridType > createSimplexGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured simplex grid.
Definition: structuredgridfactory.hh:330
static shared_ptr< GridType > createSimplexGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured simplex grid.
Definition: structuredgridfactory.hh:382
yaspgrid_ctype ctype
Type used for coordinates.
Definition: yaspgrid.hh:196
static shared_ptr< GridType > createCubeGrid(const FieldVector< ctype, dimworld > &lowerLeft, const FieldVector< ctype, dimworld > &upperRight, const array< unsigned int, dim > &elements)
Create a structured cube grid.
Definition: structuredgridfactory.hh:136
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:332