3 #ifndef DUNE_STRUCTURED_GRID_FACTORY_HH
4 #define DUNE_STRUCTURED_GRID_FACTORY_HH
16 #include <dune/common/classname.hh>
17 #include <dune/common/exceptions.hh>
18 #include <dune/common/fvector.hh>
19 #include <dune/common/parallel/mpihelper.hh>
28 template <
class Gr
idType>
31 typedef typename GridType::ctype ctype;
33 static const int dim = GridType::dimension;
35 static const int dimworld = GridType::dimensionworld;
39 const FieldVector<ctype,dimworld>& lowerLeft,
40 const FieldVector<ctype,dimworld>& upperRight,
41 const std::array<unsigned int,dim>& vertices)
46 int numVertices = index.
cycle();
49 for (
int i=0; i<numVertices; i++, ++index) {
52 FieldVector<double,dimworld> pos(0);
53 for (
int j=0; j<dim; j++)
54 pos[j] = lowerLeft[j] + index[j] * (upperRight[j]-lowerLeft[j])/(vertices[j]-1);
55 for (
int j=dim; j<dimworld; j++)
56 pos[j] = lowerLeft[j];
66 static std::array<unsigned int, dim> computeUnitOffsets(
const std::array<unsigned int,dim>& vertices)
68 std::array<unsigned int, dim> unitOffsets;
72 for (
int i=1; i<dim; i++)
73 unitOffsets[i] = unitOffsets[i-1] * vertices[i-1];
89 static std::shared_ptr<GridType>
createCubeGrid(
const FieldVector<ctype,dimworld>& lowerLeft,
90 const FieldVector<ctype,dimworld>& upperRight,
91 const std::array<unsigned int,dim>& elements)
96 if (MPIHelper::getCollectiveCommunication().rank() == 0)
99 std::array<unsigned int,dim> vertices = elements;
100 for(
size_t i = 0; i < vertices.size(); ++i )
104 insertVertices(factory, lowerLeft, upperRight, vertices);
108 std::array<unsigned int, dim> unitOffsets =
109 computeUnitOffsets(vertices);
113 unsigned int nCorners = 1<<dim;
115 std::vector<unsigned int> cornersTemplate(nCorners,0);
117 for (
size_t i=0; i<nCorners; i++)
118 for (
int j=0; j<dim; j++)
120 cornersTemplate[i] += unitOffsets[j];
126 int numElements = index.
cycle();
128 for (
int i=0; i<numElements; i++, ++index) {
131 unsigned int base = 0;
132 for (
int j=0; j<dim; j++)
133 base += index[j] * unitOffsets[j];
136 std::vector<unsigned int> corners = cornersTemplate;
137 for (
size_t j=0; j<corners.size(); j++)
148 return std::shared_ptr<GridType>(factory.
createGrid());
165 static std::shared_ptr<GridType>
createSimplexGrid(
const FieldVector<ctype,dimworld>& lowerLeft,
166 const FieldVector<ctype,dimworld>& upperRight,
167 const std::array<unsigned int,dim>& elements)
172 if(MPIHelper::getCollectiveCommunication().rank() == 0)
175 std::array<unsigned int,dim> vertices = elements;
176 for (std::size_t i=0; i<vertices.size(); i++)
179 insertVertices(factory, lowerLeft, upperRight, vertices);
183 std::array<unsigned int, dim> unitOffsets =
184 computeUnitOffsets(vertices);
189 size_t cycle = elementsIndex.
cycle();
191 for (
size_t i=0; i<cycle; ++elementsIndex, i++) {
194 unsigned int base = 0;
195 for (
int j=0; j<dim; j++)
196 base += elementsIndex[j] * unitOffsets[j];
199 std::vector<unsigned int> permutation(dim);
200 for (
int j=0; j<dim; j++)
206 std::vector<unsigned int> corners(dim+1);
209 for (
int j=0; j<dim; j++)
211 corners[j] + unitOffsets[permutation[j]];
217 }
while (std::next_permutation(permutation.begin(),
225 return std::shared_ptr<GridType>(factory.
createGrid());