dune-grid  2.9.0
tensorgridfactory.hh
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3 #ifndef DUNE_GRID_UTILITY_TENSORGRIDFACTORY_HH
4 #define DUNE_GRID_UTILITY_TENSORGRIDFACTORY_HH
5 
20 #include<array>
21 #include<memory>
22 #include<vector>
23 
24 #include <dune/common/fvector.hh>
26 #include <dune/grid/yaspgrid.hh>
28 
29 namespace Dune
30 {
31  // forward declaration of TensorGridFactoryCreator, which is the real factory
32  // that should be specialized for each grid.
33  template<typename Grid>
34  class TensorGridFactoryCreator;
35 
40  template<typename Grid>
42  {
43  public:
44  typedef typename Grid::Traits::Communication Comm;
45  typedef typename Grid::ctype ctype;
46  static const int dim = Grid::dimension;
47 
48  std::unique_ptr<Grid> createGrid(Comm comm = Comm())
49  {
50  TensorGridFactoryCreator<Grid> creator(*this);
51  return creator.createGrid(comm);
52  }
53 
54  std::array<std::vector<ctype> , dim> coords() const
55  {
56  return _coords;
57  }
58 
60  std::vector<ctype>& operator[](std::size_t d)
61  {
62  return _coords[d];
63  }
64 
66  const std::vector<ctype>& operator[](std::size_t d) const
67  {
68  return _coords[d];
69  }
70 
79  void setStart (int d, ctype value)
80  {
81  _coords[d].resize(1);
82  _coords[d][0] = value;
83  }
84 
93  void fillIntervals (int d, int n, ctype h)
94  {
95  emptyCheck (d);
96  for (int i = 0; i < n; i++)
97  _coords[d].push_back (_coords[d].back () + h);
98  }
99 
108  void fillRange (int d, int n, ctype end)
109  {
110  emptyCheck (d);
111  const ctype h = (end - _coords[d].back ()) / n;
112  for (int i = 0; i < n - 1; i++)
113  _coords[d].push_back (_coords[d].back () + h);
114  _coords[d].push_back (end);
115  }
116 
125  void fillUntil (int d, ctype h, ctype end)
126  {
127  emptyCheck (d);
128  while (_coords[d].back () < end)
129  _coords[d].push_back (_coords[d].back () + h);
130  }
131 
144  void geometricFillIntervals (int d, int n, ctype ratio, ctype h0 =
145  static_cast<ctype> (0))
146  {
147  emptyCheck (d);
148  ctype h = h0;
149  if (h0 == static_cast<ctype>(0))
150  h = lastInterval (d) * ratio;
151  for (int i = 0; i < n; i++)
152  {
153  _coords[d].push_back (_coords[d].back () + h);
154  h *= ratio;
155  }
156  }
157 
170  void geometricFillUntil (int d, ctype ratio, ctype end, ctype h0 = static_cast<ctype> (0))
171  {
172  emptyCheck (d);
173  ctype h = h0;
174  if (h0 == static_cast<ctype>(0))
175  h = lastInterval (d) * ratio;
176  while (_coords[d].back () < end)
177  {
178  _coords[d].push_back (_coords[d].back () + h);
179  h *= ratio;
180  }
181  }
182 
197  void geometricFillRange (int d, int n, ctype end, ctype h =
198  static_cast<ctype> (0),
199  bool first = true)
200  {
201  emptyCheck (d);
202  if (h < 1e-8)
203  h = lastInterval (d);
204  ctype ratio = newton (n, _coords[d].back (), end, h);
205  if (!first)
206  {
207  h = h * pow (ratio, n - 1);
208  ratio = 1 / ratio;
209  }
210  for (int i = 0; i < n - 1; i++)
211  {
212  _coords[d].push_back (_coords[d].back () + h);
213  h *= ratio;
214  }
215  _coords[d].push_back (end);
216  }
217 
219  void print()
220  {
221  for (int i=0; i<dim; i++)
222  {
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: ";
227 
228  std::vector<ctype> meshsize;
229  for (auto it = _coords[i].begin(); it != _coords[i].end()-1;)
230  {
231  meshsize.push_back(-1.*(*it));
232  ++it;
233  meshsize.back() += *it;
234  }
235 
236  for (auto it = meshsize.begin(); it != meshsize.end(); ++it)
237  std::cout << *it << " ";
238  std::cout << std::endl << "Ratios between interval lengths: ";
239 
240  std::vector<ctype> ratios;
241  for (auto it = meshsize.begin(); it != meshsize.end() - 1 ;)
242  ratios.push_back((1./(*it)) * *(++it));
243 
244  for (auto it = ratios.begin(); it != ratios.end(); ++it)
245  std::cout << *it << " ";
246  std::cout << std::endl << std::endl << std::endl;
247  }
248  }
249 
250  private:
251  // check whether the ith component is empty and add a 0.0 entry if so
252  void emptyCheck (int i)
253  {
254  if (_coords[i].empty ())
255  _coords[i].push_back (static_cast<ctype> (0));
256  }
257 
258  // returns the last interval length in direction d
259  ctype lastInterval (int d)
260  {
261  if (_coords[d].size () < 2)
262  DUNE_THROW(
263  GridError,
264  "Not enough elements in coordinate container to deduce interval length in TensorYaspFactory");
265  else
266  return _coords[d].back () - _coords[d][_coords[d].size () - 2];
267  }
268 
272  ctype newton (int n, ctype x_s, ctype x_e, ctype h)
273  {
274  ctype m = (x_e - x_s) / h;
275  ctype xold = 0.0;
276  ctype xnew = x_e - x_s;
277  while (std::abs (xnew - xold) > 1E-8)
278  {
279  xold = xnew;
280  xnew = xold
281  - (-pow (xold, n) + m * xold - m + 1)
282  / (-n * pow (xold, n - 1) + m);
283  }
284  if (std::abs (xnew - 1) < 1E-6)
285  {
286  xold = x_e - x_s;
287  xnew = 0.0;
288  while (std::abs (xnew - xold) > 1E-8)
289  {
290  xold = xnew;
291  xnew = xold
292  - (-pow (xold, n) + m * xold - m + 1)
293  / (-n * pow (xold, n - 1) + m);
294  }
295  }
296  return xnew;
297  }
298 
299  std::array<std::vector<ctype>, dim> _coords;
300  };
301 
302  // class that implements the actual grid creation process. The default is implementing
303  // standard creation for unstructured grids. Provide a specialization for other grids.
304  template<typename Grid>
306  {
307  public:
308  typedef typename Grid::Traits::Communication Comm;
309  typedef typename Grid::ctype ctype;
310  static const int dim = Grid::dimension;
311 
312  TensorGridFactoryCreator(const TensorGridFactory<Grid>& factory) : _factory(factory) {}
313 
314  std::unique_ptr<Grid> createGrid(Comm comm)
315  {
316  // The grid factory
317  GridFactory<Grid> fac;
318 
319  if (comm.rank() == 0)
320  {
321  // determine the size of the grid
322  std::array<unsigned int, dim> vsizes, esizes;
323  std::size_t size = 1;
324  for (std::size_t i = 0; i<dim; ++i)
325  {
326  vsizes[i] = _factory[i].size();
327  esizes[i] = vsizes[i] - 1;
328  size *= vsizes[i];
329  }
330 
331  // insert all vertices
332  FactoryUtilities::MultiIndex<dim> index(vsizes);
333  for (std::size_t i=0; i<size; ++i, ++index)
334  {
335  Dune::FieldVector<ctype, dim> position;
336  for (std::size_t j = 0; j<dim; ++j)
337  position[j] = _factory[j][index[j]];
338  fac.insertVertex(position);
339  }
340 
341  // compute the offsets
342  std::array<std::size_t, dim> offsets;
343  offsets[0] = 1;
344  for (std::size_t i=1; i<dim; i++)
345  offsets[i] = offsets[i-1] * vsizes[i-1];
346 
347  // Compute an element template (the cube at (0,...,0). All
348  // other cubes are constructed by moving this template around
349  unsigned int nCorners = 1<<dim;
350 
351  std::vector<unsigned int> cornersTemplate(nCorners,0);
352 
353  for (size_t i=0; i<nCorners; i++)
354  for (int j=0; j<dim; j++)
355  if ( i & (1<<j) )
356  cornersTemplate[i] += offsets[j];
357 
358  // Insert elements
359  FactoryUtilities::MultiIndex<dim> eindex(esizes);
360 
361  // Compute the total number of elements to be created
362  int numElements = eindex.cycle();
363 
364  for (int i=0; i<numElements; i++, ++eindex)
365  {
366  // 'base' is the index of the lower left element corner
367  unsigned int base = 0;
368  for (int j=0; j<dim; j++)
369  base += eindex[j] * offsets[j];
370 
371  // insert new element
372  std::vector<unsigned int> corners = cornersTemplate;
373  for (size_t j=0; j<corners.size(); j++)
374  corners[j] += base;
375 
376  fac.insertElement(GeometryTypes::cube(dim), corners);
377  }
378  }
379 
380  return std::unique_ptr<Grid>(fac.createGrid());
381  }
382 
383  private:
384  const TensorGridFactory<Grid>& _factory;
385  };
386 
387  template<typename ctype, int dim>
389  {
390  public:
392  typedef typename Grid::Communication Comm;
393 
394  TensorGridFactoryCreator(const TensorGridFactory<Grid>& factory) : _factory(factory) {}
395 
396  std::unique_ptr<Grid> createGrid(Comm comm)
397  {
398  return std::make_unique<Grid>(_factory.coords(), std::bitset<dim>(0ULL), 1, comm);
399  }
400  private:
401  const TensorGridFactory<Grid>& _factory;
402  };
403 }
404 
405 #endif
Implements a multiindex with arbitrary dimension and fixed index ranges This is used by various facto...
Include standard header files.
Definition: agrid.hh:60
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:328
constexpr static int dimension
The dimension of the grid.
Definition: common/grid.hh:387
ct ctype
Define type used for coordinates in grid module.
Definition: common/grid.hh:532
virtual void insertVertex([[maybe_unused]] const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: common/gridfactory.hh:335
virtual void insertElement([[maybe_unused]] const GeometryType &type, [[maybe_unused]] const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: common/gridfactory.hh:346
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: common/gridfactory.hh:372
Definition: multiindex.hh:19
size_t cycle() const
Compute how many times you can call operator++ before getting to (0,...,0) again.
Definition: multiindex.hh:48
Definition: tensorgridfactory.hh:306
Grid::ctype ctype
Definition: tensorgridfactory.hh:309
Grid::Traits::Communication Comm
Definition: tensorgridfactory.hh:308
TensorGridFactoryCreator(const TensorGridFactory< Grid > &factory)
Definition: tensorgridfactory.hh:312
static const int dim
Definition: tensorgridfactory.hh:310
std::unique_ptr< Grid > createGrid(Comm comm)
Definition: tensorgridfactory.hh:314
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
Grid::ctype ctype
Definition: tensorgridfactory.hh:45
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
static const int dim
Definition: tensorgridfactory.hh:46
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
std::unique_ptr< Grid > createGrid(Comm comm=Comm())
Definition: tensorgridfactory.hh:48
void print()
print the coordinate information given to the factory so far
Definition: tensorgridfactory.hh:219
std::array< std::vector< ctype >, dim > coords() const
Definition: tensorgridfactory.hh:54
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
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
Grid::Traits::Communication Comm
Definition: tensorgridfactory.hh:44
std::unique_ptr< Grid > createGrid(Comm comm)
Definition: tensorgridfactory.hh:396
TensorGridFactoryCreator(const TensorGridFactory< Grid > &factory)
Definition: tensorgridfactory.hh:394
YaspGrid< dim, TensorProductCoordinates< ctype, dim > > Grid
Definition: tensorgridfactory.hh:391
[ provides Dune::Grid ]
Definition: yaspgrid.hh:163
typename Base::Communication Communication
Definition: yaspgrid.hh:178
Coordinate container for a tensor product YaspGrid.
Definition: coordinates.hh:245
Provide a generic factory class for unstructured grids.