utilities.h

Go to the documentation of this file.

// ------------------------------------------------------
//
// SPDX-License-Identifier: GPL-2.0-or-later
//
// SPDX-FileCopyrightText: Copyright (C) 2014-2026 by the ASPECT authors.
//
// This file is part of ASPECT.
//
// Detailed license information governing the source code and
// contributions can be found in the folder LICENSES and in
// LICENSE and CONTRIBUTING.md at the top level directory.
//
// ------------------------------------------------------


#ifndef _aspect_utilities_h
#define _aspect_utilities_h

#include <aspect/global.h>

#include <array>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/thread_local_storage.h>
#include <random>
#include <deal.II/base/point.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/table_indices.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/physics/notation.h>

#include <aspect/coordinate_systems.h>
#include <aspect/structured_data.h>

#include <mpi.h>


namespace aspect
{
  template <int dim> class SimulatorAccess;
  namespace MaterialModel
  {
    template <int dim> class MeltOutputs;
  }

  namespace GeometryModel
  {
    template <int dim> class Interface;
  }

  namespace Utilities
  {
    using namespace ::Utilities;


    template <typename T>
    class ScratchSpace
    {
      public:
        class ScopedScratchObject
        {
          public:
            ScopedScratchObject (const ScratchSpace<T> &/*space_*/);

            ~ScopedScratchObject();

            operator T &() const;

          private:
            const ScratchSpace &space;
            T &t;
        };

        T &get_object_from_pool() const;

        ~ScratchSpace() = default;

        void return_object_to_pool (T &t) const;

      private:
        mutable ::Threads::ThreadLocalStorage<std::list<std::pair<T,bool>>>  object_list;
    };

    template <typename T>
    std::vector<T>
    possibly_extend_from_1_to_N (const std::vector<T> &values,
                                 const unsigned int N,
                                 const std::string &id_text);

    namespace MapParsing
    {
      struct Options
      {
        std::vector<std::string> list_of_allowed_keys;

        std::vector<std::string> list_of_required_keys;

        std::string property_name;

        bool allow_multiple_values_per_key;

        /*
         * Whether to allow for some keys in list_of_required_keys to be
         * not set to any values, i.e. they do not appear at all.
         * This also allows a completely empty map.
         */
        bool allow_missing_keys;

        bool store_values_per_key;

        bool check_values_per_key;

        std::vector<unsigned int> n_values_per_key;

        Options() = delete;

        Options(const std::vector<std::string> &list_of_required_keys,
                const std::string &property_name)
          :
          list_of_allowed_keys(list_of_required_keys),
          list_of_required_keys(list_of_required_keys),
          property_name(property_name),
          allow_multiple_values_per_key(false),
          allow_missing_keys(false),
          store_values_per_key(false),
          check_values_per_key(false),
          n_values_per_key()
        {}
      };

      std::vector<double>
      parse_map_to_double_array(const std::string &input_string,
                                Options &options);
    }

    DEAL_II_DEPRECATED
    std::vector<double>
    parse_map_to_double_array (const std::string &key_value_map,
                               const std::vector<std::string> &list_of_keys,
                               const bool expects_background_field,
                               const std::string &property_name,
                               const bool allow_multiple_values_per_key = false,
                               const std::unique_ptr<std::vector<unsigned int>> &n_values_per_key = nullptr,
                               const bool allow_missing_keys = false);

    template <typename T>
    Table<2,T>
    parse_input_table (const std::string &input_string,
                       const unsigned int n_rows,
                       const unsigned int n_columns,
                       const std::string &property_name);

    template <int dim>
    std::vector<std::string>
    expand_dimensional_variable_names (const std::vector<std::string> &var_declarations);


    template <int dim>
    Tensor<1, dim>
    calculate_approximate_darcy_velocity (const MaterialModel::MaterialModelInputs<dim> &in,
                                          const MaterialModel::MaterialModelOutputs<dim> &out,
                                          const std::shared_ptr<const MaterialModel::MeltOutputs<dim>> fluid_out,
                                          const Tensor<1, dim> &solid_velocity,
                                          const Tensor<1, dim> &gravity,
                                          const unsigned int porosity_idx,
                                          const unsigned int q,
                                          const bool use_pressure_gradient_for_darcy_field);
    template <int dim>
    IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &dof_handler,
                                                        const ComponentMask &component_mask);

    template <int dim>
    std::vector<Point<dim>> get_unit_support_points(const SimulatorAccess<dim> &simulator_access);



    template <int dim>
    bool
    point_is_in_triangulation(const Mapping<dim> &mapping,
                              const parallel::distributed::Triangulation<dim> &triangulation,
                              const Point<dim> &point,
                              const MPI_Comm mpi_communicator);


