3#ifndef DUNE_GRID_UTILITY_TENSORGRIDFACTORY_HH
4#define DUNE_GRID_UTILITY_TENSORGRIDFACTORY_HH
33 template<
typename Gr
id>
34 class TensorGridFactoryCreator;
40 template<
typename Gr
id>
44 typedef typename Grid::Traits::Communication Comm;
48 std::unique_ptr<Grid> createGrid(Comm comm = Comm())
50 TensorGridFactoryCreator<Grid> creator(*
this);
51 return creator.createGrid(comm);
54 std::array<std::vector<ctype> , dim> coords()
const
66 const std::vector<ctype>&
operator[](std::size_t d)
const
82 _coords[d][0] = value;
96 for (
int i = 0; i < n; i++)
111 const ctype h = (end - _coords[d].back ()) / n;
112 for (
int i = 0; i < n - 1; i++)
114 _coords[d].push_back (end);
128 while (_coords[d].
back () < end)
129 _coords[d].push_back (_coords[d].
back () + h);
145 static_cast<ctype
> (0))
149 if (h0 ==
static_cast<ctype
>(0))
150 h = lastInterval (d) * ratio;
151 for (
int i = 0; i < n; i++)
153 _coords[d].push_back (_coords[d].
back () + h);
174 if (h0 ==
static_cast<ctype
>(0))
175 h = lastInterval (d) * ratio;
176 while (_coords[d].
back () < end)
178 _coords[d].push_back (_coords[d].
back () + h);
198 static_cast<ctype
> (0),
203 h = lastInterval (d);
204 ctype ratio = newton (n, _coords[d].
back (), end, h);
207 h = h * pow (ratio, n - 1);
210 for (
int i = 0; i < n - 1; i++)
212 _coords[d].push_back (_coords[d].
back () + h);
215 _coords[d].push_back (end);
221 for (
int i=0; i<dim; i++)
223 std::cout <<
"Container in direction " << i <<
":" << std::endl <<
"Coordinates: ";
224 for (
auto it = _coords[i].begin(); it != _coords[i].end(); ++it)
225 std::cout << *it <<
" ";
226 std::cout << std::endl <<
"Interval lengths: ";
228 std::vector<ctype> meshsize;
229 for (
auto it = _coords[i].begin(); it != _coords[i].end()-1;)
231 meshsize.push_back(-1.*(*it));
233 meshsize.back() += *it;
236 for (
auto it = meshsize.begin(); it != meshsize.end(); ++it)
237 std::cout << *it <<
" ";
238 std::cout << std::endl <<
"Ratios between interval lengths: ";
240 std::vector<ctype> ratios;
241 for (
auto it = meshsize.begin(); it != meshsize.end() - 1 ;)
242 ratios.push_back((1./(*it)) * *(++it));
244 for (
auto it = ratios.begin(); it != ratios.end(); ++it)
245 std::cout << *it <<
" ";
246 std::cout << std::endl << std::endl << std::endl;
252 void emptyCheck (
int i)
254 if (_coords[i].
empty ())
255 _coords[i].push_back (
static_cast<ctype
> (0));
259 ctype lastInterval (
int d)
261 if (_coords[d].
size () < 2)
264 "Not enough elements in coordinate container to deduce interval length in TensorYaspFactory");
266 return _coords[d].back () - _coords[d][_coords[d].size () - 2];
272 ctype newton (
int n, ctype x_s, ctype x_e, ctype h)
274 ctype m = (x_e - x_s) / h;
276 ctype xnew = x_e - x_s;
277 while (std::abs (xnew - xold) > 1E-8)
281 - (-pow (xold, n) + m * xold - m + 1)
282 / (-n * pow (xold, n - 1) + m);
284 if (std::abs (xnew - 1) < 1E-6)
288 while (std::abs (xnew - xold) > 1E-8)
292 - (-pow (xold, n) + m * xold - m + 1)
293 / (-n * pow (xold, n - 1) + m);
299 std::array<std::vector<ctype>, dim> _coords;
304 template<
typename Gr
id>
305 class TensorGridFactoryCreator
308 typedef typename Grid::Traits::Communication Comm;
312 TensorGridFactoryCreator(
const TensorGridFactory<Grid>& factory) : _factory(factory) {}
314 std::unique_ptr<Grid> createGrid(Comm comm)
317 GridFactory<Grid> fac;
319 if (comm.rank() == 0)
322 std::array<unsigned int, dim> vsizes, esizes;
323 std::size_t
size = 1;
324 for (std::size_t i = 0; i<dim; ++i)
326 vsizes[i] = _factory[i].size();
327 esizes[i] = vsizes[i] - 1;
332 FactoryUtilities::MultiIndex<dim> index(vsizes);
333 for (std::size_t i=0; i<
size; ++i, ++index)
336 for (std::size_t j = 0; j<dim; ++j)
337 position[j] = _factory[j][index[j]];
338 fac.insertVertex(position);
342 std::array<std::size_t, dim> offsets;
344 for (std::size_t i=1; i<dim; i++)
345 offsets[i] = offsets[i-1] * vsizes[i-1];
349 unsigned int nCorners = 1<<dim;
351 std::vector<unsigned int> cornersTemplate(nCorners,0);
353 for (
size_t i=0; i<nCorners; i++)
354 for (
int j=0; j<dim; j++)
356 cornersTemplate[i] += offsets[j];
359 FactoryUtilities::MultiIndex<dim> eindex(esizes);
362 int numElements = eindex.cycle();
364 for (
int i=0; i<numElements; i++, ++eindex)
367 unsigned int base = 0;
368 for (
int j=0; j<dim; j++)
369 base += eindex[j] * offsets[j];
372 std::vector<unsigned int> corners = cornersTemplate;
373 for (
size_t j=0; j<corners.size(); j++)
380 return std::unique_ptr<Grid>(fac.createGrid());
384 const TensorGridFactory<Grid>& _factory;
387 template<
typename ctype,
int dim>
388 class TensorGridFactoryCreator<YaspGrid<dim, TensorProductCoordinates<ctype, dim> > >
391 typedef YaspGrid<dim, TensorProductCoordinates<ctype, dim> > Grid;
394 TensorGridFactoryCreator(
const TensorGridFactory<Grid>& factory) : _factory(factory) {}
396 std::unique_ptr<Grid> createGrid(Comm comm)
398 return std::make_unique<Grid>(_factory.coords(), std::bitset<dim>(0ULL), 1, comm);
401 const TensorGridFactory<Grid>& _factory;
static constexpr int dimension
The dimension of the grid.
Definition: grid.hh:387
typename GridFamily::Traits::Communication Communication
A type that is a model of Dune::Communication. It provides a portable way for communication on the se...
Definition: grid.hh:515
ct ctype
Define type used for coordinates in grid module.
Definition: grid.hh:518
A factory class for conveniently creating tensorproduct grids.
Definition: tensorgridfactory.hh:42
void fillRange(int d, int n, ctype end)
fills the range to end with n intervals in direction d
Definition: tensorgridfactory.hh:108
void geometricFillRange(int d, int n, ctype end, ctype h=static_cast< ctype >(0), bool first=true)
fills a coordinate range in direction d with n intervals according to a geometric series
Definition: tensorgridfactory.hh:197
void fillUntil(int d, ctype h, ctype end)
adds intervals in direction d until a given coordinate is reached
Definition: tensorgridfactory.hh:125
void geometricFillUntil(int d, ctype ratio, ctype end, ctype h0=static_cast< ctype >(0))
adds intervals in direction d according with a given length ratio until a given coordinate is reached
Definition: tensorgridfactory.hh:170
void print()
print the coordinate information given to the factory so far
Definition: tensorgridfactory.hh:219
void fillIntervals(int d, int n, ctype h)
pushs n intervals of length h in direction d
Definition: tensorgridfactory.hh:93
void setStart(int d, ctype value)
set a starting value in a given direction d
Definition: tensorgridfactory.hh:79
void geometricFillIntervals(int d, int n, ctype ratio, ctype h0=static_cast< ctype >(0))
adds n intervals in direction d with a given length ratio and a given starting interval length.
Definition: tensorgridfactory.hh:144
std::vector< ctype > & operator[](std::size_t d)
allow to manually tune the factory by overloading operator[] to export the coordinate vectors in the ...
Definition: tensorgridfactory.hh:60
const std::vector< ctype > & operator[](std::size_t d) const
allow to manually tune the factory by overloading operator[] to export the coordinate vectors in the ...
Definition: tensorgridfactory.hh:66
Provide a generic factory class for unstructured grids.
Implements a multiindex with arbitrary dimension and fixed index ranges This is used by various facto...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integer_sequence< T, II..., T(IN)> push_back(std::integer_sequence< T, II... >, std::integral_constant< T, IN >={})
Append an index IN to the back of the sequence.
Definition: integersequence.hh:69
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::bool_constant<(sizeof...(II)==0)> empty(std::integer_sequence< T, II... >)
Checks whether the sequence is empty.
Definition: integersequence.hh:80
constexpr auto back(std::integer_sequence< T, II... > seq)
Return the last entry of the sequence.
Definition: integersequence.hh:44