DUNE-FEM (unstable)

tensorgridfactory.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © 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
26#include <dune/grid/yaspgrid.hh>
28
29namespace 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)
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
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)