    namespace Coordinates
    {

      template <int dim>
      std::array<double,dim>
      WGS84_coordinates(const ::Point<dim> &position);

      template <int dim>
      std::array<double,dim>
      cartesian_to_spherical_coordinates(const ::Point<dim> &position);

      template <int dim>
      ::Point<dim>
      spherical_to_cartesian_coordinates(const std::array<double,dim> &scoord);

      template <int dim>
      Tensor<1,dim>
      spherical_to_cartesian_vector(const Tensor<1,dim> &spherical_vector,
                                    const ::Point<dim> &position);


      template <int dim>
      std::array<double,3>
      cartesian_to_ellipsoidal_coordinates(const ::Point<3> &position,
                                           const double semi_major_axis_a,
                                           const double eccentricity);

      template <int dim>
      ::Point<3>
      ellipsoidal_to_cartesian_coordinates(const std::array<double,3> &phi_theta_d,
                                           const double semi_major_axis_a,
                                           const double eccentricity);


      CoordinateSystem
      string_to_coordinate_system (const std::string &);
    }


    template <int dim>
    bool
    polygon_contains_point(const std::vector<Point<2>> &point_list,
                           const ::Point<2> &point);

    template <int dim>
    double
    signed_distance_to_polygon(const std::vector<Point<2>> &point_list,
                               const ::Point<2> &point);


    double
    distance_to_line(const std::array<::Point<2>,2> &point_list,
                     const ::Point<2> &point);

    template <int dim>
    std::array<Tensor<1,dim>,dim-1>
    orthogonal_vectors (const Tensor<1,dim> &v);

    Tensor<2,3>
    rotation_matrix_from_axis (const Tensor<1,3> &rotation_axis,
                               const double rotation_angle);

    Tensor<2,3>
    compute_rotation_matrix_for_slice (const Tensor<1,3> &point_one,
                                       const Tensor<1,3> &point_two);

    std::pair<double,double> real_spherical_harmonic( unsigned int l, // degree
                                                      unsigned int m, // order
                                                      double theta,   // colatitude (radians)
                                                      double phi );   // longitude (radians)

    struct ThousandSep : std::numpunct<char>
    {
      protected:
        char do_thousands_sep() const override
        {
          return ',';
        }

        std::string do_grouping() const override
        {
          return "\003";  // groups of 3 digits (this string is in octal format)
        }

    };

    bool fexists(const std::string &filename);

    bool fexists(const std::string &filename,
                 const MPI_Comm comm);

    bool filename_is_url(const std::string &filename);

    std::string
    read_and_distribute_file_content(const std::string &filename,
                                     const MPI_Comm comm);

    void
    collect_and_write_file_content(const std::string &filename,
                                   const std::string &file_content,
                                   const MPI_Comm comm);

    int
    mkdirp(std::string pathname, const mode_t mode = 0755);

    void create_directory(const std::string &pathname,
                          const MPI_Comm comm,
                          const bool silent);

    namespace tk
    {
      class spline
      {
        public:
          void set_points(const std::vector<double> &x,
                          const std::vector<double> &y,
                          const bool cubic_spline = true,
                          const bool monotone_spline = false);
          double operator() (double x) const;

        private:
          std::vector<double> m_x;

          std::vector<double> m_a, m_b, m_c, m_y;
      };
    }

