compat.h

Go to the documentation of this file.

/*
  Copyright (C) 2015 - 2024 by the authors of the ASPECT code.
 
  This file is part of ASPECT.
 
  ASPECT is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2, or (at your option)
  any later version.
 
  ASPECT is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with ASPECT; see the file LICENSE.  If not see
  <http://www.gnu.org/licenses/>.
*/
 
#ifndef _aspect_compat_h
#define _aspect_compat_h
 
#include <deal.II/base/config.h>
#include <deal.II/base/mpi.h>
 
// C++11 related includes.
#include <functional>
#include <memory>
 
#if !DEAL_II_VERSION_GTE(9,6,0)
#  include <deal.II/grid/manifold.h>
#  include <deal.II/grid/manifold_lib.h>
#else
// We have a 'using' declaration for new enough deal.II versions. When
// we remove the backward compatibility, we should also remove the
// 'using' declaration and then can also remove the following
// #include:
#  include <deal.II/grid/manifold_lib.h>
#endif
 
#if !DEAL_II_VERSION_GTE(9,7,0)
#  include <deal.II/grid/grid_generator.h>
#endif
 
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_transfer_global_coarsening.h>
 
namespace aspect
{
  namespace big_mpi
  {
 
    using ::Utilities::MPI::broadcast;
 
  }
 
  // deal.II 9.6 introduced MGTransferMF as a replacement for
  // MGTransferMatrixFree; deal.II 9.9 renamed it back.
#if DEAL_II_VERSION_GTE(9,9,0)
  template <int dim, typename NumberType>
  using MGTransferType = ::MGTransferMatrixFree<dim, NumberType>;
#else
  template <int dim, typename NumberType>
  using MGTransferType = ::MGTransferMF<dim, NumberType>;
#endif
 
  // Global-coarsening transfers share the same class as local smoothing
  // from deal.II 9.7 onward; 9.6 still needs MGTransferGlobalCoarsening.
#if DEAL_II_VERSION_GTE(9,7,0)
  template <int dim, typename NumberType>
  using GCMGTransferType = MGTransferType<dim, NumberType>;
#else
  template <int dim, typename NumberType>
  using GCMGTransferType = ::MGTransferGlobalCoarsening<dim, ::LinearAlgebra::distributed::Vector<NumberType>>;
#endif
 
 
#if !DEAL_II_VERSION_GTE(9,7,0)
  template <typename T, typename P = void>
  using ObserverPointer = ::SmartPointer<T, P>;
 
  using EnableObserverPointer = ::Subscriptor;
#endif
 
 
// deal.II versions up to 9.5 had a poorly designed interface of the
// SphericalManifold class that made it impossible for us to use.
// This file thus contains a copy of it.
#if !DEAL_II_VERSION_GTE(9,6,0)
  using namespace dealii;
 
  template <int dim, int spacedim = dim>
  class SphericalManifold : public Manifold<dim, spacedim>
  {
    public:
      SphericalManifold(const Point<spacedim> center = Point<spacedim>());
 
      virtual std::unique_ptr<Manifold<dim, spacedim>>
      clone() const override;
 
      virtual Point<spacedim>
      get_intermediate_point(const Point<spacedim> &p1,
                             const Point<spacedim> &p2,
                             const double           w) const override;
 
      virtual Tensor<1, spacedim>
      get_tangent_vector(const Point<spacedim> &x1,
                         const Point<spacedim> &x2) const override;
 
      virtual Tensor<1, spacedim>
      normal_vector(
        const typename Triangulation<dim, spacedim>::face_iterator &face,
        const Point<spacedim> &p) const override;
 
      virtual void
      get_normals_at_vertices(
        const typename Triangulation<dim, spacedim>::face_iterator &face,
        typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
      const override;
 
      virtual void
      get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                     const Table<2, double>                 &weights,
                     ArrayView<Point<spacedim>> new_points) const override;
 
      virtual Point<spacedim>
      get_new_point(const ArrayView<const Point<spacedim>> &vertices,
                    const ArrayView<const double>          &weights) const override;
 
      const Point<spacedim> center;
 
    private:
      std::pair<double, Tensor<1, spacedim>>
      guess_new_point(const ArrayView<const Tensor<1, spacedim>> &directions,
                      const ArrayView<const double>              &distances,
                      const ArrayView<const double>              &weights) const;
 
      void
      do_get_new_points(const ArrayView<const Point<spacedim>> &surrounding_points,
                        const ArrayView<const double>          &weights,
                        ArrayView<Point<spacedim>>              new_points) const;
 
      const PolarManifold<spacedim> polar_manifold;
  };
 
#else
 
// For sufficiently new deal.II versions, we can use the deal.II class, but to
// avoid name clashes, we have to import the class into namespace aspect. Once
// we rely on these sufficiently new versions of deal.II, we can not only remove
// the code above, but also the following lines, and in all places where we
// reference 'aspect::SphericalManifold' simply use 'SphericalManifold' instead
// (which then refers to the deal.II class).
 
