3 #ifndef DUNE_STRUCTURED_GRID_FACTORY_HH
4 #define DUNE_STRUCTURED_GRID_FACTORY_HH
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/mpihelper.hh>
19 #include <dune/common/shared_ptr.hh>
29 template <
class Gr
idType>
32 typedef typename GridType::ctype ctype;
34 static const int dim = GridType::dimension;
36 static const int dimworld = GridType::dimensionworld;
41 :
public array<unsigned int,dim>
45 array<unsigned int,dim> limits_;
49 MultiIndex(
const array<unsigned int,dim>& limits)
52 std::fill(this->begin(), this->end(), 0);
56 MultiIndex& operator++() {
58 for (
int i=0; i<dim; i++) {
64 if ((*
this)[i]<limits_[i])
74 size_t cycle()
const {
76 for (
int i=0; i<dim; i++)
85 const FieldVector<ctype,dimworld>& lowerLeft,
86 const FieldVector<ctype,dimworld>& upperRight,
87 const array<unsigned int,dim>& vertices)
90 MultiIndex index(vertices);
93 int numVertices = index.cycle();
96 for (
int i=0; i<numVertices; i++, ++index) {
99 FieldVector<double,dimworld> pos(0);
100 for (
int j=0; j<dimworld; j++)
101 pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
111 static array<unsigned int, dim> computeUnitOffsets(
const array<unsigned int,dim>& vertices)
113 array<unsigned int, dim> unitOffsets;
117 for (
int i=1; i<dim; i++)
118 unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
130 static shared_ptr<GridType>
createCubeGrid(
const FieldVector<ctype,dimworld>& lowerLeft,
131 const FieldVector<ctype,dimworld>& upperRight,
132 const array<unsigned int,dim>& elements)
137 if (MPIHelper::getCollectiveCommunication().rank() == 0)
140 array<unsigned int,dim> vertices = elements;
141 for(
size_t i = 0; i < vertices.size(); ++i )
145 insertVertices(factory, lowerLeft, upperRight, vertices);
149 array<unsigned int, dim> unitOffsets =
150 computeUnitOffsets(vertices);
154 unsigned int nCorners = 1<<dim;
156 std::vector<unsigned int> cornersTemplate(nCorners,0);
158 for (
size_t i=0; i<nCorners; i++)
159 for (
int j=0; j<dim; j++)
161 cornersTemplate[i] += unitOffsets[j];
164 MultiIndex index(elements);
167 int numElements = index.cycle();
169 for (
int i=0; i<numElements; i++, ++index) {
172 unsigned int base = 0;
173 for (
int j=0; j<dim; j++)
174 base += index[j] * unitOffsets[j];
177 std::vector<unsigned int> corners = cornersTemplate;
178 for (
size_t j=0; j<corners.size(); j++)
189 return shared_ptr<GridType>(factory.
createGrid());
200 const FieldVector<ctype,dimworld>& upperRight,
201 const array<unsigned int,dim>& elements)
206 if(MPIHelper::getCollectiveCommunication().rank() == 0)
209 array<unsigned int,dim> vertices = elements;
210 for (std::size_t i=0; i<vertices.size(); i++)
213 insertVertices(factory, lowerLeft, upperRight, vertices);
217 array<unsigned int, dim> unitOffsets =
218 computeUnitOffsets(vertices);
221 std::vector<unsigned int> corners(dim+1);
225 MultiIndex elementsIndex(elements);
226 size_t cycle = elementsIndex.cycle();
228 for (
size_t i=0; i<cycle; ++elementsIndex, i++) {
231 unsigned int base = 0;
232 for (
int j=0; j<dim; j++)
233 base += elementsIndex[j] * unitOffsets[j];
236 std::vector<unsigned int> permutation(dim);
237 for (
int j=0; j<dim; j++)
243 std::vector<unsigned int> corners(dim+1);
246 for (
int j=0; j<dim; j++)
248 corners[j] + unitOffsets[permutation[j]];
254 }
while (std::next_permutation(permutation.begin(),
262 return shared_ptr<GridType>(factory.
createGrid());
280 static const int dimworld = GridType::dimensionworld;
292 static shared_ptr<GridType>
294 const FieldVector<ctype,dimworld>& upperRight,
295 const array<unsigned int,dim>& elements)
297 for(
int d = 0; d < dimworld; ++d)
299 DUNE_THROW(
GridError, className<StructuredGridFactory>()
300 <<
"::createCubeGrid(): The lower coordinates "
301 "must be at the origin for YaspGrid.");
303 FieldVector<int, dim> elements_;
304 std::copy(elements.begin(), elements.end(), elements_.begin());
306 return shared_ptr<GridType>
307 (
new GridType(upperRight, elements_,
308 FieldVector<bool,dim>(
false), 0));
316 static shared_ptr<GridType>
318 const FieldVector<ctype,dimworld>& upperRight,
319 const array<unsigned int,dim>& elements)
321 DUNE_THROW(
GridError, className<StructuredGridFactory>()
322 <<
"::createSimplexGrid(): Simplices are not supported "
338 static const int dimworld = GridType::dimensionworld;
347 static shared_ptr<GridType>
349 const FieldVector<ctype,dimworld>& upperRight,
350 const array<unsigned int,dim>& elements)
352 FieldVector<int, dim> elements_;
353 std::copy(elements.begin(), elements.end(), elements_.begin());
355 return shared_ptr<GridType>
356 (
new GridType(elements_, lowerLeft, upperRight));
368 static shared_ptr<GridType>
370 const FieldVector<ctype,dimworld>& upperRight,
371 const array<unsigned int,dim>& elements)
373 DUNE_THROW(
GridError, className<StructuredGridFactory>()
374 <<
"::createSimplexGrid(): Simplices are not supported "