    inline
    void
    extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
                                           const unsigned int q,
                                           std::vector<double> &composition_values_at_q_point);

    std::string
    expand_ASPECT_SOURCE_DIR (const std::string &location);

    std::string parenthesize_if_nonempty (const std::string &s);

    bool
    string_to_bool(const std::string &s);

    std::vector<bool>
    string_to_bool(const std::vector<std::string> &s);

    unsigned int
    string_to_unsigned_int(const std::string &s);

    std::vector<unsigned int>
    string_to_unsigned_int(const std::vector<std::string> &s);

    bool has_unique_entries (const std::vector<std::string> &strings);



    double weighted_p_norm_average (const std::vector<double> &weights,
                                    const std::vector<double> &values,
                                    const double p);


    template <typename T>
    T derivative_of_weighted_p_norm_average (const double averaged_parameter,
                                             const std::vector<double> &weights,
                                             const std::vector<double> &values,
                                             const std::vector<T> &derivatives,
                                             const double p);
    template <int dim>
    double compute_spd_factor(const double eta,
                              const SymmetricTensor<2,dim> &strain_rate,
                              const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
                              const double SPD_safety_factor);

    template <int dim>
    Point<dim> convert_array_to_point(const std::array<double,dim> &array);

    template <int dim>
    std::array<double,dim> convert_point_to_array(const Point<dim> &point);

    class Operator
    {
      public:
        enum operation
        {
          uninitialized,
          add,
          subtract,
          minimum,
          maximum,
          replace_if_valid
        };

        Operator();

        Operator(const operation op);

        double operator() (const double x, const double y) const;

        bool operator== (const operation op) const;

      private:
        operation op;
    };

    std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);

    const std::string get_model_operator_options();

    template <int dim>
    SymmetricTensor<2,dim> nth_basis_for_symmetric_tensors (const unsigned int k);

    template <int dim>
    class NaturalCoordinate
    {
      public:
        NaturalCoordinate(Point<dim> &position,
                          const GeometryModel::Interface<dim> &geometry_model);

        NaturalCoordinate(const std::array<double, dim> &coord,
                          const Utilities::Coordinates::CoordinateSystem &coord_system);

        std::array<double,dim> &get_coordinates();

        const std::array<double,dim> &get_coordinates() const;

        std::array<double,dim-1> get_surface_coordinates() const;

        double get_depth_coordinate() const;

      private:
        Utilities::Coordinates::CoordinateSystem coordinate_system;

        std::array<double,dim> coordinates;
    };


    template <int dim, typename VectorType>
    void
    project_cellwise(const Mapping<dim>                                        &mapping,
                     const DoFHandler<dim>                                     &dof_handler,
                     const unsigned int                                         component_index,
                     const Quadrature<dim>                                     &quadrature,
                     const std::function<void(
                       const typename DoFHandler<dim>::active_cell_iterator &,
                       const typename std_cxx20::type_identity<std::vector<Point<dim>>>::type &,
                       std::vector<double> &)>                                 &function,
                     VectorType                                                &vec_result);

    void throw_linear_solver_failure_exception(const std::string &solver_name,
                                               const std::string &function_name,
                                               const std::vector<SolverControl> &solver_controls,
                                               const std::exception &exc,
                                               const MPI_Comm mpi_communicator,
                                               const std::string &output_filename = "");

    template <int dim>
    class VectorFunctionFromVelocityFunctionObject : public VectorFunctionFromTensorFunctionObject<dim>
    {
      public:
        VectorFunctionFromVelocityFunctionObject (const unsigned int n_components,
                                                  const std::function<Tensor<1,dim> (const Point<dim> &)> &function_object);
    };

    template <typename T>
    inline
    std::vector<std::size_t>
    compute_sorting_permutation(const std::vector<T> &vector);

    template <typename T>
    inline
    std::vector<T>
    apply_permutation(
      const std::vector<T> &vector,
      const std::vector<std::size_t> &permutation_vector);

    std::vector<Tensor<2,3>>
    rotation_matrices_random_draw_volume_weighting(const std::vector<double> &volume_fractions,
                                                   const std::vector<Tensor<2,3>> &rotation_matrices,
                                                   const unsigned int n_output_matrices,
                                                   std::mt19937 &random_number_generator);

    double wrap_angle(const double angle);

    std::array<double,3> zxz_euler_angles_from_rotation_matrix(const Tensor<2,3> &rotation_matrix);

    Tensor<2,3> zxz_euler_angles_to_rotation_matrix(const double phi1,
                                                    const double theta,
                                                    const double phi2);

  }