  using ::SphericalManifold;
 
#endif
 
 
// deal.II version 9.7 introduces a new class VectorFunctionFromTensorFunctionObject
// that we would like to use also for earlier versions
#if !DEAL_II_VERSION_GTE(9,7,0)
 
  using namespace dealii;
 
  template <int dim, typename RangeNumberType = double>
  class VectorFunctionFromTensorFunctionObject
    : public Function<dim, RangeNumberType>
  {
    public:
      explicit 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);
 
      virtual ~VectorFunctionFromTensorFunctionObject() override = default;
 
      virtual RangeNumberType
      value(const Point<dim> &p, const unsigned int component = 0) const override;
 
      virtual void
      vector_value(const Point<dim>        &p,
                   Vector<RangeNumberType> &values) const override;
 
      virtual void
      vector_value_list(
        const std::vector<Point<dim>>        &points,
        std::vector<Vector<RangeNumberType>> &value_list) const override;
 
    private:
      const std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>
      tensor_function_object;
 
      const unsigned int selected_component;
  };
 
 
  template <int dim, typename RangeNumberType>
  VectorFunctionFromTensorFunctionObject<dim, RangeNumberType>::
  VectorFunctionFromTensorFunctionObject(
    const std::function<Tensor<1, dim, RangeNumberType>(const Point<dim> &)>
    &tensor_function_object,
    const unsigned int selected_component,
    const unsigned int n_components)
    : Function<dim, RangeNumberType>(n_components)
    , tensor_function_object(tensor_function_object)
    , selected_component(selected_component)
  {
    // Verify that the Tensor<1,dim,RangeNumberType> will fit in the given length
    // selected_components and not hang over the end of the vector.
    AssertIndexRange(selected_component + dim - 1, this->n_components);
  }
 
 
 
  template <int dim, typename RangeNumberType>
  inline RangeNumberType
  VectorFunctionFromTensorFunctionObject<dim, RangeNumberType>::value(
    const Point<dim>  &p,
    const unsigned int component) const
  {
    AssertIndexRange(component, this->n_components);
 
    // if the requested component is out of the range selected, then we can
    // return early
    if ((component < selected_component) ||
        (component >= selected_component + dim))
      return 0;
 
    // otherwise retrieve the values from the <tt>tensor_function</tt> to be
    // placed at the <tt>selected_component</tt> to
    // <tt>selected_component + dim - 1</tt> elements of the <tt>Vector</tt>
    // values and pick the correct one
    const Tensor<1, dim, RangeNumberType> tensor_value =
      tensor_function_object(p);
 
    return tensor_value[component - selected_component];
  }
 
 
  template <int dim, typename RangeNumberType>
  inline void
  VectorFunctionFromTensorFunctionObject<dim, RangeNumberType>::vector_value(
    const Point<dim>        &p,
    Vector<RangeNumberType> &values) const
  {
    Assert(values.size() == this->n_components,
           ExcDimensionMismatch(values.size(), this->n_components));
 
    // Retrieve the values from the <tt>tensor_function</tt> to be placed at
    // the <tt>selected_component</tt> to
    // <tt>selected_component + dim - 1</tt> elements of the <tt>Vector</tt>
    // values.
    const Tensor<1, dim, RangeNumberType> tensor_value =
      tensor_function_object(p);
 
    // First we make all elements of values = 0
    values = 0;
 
    // Second we adjust the desired components to take on the values in
    // <tt>tensor_value</tt>.
    for (unsigned int i = 0; i < dim; ++i)
      values(i + selected_component) = tensor_value[i];
  }
 
 
  template <int dim, typename RangeNumberType>
  void
  VectorFunctionFromTensorFunctionObject<dim, RangeNumberType>::vector_value_list(
    const std::vector<Point<dim>>        &points,
    std::vector<Vector<RangeNumberType>> &value_list) const
  {
    Assert(value_list.size() == points.size(),
           ExcDimensionMismatch(value_list.size(), points.size()));
 
    const unsigned int n_points = points.size();
 
    for (unsigned int p = 0; p < n_points; ++p)
      VectorFunctionFromTensorFunctionObject<dim, RangeNumberType>::vector_value(
        points[p], value_list[p]);
  }
 
#endif
 
// deal.II versions up to 9.6 had a bug for very thin shell geometries.
// This function contains a fixed version.
#if !DEAL_II_VERSION_GTE(9,7,0)
 
  template <int dim>
  void
  colorize_quarter_hyper_shell(Triangulation<dim>  &tria,
                               const Point<dim>   &center,
                               const double      inner_radius,
                               const double      outer_radius);
#endif
 
  // deal.II 9.8 made ReferenceCell a template class, whereas older versions
  // had it as a non-template class. This is a problem.
  // Rather than litter our own code base with #ifdefs, we can just define the
  // templated class variant here for older deal.II versions, and then we can
  // use the same code in all versions.
#if !DEAL_II_VERSION_GTE(9,8,0)
  template <int dim> using ReferenceCell = ::ReferenceCell;
#endif
 
}
 
#endif