utilities.h

Go to the documentation of this file.

/*
  Copyright (C) 2014 - 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_utilities_h
#define _aspect_utilities_h
 
#include <aspect/global.h>
 
#include <array>
#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 GeometryModel
  {
    template <int dim> class Interface;
  }
 
  namespace Utilities
  {
    using namespace ::Utilities;
 
 
    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>
    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 std::vector<Point<dim>> &,
                       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>
    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);
    }
 
  }
#endif
}
 
#endif