// inline implementations:
#ifndef DOXYGEN
  namespace Utilities
  {
    template<typename T>
    T &ScratchSpace<T>::get_object_from_pool() const
    {
      for (auto &pair : object_list.get())
        if (pair.second == false)
          {
            pair.second = true;
            return pair.first;
          }

      object_list.get().emplace_back (T(), true);
      return object_list.get().back().first;
    }

    template<typename T>
    void ScratchSpace<T>::return_object_to_pool (T &t) const
    {
      for (auto &pair : object_list.get())
        if (&pair.first == &t)
          {
            pair.second = false;
            return;
          }
      AssertThrow(false, ExcMessage("You are tying to return an object to the pool which has apparently not been allocated by this pool."));
    }

    template<typename T>
    ScratchSpace<T>::ScopedScratchObject::ScopedScratchObject(const ScratchSpace<T> &space_)
      : space (space_),
        t (space.get_object_from_pool())
    {}

    template<typename T>
    ScratchSpace<T>::ScopedScratchObject::~ScopedScratchObject()
    {
      space.return_object_to_pool(t);
    }

    template<typename T>
    ScratchSpace<T>::ScopedScratchObject::operator T &() const
    {
      return t;
    }

    template <typename T>
    inline
    std::vector<T>
    possibly_extend_from_1_to_N (const std::vector<T> &values,
                                 const unsigned int N,
                                 const std::string &id_text)
    {
      if (values.size() == 1)
        {
          return std::vector<T> (N, values[0]);
        }
      else if (values.size() == N)
        {
          return values;
        }
      else
        {
          // Non-specified behavior
          AssertThrow(false,
                      ExcMessage("Length of " + id_text + " list must be " +
                                 "either one or " + Utilities::to_string(N) +
                                 ". Currently it is " + Utilities::to_string(values.size()) + "."));
        }

      // This should never happen, but return an empty vector so the compiler
      // will be happy
      return std::vector<T> ();
    }

    inline
    void
    extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
                                           const unsigned int q,
                                           std::vector<double> &composition_values_at_q_point)
    {
      Assert(q<composition_values.size(), ExcInternalError());
      Assert(composition_values_at_q_point.size() > 0,
             ExcInternalError());

      for (unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
        {
          Assert(composition_values[k].size() == composition_values_at_q_point.size(),
                 ExcInternalError());
          composition_values_at_q_point[k] = composition_values[k][q];
        }
    }

    template <typename T>
    inline
    std::vector<std::size_t>
    compute_sorting_permutation(const std::vector<T> &vector)
    {
      std::vector<std::size_t> p(vector.size());
      std::iota(p.begin(), p.end(), 0);
      std::sort(p.begin(), p.end(),
                [&](std::size_t i, std::size_t j)
      {
        return vector[i] < vector[j];
      });
      return p;
    }

    template <typename T>
    inline
    std::vector<T>
    apply_permutation(
      const std::vector<T> &vector,
      const std::vector<std::size_t> &permutation_vector)
    {
      std::vector<T> sorted_vec(vector.size());
      std::transform(permutation_vector.begin(), permutation_vector.end(), sorted_vec.begin(),
                     [&](std::size_t i)
      {
        return vector[i];
      });
      return sorted_vec;
    }

    namespace Tensors
    {
      template <int dim, class Iterator>
      inline
      SymmetricTensor<2,dim>
      to_symmetric_tensor(const Iterator begin,
                          const Iterator end)
      {
        AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
        (void) end;

        SymmetricTensor<2,dim> output;

        Iterator next = begin;
        for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
          output[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)] = *next;

        return output;
      }

      template <int dim, class Iterator>
      inline
      void
      unroll_symmetric_tensor_into_array(const SymmetricTensor<2,dim> &tensor,
                                         const Iterator begin,
                                         const Iterator end)
      {
        AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
        (void) end;

        Iterator next = begin;
        for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
          *next = tensor[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)];
      }

      SymmetricTensor<4,3>
      rotate_full_stiffness_tensor(const Tensor<2,3> &rotation_tensor, const SymmetricTensor<4,3> &input_tensor);

      SymmetricTensor<2,6>
      rotate_voigt_stiffness_matrix(const Tensor<2,3> &rotation_tensor, const SymmetricTensor<2,6> &input_tensor);

      SymmetricTensor<2,6>
      rotate_kelvin_tensor(const Tensor<2,3> &rotation_tensor, const SymmetricTensor<2,6> &input_tensor);

      SymmetricTensor<2,6>
      to_voigt_stiffness_matrix(const SymmetricTensor<4,3> &input_tensor);

      SymmetricTensor<4,3>
      to_full_stiffness_tensor(const SymmetricTensor<2,6> &input_tensor);

      Tensor<1,21>
      to_voigt_stiffness_vector(const SymmetricTensor<2,6> &input_tensor);

      SymmetricTensor<2,6>
      to_voigt_stiffness_matrix(const Tensor<1,21> &input_tensor);

      Tensor<1,21>
      to_voigt_stiffness_vector(const SymmetricTensor<4,3> &input);


      template <int dim>
      const Tensor<dim,dim> &levi_civita();

      // Declare the existence of a specialization:
      template <>
      const Tensor<3,3> &levi_civita<3>();

      template <int dim>
      SymmetricTensor<2,dim>
      consistent_deviator(const SymmetricTensor<2,dim> &input);

      template <int dim>
      double
      consistent_second_invariant_of_deviatoric_tensor(const SymmetricTensor<2,dim> &input);
    }

    namespace Quaternions
    {
      std::array<double,4> rotation_matrix_to_quaternion(const Tensor<2,3> &rotation_matrix,
                                                         const double tolerance=1e-12);

      Tensor<2,3> quaternion_to_rotation_matrix(const std::array<double,4> &quaternion);
    }

  }
#endif
}

#endif