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::unique_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++)
147 return std::unique_ptr<GridType>(factory.
createGrid());
164 static std::unique_ptr<GridType>
createSimplexGrid(
const FieldVector<ctype,dimworld>& lowerLeft,
165 const FieldVector<ctype,dimworld>& upperRight,
166 const std::array<unsigned int,dim>& elements)
171 if(MPIHelper::getCollectiveCommunication().rank() == 0)
174 std::array<unsigned int,dim> vertices = elements;
175 for (std::size_t i=0; i<vertices.size(); i++)
178 insertVertices(factory, lowerLeft, upperRight, vertices);
182 std::array<unsigned int, dim> unitOffsets =
183 computeUnitOffsets(vertices);
188 size_t cycle = elementsIndex.
cycle();
190 for (
size_t i=0; i<cycle; ++elementsIndex, i++) {
193 unsigned int base = 0;
194 for (
int j=0; j<dim; j++)
195 base += elementsIndex[j] * unitOffsets[j];
198 std::vector<unsigned int> permutation(dim);
199 for (
int j=0; j<dim; j++)
205 std::vector<unsigned int> corners(dim+1);
208 for (
int j=0; j<dim; j++)
210 corners[j] + unitOffsets[permutation[j]];
214 }
while (std::next_permutation(permutation.begin(),
222 return std::unique_ptr<GridType>(factory.
createGrid());