ASPECT
compat.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2015 - 2020 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef _aspect_compat_h
22 #define _aspect_compat_h
23 
24 #include <aspect/global.h>
25 
26 // C++11 related includes.
27 #include <array>
28 #include <functional>
29 #include <memory>
30 
31 // for std_cxx14::make_unique:
32 #if DEAL_II_VERSION_GTE(9,3,0)
33 // avoid deprecated std_cxx14 inside deal.II
34 namespace std_cxx14
35 {
36  using std::make_unique;
37 }
38 #else
40 #endif
41 
42 #if DEAL_II_VERSION_GTE(9,1,0)
44 namespace aspect
45 {
54  using ConstraintMatrix = class ::AffineConstraints<double>;
55 }
56 #endif
57 
58 #if !DEAL_II_VERSION_GTE(9,2,0)
59 #include <deal.II/base/table.h>
61 namespace aspect
62 {
63  namespace Functions
64  {
65  using namespace ::Functions;
66 
67  // Pull in a couple of symbols from deal.II so that the implementation below
68  // does not need to qualify everything with ::. Note, we do not want
69  // to "using namespace dealii" because this will pull in many symbols, including
70  // deprecated ZeroFunction, into this namespace.
71  using ::Point;
72  using ::Tensor;
73  using ::Table;
74  using ::TableIndices;
76 
83  template <int dim>
85  {
86  public:
102  const std::array<std::pair<double, double>, dim> &interval_endpoints,
103  const std::array<unsigned int, dim> &n_subintervals,
104  const Table<dim, double> &data_values);
105 
117  gradient(const Point<dim> &p,
118  const unsigned int component = 0) const override;
119 
130  double
131  value(const Point<dim> &p,
132  const unsigned int component = 0) const override;
133 
134 
135  private:
139  const std::array<std::pair<double, double>, dim> interval_endpoints;
140 
144  const std::array<unsigned int, dim> n_subintervals;
145 
150  };
151  }
152 }
153 
154 namespace aspect
155 {
156  namespace Functions
157  {
158  namespace
159  {
160  // interpolate a data value from a table where ix denotes
161  // the (lower) left endpoint of the interval to interpolate
162  // in, and p_unit denotes the point in unit coordinates to do so.
163  inline
164  double
165  interpolate(const Table<1, double> &data_values,
166  const TableIndices<1> &ix,
167  const Point<1> &xi)
168  {
169  return ((1 - xi[0]) * data_values[ix[0]] +
170  xi[0] * data_values[ix[0] + 1]);
171  }
172 
173  inline
174  double
175  interpolate(const Table<2, double> &data_values,
176  const TableIndices<2> &ix,
177  const Point<2> &p_unit)
178  {
179  return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
180  p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
181  (1 - p_unit[1]) +
182  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
183  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
184  p_unit[1]);
185  }
186 
187  inline
188  double
189  interpolate(const Table<3, double> &data_values,
190  const TableIndices<3> &ix,
191  const Point<3> &p_unit)
192  {
193  return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
194  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
195  (1 - p_unit[1]) +
196  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
197  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
198  p_unit[1]) *
199  (1 - p_unit[2]) +
200  (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
201  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
202  (1 - p_unit[1]) +
203  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
204  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
205  p_unit[1]) *
206  p_unit[2]);
207  }
208 
209 
210  // Interpolate the gradient of a data value from a table where ix
211  // denotes the lower left endpoint of the interval to interpolate
212  // in, p_unit denotes the point in unit coordinates, and dx
213  // denotes the width of the interval in each dimension.
214  inline
216  gradient_interpolate(const Table<1, double> &data_values,
217  const TableIndices<1> &ix,
218  const Point<1> &p_unit,
219  const Point<1> &dx)
220  {
221  (void)p_unit;
222  Tensor<1, 1> grad;
223  grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
224  return grad;
225  }
226 
227 
228  inline
230  gradient_interpolate(const Table<2, double> &data_values,
231  const TableIndices<2> &ix,
232  const Point<2> &p_unit,
233  const Point<2> &dx)
234  {
235  Tensor<1, 2> grad;
236  double u00 = data_values[ix[0]][ix[1]],
237  u01 = data_values[ix[0] + 1][ix[1]],
238  u10 = data_values[ix[0]][ix[1] + 1],
239  u11 = data_values[ix[0] + 1][ix[1] + 1];
240 
241  grad[0] =
242  ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
243  grad[1] =
244  ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
245  return grad;
246  }
247 
248  inline
250  gradient_interpolate(const Table<3, double> &data_values,
251  const TableIndices<3> &ix,
252  const Point<3> &p_unit,
253  const Point<3> &dx)
254  {
255  Tensor<1, 3> grad;
256  double u000 = data_values[ix[0]][ix[1]][ix[2]],
257  u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
258  u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
259  u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
260  u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
261  u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
262  u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
263  u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
264 
265  grad[0] =
266  ((1 - p_unit[2]) *
267  ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
268  p_unit[2] *
269  ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
270  dx[0];
271  grad[1] =
272  ((1 - p_unit[2]) *
273  ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
274  p_unit[2] *
275  ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
276  dx[1];
277  grad[2] =
278  ((1 - p_unit[1]) *
279  ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
280  p_unit[1] *
281  ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
282  dx[2];
283 
284  return grad;
285  }
286  } // namespace internal
287 
288  template <int dim>
289  inline
291  const std::array<std::pair<double, double>, dim> &interval_endpoints,
292  const std::array<unsigned int, dim> &n_subintervals,
293  const Table<dim, double> &data_values)
294  :
295  interval_endpoints(interval_endpoints),
296  n_subintervals(n_subintervals),
297  data_values(data_values)
298  {
299  for (unsigned int d = 0; d < dim; ++d)
300  {
301  Assert(n_subintervals[d] >= 1,
302  ExcMessage("There needs to be at least one subinterval in each "
303  "coordinate direction."));
305  ExcMessage("The interval in each coordinate direction needs "
306  "to have positive size"));
307  Assert(data_values.size()[d] == n_subintervals[d] + 1,
308  ExcMessage("The data table does not have the correct size."));
309  }
310  }
311 
312 
317  template <int dim>
318  inline
321  const Point<dim> &p,
322  const unsigned int component) const
323  {
324  (void)component;
325  Assert(
326  component == 0,
327  ExcMessage(
328  "This is a scalar function object, the component can only be zero."));
329 
330  // find out where this data point lies, relative to the given
331  // subdivision points
333  for (unsigned int d = 0; d < dim; ++d)
334  {
335  const double delta_x =
336  ((this->interval_endpoints[d].second - this->interval_endpoints[d].first) /
337  this->n_subintervals[d]);
338  if (p[d] <= this->interval_endpoints[d].first)
339  ix[d] = 0;
340  else if (p[d] >= this->interval_endpoints[d].second - delta_x)
341  ix[d] = this->n_subintervals[d] - 1;
342  else
343  ix[d] = static_cast<unsigned int>(
344  (p[d] - this->interval_endpoints[d].first) / delta_x);
345  }
346 
347  // now compute the relative point within the interval/rectangle/box
348  // defined by the point coordinates found above. truncate below and
349  // above to accommodate points that may lie outside the range
350  Point<dim> p_unit;
351  Point<dim> delta_x;
352  for (unsigned int d = 0; d < dim; ++d)
353  {
354  delta_x[d] =
355  ((this->interval_endpoints[d].second - this->interval_endpoints[d].first) /
356  this->n_subintervals[d]);
357  p_unit[d] = std::max(std::min((p[d] - this->interval_endpoints[d].first -
358  ix[d] * delta_x[d]) /
359  delta_x[d],
360  1.),
361  0.);
362  }
363 
364  return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
365  }
366 
367  template <int dim>
368  inline
369  double
371  const unsigned int component) const
372  {
373  (void)component;
374  Assert(
375  component == 0,
376  ExcMessage(
377  "This is a scalar function object, the component can only be zero."));
378 
379  // find out where this data point lies, relative to the given
380  // subdivision points
382  for (unsigned int d = 0; d < dim; ++d)
383  {
384  const double delta_x =
385  ((interval_endpoints[d].second - interval_endpoints[d].first) /
386  n_subintervals[d]);
387  if (p[d] <= interval_endpoints[d].first)
388  ix[d] = 0;
389  else if (p[d] >= interval_endpoints[d].second - delta_x)
390  ix[d] = n_subintervals[d] - 1;
391  else
392  ix[d] = static_cast<unsigned int>(
393  (p[d] - interval_endpoints[d].first) / delta_x);
394  }
395 
396  // now compute the relative point within the interval/rectangle/box
397  // defined by the point coordinates found above. truncate below and
398  // above to accommodate points that may lie outside the range
399  Point<dim> p_unit;
400  for (unsigned int d = 0; d < dim; ++d)
401  {
402  const double delta_x =
403  ((interval_endpoints[d].second - interval_endpoints[d].first) /
404  n_subintervals[d]);
405  p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
406  ix[d] * delta_x) /
407  delta_x,
408  1.),
409  0.);
410  }
411 
412  return interpolate(data_values, ix, p_unit);
413  }
414  }
415 
416 }
417 #endif
418 
419 #endif
const Table< dim, double > data_values
Definition: compat.h:149
dx
Point< 2 > second
double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: compat.h:370
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
static const bool value
const std::array< std::pair< double, double >, dim > interval_endpoints
Definition: compat.h:139
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Point< 2 > first
Definition: compat.h:61
const std::array< unsigned int, dim > n_subintervals
Definition: compat.h:144
size_type size(const unsigned int i) const
Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: compat.h:320