ASPECT
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
compat.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2015 - 2024 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 <deal.II/base/config.h>
25 #include <deal.II/base/mpi.h>
26 
27 // C++11 related includes.
28 #include <array>
29 #include <functional>
30 #include <memory>
31 
32 #if !DEAL_II_VERSION_GTE(9,6,0)
33 # include <deal.II/multigrid/mg_transfer_matrix_free.h>
34 # include <deal.II/grid/manifold.h>
35 # include <deal.II/grid/manifold_lib.h>
36 #else
37 // We have a 'using' declaration for new enough deal.II versions. When
38 // we remove the backward compatibility, we should also remove the
39 // 'using' declaration and then can also remove the following
40 // #include:
41 # include <deal.II/grid/manifold_lib.h>
42 #endif
43 
44 #if !DEAL_II_VERSION_GTE(9,7,0)
45 # include <deal.II/grid/grid_generator.h>
46 #endif
47 
48 
49 namespace aspect
50 {
51  namespace big_mpi
52  {
53 
54  using ::Utilities::MPI::broadcast;
55 
56  }
57 
58 
59 // deal.II 9.6 introduces the new MGTransferMF class as a replacement
60 // for MGTransferMatrixFree. Instead of putting an ifdef in every place,
61 // do this in one central location:
62 #if !DEAL_II_VERSION_GTE(9,6,0)
63  template <int dim, class NumberType>
64  using MGTransferMF = ::MGTransferMatrixFree<dim,NumberType>;
65 #endif
66 
67 
68 
69 // deal.II versions up to 9.5 had a poorly designed interface of the
70 // SphericalManifold class that made it impossible for us to use.
71 // This file thus contains a copy of it.
72 #if !DEAL_II_VERSION_GTE(9,6,0)
73  using namespace dealii;
74 
83  template <int dim, int spacedim = dim>
84  class SphericalManifold : public Manifold<dim, spacedim>
85  {
86  public:
93  SphericalManifold(const Point<spacedim> center = Point<spacedim>());
94 
98  virtual std::unique_ptr<Manifold<dim, spacedim>>
99  clone() const override;
100 
108  virtual Point<spacedim>
109  get_intermediate_point(const Point<spacedim> &p1,
110  const Point<spacedim> &p2,
111  const double w) const override;
112 
117  virtual Tensor<1, spacedim>
118  get_tangent_vector(const Point<spacedim> &x1,
119  const Point<spacedim> &x2) const override;
120 
124  virtual Tensor<1, spacedim>
125  normal_vector(
126  const typename Triangulation<dim, spacedim>::face_iterator &face,
127  const Point<spacedim> &p) const override;
128 
132  virtual void
133  get_normals_at_vertices(
134  const typename Triangulation<dim, spacedim>::face_iterator &face,
135  typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
136  const override;
137 
152  virtual void
153  get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
154  const Table<2, double> &weights,
155  ArrayView<Point<spacedim>> new_points) const override;
156 
161  virtual Point<spacedim>
162  get_new_point(const ArrayView<const Point<spacedim>> &vertices,
163  const ArrayView<const double> &weights) const override;
164 
168  const Point<spacedim> center;
169 
170  private:
177  std::pair<double, Tensor<1, spacedim>>
178  guess_new_point(const ArrayView<const Tensor<1, spacedim>> &directions,
179  const ArrayView<const double> &distances,
180  const ArrayView<const double> &weights) const;
181 
199  void
200  do_get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
201  const ArrayView<const double> &weights,
202  ArrayView<Point<spacedim>> new_points) const;
203 
207  const PolarManifold<spacedim> polar_manifold;
208  };
209 
210 #else
211 
212 // For sufficiently new deal.II versions, we can use the deal.II class, but to
213 // avoid name clashes, we have to import the class into namespace aspect. Once
214 // we rely on these sufficiently new versions of deal.II, we can not only remove
215 // the code above, but also the following lines, and in all places where we
216 // reference 'aspect::SphericalManifold' simply use 'SphericalManifold' instead
217 // (which then refers to the deal.II class).
218 
219  using ::SphericalManifold;
220 
221 #endif
222 
223 
224 // deal.II version 9.7 introduces a new class VectorFunctionFromTensorFunctionObject
225 // that we would like to use also for earlier versions
226 #if !DEAL_II_VERSION_GTE(9,7,0)
227 
228  using namespace dealii;
229 
264  template <int dim, typename RangeNumberType = double>
266  : public Function<dim, RangeNumberType>
267  {
268  public:
286  const std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>
287  &tensor_function_object,
288  const unsigned int selected_component = 0,
289  const unsigned int n_components = dim);
290 
295  virtual ~VectorFunctionFromTensorFunctionObject() override = default;
296 
300  virtual RangeNumberType
301  value(const Point<dim> &p, const unsigned int component = 0) const override;
302 
308  virtual void
309  vector_value(const Point<dim> &p,
310  Vector<RangeNumberType> &values) const override;
311 
319  virtual void
320  vector_value_list(
321  const std::vector<Point<dim>> &points,
322  std::vector<Vector<RangeNumberType>> &value_list) const override;
323 
324  private:
329  const std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>
330  tensor_function_object;
331 
338  const unsigned int selected_component;
339  };
340 
341 
342  template <int dim, typename RangeNumberType>
345  const std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>
346  &tensor_function_object,
347  const unsigned int selected_component,
348  const unsigned int n_components)
349  : Function<dim, RangeNumberType>(n_components)
350  , tensor_function_object(tensor_function_object)
351  , selected_component(selected_component)
352  {
353  // Verify that the Tensor<1,dim,RangeNumberType> will fit in the given length
354  // selected_components and not hang over the end of the vector.
355  AssertIndexRange(selected_component + dim - 1, this->n_components);
356  }
357 
358 
359 
360  template <int dim, typename RangeNumberType>
361  inline RangeNumberType
363  const Point<dim> &p,
364  const unsigned int component) const
365  {
366  AssertIndexRange(component, this->n_components);
367 
368  // if the requested component is out of the range selected, then we can
369  // return early
370  if ((component < selected_component) ||
371  (component >= selected_component + dim))
372  return 0;
373 
374  // otherwise retrieve the values from the <tt>tensor_function</tt> to be
375  // placed at the <tt>selected_component</tt> to
376  // <tt>selected_component + dim - 1</tt> elements of the <tt>Vector</tt>
377  // values and pick the correct one
378  const Tensor<1, dim, RangeNumberType> tensor_value =
379  tensor_function_object(p);
380 
381  return tensor_value[component - selected_component];
382  }
383 
384 
385  template <int dim, typename RangeNumberType>
386  inline void
388  const Point<dim> &p,
389  Vector<RangeNumberType> &values) const
390  {
391  Assert(values.size() == this->n_components,
392  ExcDimensionMismatch(values.size(), this->n_components));
393 
394  // Retrieve the values from the <tt>tensor_function</tt> to be placed at
395  // the <tt>selected_component</tt> to
396  // <tt>selected_component + dim - 1</tt> elements of the <tt>Vector</tt>
397  // values.
398  const Tensor<1, dim, RangeNumberType> tensor_value =
399  tensor_function_object(p);
400 
401  // First we make all elements of values = 0
402  values = 0;
403 
404  // Second we adjust the desired components to take on the values in
405  // <tt>tensor_value</tt>.
406  for (unsigned int i = 0; i < dim; ++i)
407  values(i + selected_component) = tensor_value[i];
408  }
409 
410 
418  template <int dim, typename RangeNumberType>
419  void
421  const std::vector<Point<dim>> &points,
422  std::vector<Vector<RangeNumberType>> &value_list) const
423  {
424  Assert(value_list.size() == points.size(),
425  ExcDimensionMismatch(value_list.size(), points.size()));
426 
427  const unsigned int n_points = points.size();
428 
429  for (unsigned int p = 0; p < n_points; ++p)
431  points[p], value_list[p]);
432  }
433 
434 #endif
435 
436 // deal.II versions up to 9.6 had a bug for very thin shell geometries.
437 // This function contains a fixed version.
438 #if !DEAL_II_VERSION_GTE(9,7,0)
439 
440  template <int dim>
441  void
442  colorize_quarter_hyper_shell(Triangulation<dim> &tria,
443  const Point<dim> &center,
444  const double inner_radius,
445  const double outer_radius);
446 #endif
447 
448 }
449 
450 #endif
VectorFunctionFromTensorFunctionObject(const std::function< Tensor< 1, dim, RangeNumberType >(const Point< dim > &)> &tensor_function_object, const unsigned int selected_component=0, const unsigned int n_components=dim)
Definition: compat.h:344
::MGTransferMatrixFree< dim, NumberType > MGTransferMF
Definition: compat.h:64
const PolarManifold< spacedim > polar_manifold
Definition: compat.h:207
void colorize_quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > &center, const double inner_radius, const double outer_radius)
const Point< spacedim > center
Definition: compat.h:168