ASPECT
compat.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2015 - 2019 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 #include <deal.II/base/std_cxx14/memory.h>
33 
34 #if DEAL_II_VERSION_GTE(9,1,0)
35 #include <deal.II/lac/affine_constraints.h>
36 namespace aspect
37 {
46  using ConstraintMatrix = class ::AffineConstraints<double>;
47 }
48 #endif
49 
50 #if !DEAL_II_VERSION_GTE(9,2,0)
51 #include <deal.II/base/table.h>
52 #include <deal.II/base/function_lib.h>
53 namespace aspect
54 {
55  namespace Functions
56  {
57  using namespace dealii;
58  using namespace ::Functions;
65  template <int dim>
67  {
68  public:
84  const std::array<std::pair<double, double>, dim> &interval_endpoints,
85  const std::array<unsigned int, dim> &n_subintervals,
86  const Table<dim, double> &data_values);
87 
99  gradient(const Point<dim> &p,
100  const unsigned int component = 0) const override;
101 
112  double
113  value(const Point<dim> &p,
114  const unsigned int component = 0) const override;
115 
116 
117  private:
121  const std::array<std::pair<double, double>, dim> interval_endpoints;
122 
126  const std::array<unsigned int, dim> n_subintervals;
127 
132  };
133  }
134 }
135 
136 namespace aspect
137 {
138  namespace Functions
139  {
140  namespace
141  {
142  // interpolate a data value from a table where ix denotes
143  // the (lower) left endpoint of the interval to interpolate
144  // in, and p_unit denotes the point in unit coordinates to do so.
145  inline
146  double
147  interpolate(const Table<1, double> &data_values,
148  const TableIndices<1> &ix,
149  const Point<1> &xi)
150  {
151  return ((1 - xi[0]) * data_values[ix[0]] +
152  xi[0] * data_values[ix[0] + 1]);
153  }
154 
155  inline
156  double
157  interpolate(const Table<2, double> &data_values,
158  const TableIndices<2> &ix,
159  const Point<2> &p_unit)
160  {
161  return (((1 - p_unit[0]) * data_values[ix[0]][ix[1]] +
162  p_unit[0] * data_values[ix[0] + 1][ix[1]]) *
163  (1 - p_unit[1]) +
164  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1] +
165  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1]) *
166  p_unit[1]);
167  }
168 
169  inline
170  double
171  interpolate(const Table<3, double> &data_values,
172  const TableIndices<3> &ix,
173  const Point<3> &p_unit)
174  {
175  return ((((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2]] +
176  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2]]) *
177  (1 - p_unit[1]) +
178  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2]] +
179  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2]]) *
180  p_unit[1]) *
181  (1 - p_unit[2]) +
182  (((1 - p_unit[0]) * data_values[ix[0]][ix[1]][ix[2] + 1] +
183  p_unit[0] * data_values[ix[0] + 1][ix[1]][ix[2] + 1]) *
184  (1 - p_unit[1]) +
185  ((1 - p_unit[0]) * data_values[ix[0]][ix[1] + 1][ix[2] + 1] +
186  p_unit[0] * data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1]) *
187  p_unit[1]) *
188  p_unit[2]);
189  }
190 
191 
192  // Interpolate the gradient of a data value from a table where ix
193  // denotes the lower left endpoint of the interval to interpolate
194  // in, p_unit denotes the point in unit coordinates, and dx
195  // denotes the width of the interval in each dimension.
196  inline
198  gradient_interpolate(const Table<1, double> &data_values,
199  const TableIndices<1> &ix,
200  const Point<1> &p_unit,
201  const Point<1> &dx)
202  {
203  (void)p_unit;
204  Tensor<1, 1> grad;
205  grad[0] = (data_values[ix[0] + 1] - data_values[ix[0]]) / dx[0];
206  return grad;
207  }
208 
209 
210  inline
212  gradient_interpolate(const Table<2, double> &data_values,
213  const TableIndices<2> &ix,
214  const Point<2> &p_unit,
215  const Point<2> &dx)
216  {
217  Tensor<1, 2> grad;
218  double u00 = data_values[ix[0]][ix[1]],
219  u01 = data_values[ix[0] + 1][ix[1]],
220  u10 = data_values[ix[0]][ix[1] + 1],
221  u11 = data_values[ix[0] + 1][ix[1] + 1];
222 
223  grad[0] =
224  ((1 - p_unit[1]) * (u01 - u00) + p_unit[1] * (u11 - u10)) / dx[0];
225  grad[1] =
226  ((1 - p_unit[0]) * (u10 - u00) + p_unit[0] * (u11 - u01)) / dx[1];
227  return grad;
228  }
229 
230  inline
232  gradient_interpolate(const Table<3, double> &data_values,
233  const TableIndices<3> &ix,
234  const Point<3> &p_unit,
235  const Point<3> &dx)
236  {
237  Tensor<1, 3> grad;
238  double u000 = data_values[ix[0]][ix[1]][ix[2]],
239  u001 = data_values[ix[0] + 1][ix[1]][ix[2]],
240  u010 = data_values[ix[0]][ix[1] + 1][ix[2]],
241  u100 = data_values[ix[0]][ix[1]][ix[2] + 1],
242  u011 = data_values[ix[0] + 1][ix[1] + 1][ix[2]],
243  u101 = data_values[ix[0] + 1][ix[1]][ix[2] + 1],
244  u110 = data_values[ix[0]][ix[1] + 1][ix[2] + 1],
245  u111 = data_values[ix[0] + 1][ix[1] + 1][ix[2] + 1];
246 
247  grad[0] =
248  ((1 - p_unit[2]) *
249  ((1 - p_unit[1]) * (u001 - u000) + p_unit[1] * (u011 - u010)) +
250  p_unit[2] *
251  ((1 - p_unit[1]) * (u101 - u100) + p_unit[1] * (u111 - u110))) /
252  dx[0];
253  grad[1] =
254  ((1 - p_unit[2]) *
255  ((1 - p_unit[0]) * (u010 - u000) + p_unit[0] * (u011 - u001)) +
256  p_unit[2] *
257  ((1 - p_unit[0]) * (u110 - u100) + p_unit[0] * (u111 - u101))) /
258  dx[1];
259  grad[2] =
260  ((1 - p_unit[1]) *
261  ((1 - p_unit[0]) * (u100 - u000) + p_unit[0] * (u101 - u001)) +
262  p_unit[1] *
263  ((1 - p_unit[0]) * (u110 - u010) + p_unit[0] * (u111 - u011))) /
264  dx[2];
265 
266  return grad;
267  }
268  } // namespace internal
269 
270  template <int dim>
271  inline
273  const std::array<std::pair<double, double>, dim> &interval_endpoints,
274  const std::array<unsigned int, dim> &n_subintervals,
275  const Table<dim, double> &data_values)
276  :
277  interval_endpoints(interval_endpoints),
278  n_subintervals(n_subintervals),
279  data_values(data_values)
280  {
281  for (unsigned int d = 0; d < dim; ++d)
282  {
283  Assert(n_subintervals[d] >= 1,
284  ExcMessage("There needs to be at least one subinterval in each "
285  "coordinate direction."));
286  Assert(interval_endpoints[d].first < interval_endpoints[d].second,
287  ExcMessage("The interval in each coordinate direction needs "
288  "to have positive size"));
289  Assert(data_values.size()[d] == n_subintervals[d] + 1,
290  ExcMessage("The data table does not have the correct size."));
291  }
292  }
293 
294 
299  template <int dim>
300  inline
303  const Point<dim> &p,
304  const unsigned int component) const
305  {
306  (void)component;
307  Assert(
308  component == 0,
309  ExcMessage(
310  "This is a scalar function object, the component can only be zero."));
311 
312  // find out where this data point lies, relative to the given
313  // subdivision points
315  for (unsigned int d = 0; d < dim; ++d)
316  {
317  const double delta_x =
318  ((this->interval_endpoints[d].second - this->interval_endpoints[d].first) /
319  this->n_subintervals[d]);
320  if (p[d] <= this->interval_endpoints[d].first)
321  ix[d] = 0;
322  else if (p[d] >= this->interval_endpoints[d].second - delta_x)
323  ix[d] = this->n_subintervals[d] - 1;
324  else
325  ix[d] = static_cast<unsigned int>(
326  (p[d] - this->interval_endpoints[d].first) / delta_x);
327  }
328 
329  // now compute the relative point within the interval/rectangle/box
330  // defined by the point coordinates found above. truncate below and
331  // above to accommodate points that may lie outside the range
332  Point<dim> p_unit;
333  Point<dim> delta_x;
334  for (unsigned int d = 0; d < dim; ++d)
335  {
336  delta_x[d] =
337  ((this->interval_endpoints[d].second - this->interval_endpoints[d].first) /
338  this->n_subintervals[d]);
339  p_unit[d] = std::max(std::min((p[d] - this->interval_endpoints[d].first -
340  ix[d] * delta_x[d]) /
341  delta_x[d],
342  1.),
343  0.);
344  }
345 
346  return gradient_interpolate(this->data_values, ix, p_unit, delta_x);
347  }
348 
349  template <int dim>
350  inline
351  double
353  const unsigned int component) const
354  {
355  (void)component;
356  Assert(
357  component == 0,
358  ExcMessage(
359  "This is a scalar function object, the component can only be zero."));
360 
361  // find out where this data point lies, relative to the given
362  // subdivision points
364  for (unsigned int d = 0; d < dim; ++d)
365  {
366  const double delta_x =
367  ((interval_endpoints[d].second - interval_endpoints[d].first) /
368  n_subintervals[d]);
369  if (p[d] <= interval_endpoints[d].first)
370  ix[d] = 0;
371  else if (p[d] >= interval_endpoints[d].second - delta_x)
372  ix[d] = n_subintervals[d] - 1;
373  else
374  ix[d] = static_cast<unsigned int>(
375  (p[d] - interval_endpoints[d].first) / delta_x);
376  }
377 
378  // now compute the relative point within the interval/rectangle/box
379  // defined by the point coordinates found above. truncate below and
380  // above to accommodate points that may lie outside the range
381  Point<dim> p_unit;
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  p_unit[d] = std::max(std::min((p[d] - interval_endpoints[d].first -
388  ix[d] * delta_x) /
389  delta_x,
390  1.),
391  0.);
392  }
393 
394  return interpolate(data_values, ix, p_unit);
395  }
396  }
397 
398 }
399 #endif
400 
401 #endif
const Table< dim, double > data_values
Definition: compat.h:131
double value(const Point< dim > &p, const unsigned int component=0) const override
Definition: compat.h:352
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:121
Definition: compat.h:53
const std::array< unsigned int, dim > n_subintervals
Definition: compat.h:126
InterpolatedUniformGridData(const std::array< std::pair< double, double >, dim > &interval_endpoints, const std::array< unsigned int, dim > &n_subintervals, const Table< dim, double > &data_values)
Definition: compat.h:272
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:302