Dune Core Modules (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>
305  class TensorGridFactoryCreator
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  {
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>
388  class TensorGridFactoryCreator<YaspGrid<dim, TensorProductCoordinates<ctype, dim> > >
389  {
390  public:
391  typedef YaspGrid<dim, TensorProductCoordinates<ctype, dim> > Grid;
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
constexpr static int dimension
The dimension of the grid.
Definition: grid.hh:387
detected_or_fallback_t< DeprecatedCollectiveCommunication_t, Communication_t, typename GridFamily::Traits > Communication
A type that is a model of Dune::Communication. It provides a portable way for communication on the se...
Definition: grid.hh:525
ct ctype
Define type used for coordinates in grid module.
Definition: grid.hh:532
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
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
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
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:472
constexpr auto back(const HybridTreePath< T... > &tp) -> decltype(treePathEntry< sizeof...(T) -1 >(tp))
Returns a copy of the last element of the HybridTreePath.
Definition: treepath.hh:257
constexpr HybridTreePath< T..., std::size_t > push_back(const HybridTreePath< T... > &tp, std::size_t i)
Appends a run time index to a HybridTreePath.
Definition: treepath.hh:281
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Mar 27, 23:31, 2024)