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>

#include <deal.II/grid/manifold_lib.h>

#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


  using ::SphericalManifold;



// 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