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