ASPECT
utilities.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2014 - 2024 by the authors of the ASPECT code.
3 
4  This file is part of ASPECT.
5 
6  ASPECT is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2, or (at your option)
9  any later version.
10 
11  ASPECT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with ASPECT; see the file LICENSE. If not see
18  <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #ifndef _aspect_utilities_h
23 #define _aspect_utilities_h
24 
25 #include <aspect/global.h>
26 
27 #include <array>
28 #include <random>
29 #include <deal.II/base/point.h>
30 #include <deal.II/base/conditional_ostream.h>
31 #include <deal.II/base/table_indices.h>
32 #include <deal.II/base/function_lib.h>
33 #include <deal.II/dofs/dof_handler.h>
34 #include <deal.II/fe/component_mask.h>
35 
37 #include <aspect/structured_data.h>
38 
39 #include <mpi.h>
40 
41 
42 namespace aspect
43 {
44  template <int dim> class SimulatorAccess;
45 
46  namespace GeometryModel
47  {
48  template <int dim> class Interface;
49  }
50 
55  namespace Utilities
56  {
71  using namespace ::Utilities;
72 
73 
87  template <typename T>
88  std::vector<T>
89  possibly_extend_from_1_to_N (const std::vector<T> &values,
90  const unsigned int N,
91  const std::string &id_text);
92 
100  namespace MapParsing
101  {
106  struct Options
107  {
117  std::vector<std::string> list_of_allowed_keys;
118 
126  std::vector<std::string> list_of_required_keys;
127 
133  std::string property_name;
134 
143 
144  /*
145  * Whether to allow for some keys in list_of_required_keys to be
146  * not set to any values, i.e. they do not appear at all.
147  * This also allows a completely empty map.
148  */
150 
160 
170 
179  std::vector<unsigned int> n_values_per_key;
180 
185  Options() = delete;
186 
197  Options(const std::vector<std::string> &list_of_required_keys,
198  const std::string &property_name)
199  :
200  list_of_allowed_keys(list_of_required_keys),
201  list_of_required_keys(list_of_required_keys),
202  property_name(property_name),
203  allow_multiple_values_per_key(false),
204  allow_missing_keys(false),
205  store_values_per_key(false),
206  check_values_per_key(false),
207  n_values_per_key()
208  {}
209  };
210 
236  std::vector<double>
237  parse_map_to_double_array(const std::string &input_string,
238  Options &options);
239  }
240 
314  DEAL_II_DEPRECATED
315  std::vector<double>
316  parse_map_to_double_array (const std::string &key_value_map,
317  const std::vector<std::string> &list_of_keys,
318  const bool expects_background_field,
319  const std::string &property_name,
320  const bool allow_multiple_values_per_key = false,
321  const std::unique_ptr<std::vector<unsigned int>> &n_values_per_key = nullptr,
322  const bool allow_missing_keys = false);
323 
337  template <typename T>
338  Table<2,T>
339  parse_input_table (const std::string &input_string,
340  const unsigned int n_rows,
341  const unsigned int n_columns,
342  const std::string &property_name);
343 
356  template <int dim>
357  std::vector<std::string>
358  expand_dimensional_variable_names (const std::vector<std::string> &var_declarations);
359 
366  template <int dim>
367  IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &dof_handler,
368  const ComponentMask &component_mask);
369 
379  template <int dim>
380  std::vector<Point<dim>> get_unit_support_points(const SimulatorAccess<dim> &simulator_access);
381 
382 
383 
389  template <int dim>
390  bool
391  point_is_in_triangulation(const Mapping<dim> &mapping,
392  const parallel::distributed::Triangulation<dim> &triangulation,
393  const Point<dim> &point,
394  const MPI_Comm mpi_communicator);
395 
396 
397  namespace Coordinates
398  {
399 
405  template <int dim>
406  std::array<double,dim>
407  WGS84_coordinates(const ::Point<dim> &position);
408 
419  template <int dim>
420  std::array<double,dim>
421  cartesian_to_spherical_coordinates(const ::Point<dim> &position);
422 
428  template <int dim>
429  ::Point<dim>
430  spherical_to_cartesian_coordinates(const std::array<double,dim> &scoord);
431 
437  template <int dim>
438  Tensor<1,dim>
439  spherical_to_cartesian_vector(const Tensor<1,dim> &spherical_vector,
440  const ::Point<dim> &position);
441 
442 
448  template <int dim>
449  std::array<double,3>
450  cartesian_to_ellipsoidal_coordinates(const ::Point<3> &position,
451  const double semi_major_axis_a,
452  const double eccentricity);
453 
458  template <int dim>
459  ::Point<3>
460  ellipsoidal_to_cartesian_coordinates(const std::array<double,3> &phi_theta_d,
461  const double semi_major_axis_a,
462  const double eccentricity);
463 
464 
471  string_to_coordinate_system (const std::string &);
472  }
473 
474 
479  template <int dim>
480  bool
481  polygon_contains_point(const std::vector<Point<2>> &point_list,
482  const ::Point<2> &point);
483 
489  template <int dim>
490  double
491  signed_distance_to_polygon(const std::vector<Point<2>> &point_list,
492  const ::Point<2> &point);
493 
494 
501  double
502  distance_to_line(const std::array<::Point<2>,2> &point_list,
503  const ::Point<2> &point);
504 
512  template <int dim>
513  std::array<Tensor<1,dim>,dim-1>
514  orthogonal_vectors (const Tensor<1,dim> &v);
515 
520  Tensor<2,3>
521  rotation_matrix_from_axis (const Tensor<1,3> &rotation_axis,
522  const double rotation_angle);
523 
530  Tensor<2,3>
531  compute_rotation_matrix_for_slice (const Tensor<1,3> &point_one,
532  const Tensor<1,3> &point_two);
533 
570  std::pair<double,double> real_spherical_harmonic( unsigned int l, // degree
571  unsigned int m, // order
572  double theta, // colatitude (radians)
573  double phi ); // longitude (radians)
574 
578  struct ThousandSep : std::numpunct<char>
579  {
580  protected:
581  char do_thousands_sep() const override
582  {
583  return ',';
584  }
585 
586  std::string do_grouping() const override
587  {
588  return "\003"; // groups of 3 digits (this string is in octal format)
589  }
590 
591  };
592 
605  bool fexists(const std::string &filename);
606 
620  bool fexists(const std::string &filename,
621  const MPI_Comm comm);
622 
628  bool filename_is_url(const std::string &filename);
629 
649  std::string
650  read_and_distribute_file_content(const std::string &filename,
651  const MPI_Comm comm);
652 
665  void
666  collect_and_write_file_content(const std::string &filename,
667  const std::string &file_content,
668  const MPI_Comm comm);
669 
683  int
684  mkdirp(std::string pathname, const mode_t mode = 0755);
685 
697  void create_directory(const std::string &pathname,
698  const MPI_Comm comm,
699  const bool silent);
700 
705  namespace tk
706  {
710  class spline
711  {
712  public:
722  void set_points(const std::vector<double> &x,
723  const std::vector<double> &y,
724  const bool cubic_spline = true,
725  const bool monotone_spline = false);
729  double operator() (double x) const;
730 
731  private:
735  std::vector<double> m_x;
736 
743  std::vector<double> m_a, m_b, m_c, m_y;
744  };
745  }
746 
754  inline
755  void
756  extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
757  const unsigned int q,
758  std::vector<double> &composition_values_at_q_point);
759 
764  std::string
765  expand_ASPECT_SOURCE_DIR (const std::string &location);
766 
771  std::string parenthesize_if_nonempty (const std::string &s);
772 
776  bool
777  string_to_bool(const std::string &s);
778 
782  std::vector<bool>
783  string_to_bool(const std::vector<std::string> &s);
784 
788  unsigned int
789  string_to_unsigned_int(const std::string &s);
790 
794  std::vector<unsigned int>
795  string_to_unsigned_int(const std::vector<std::string> &s);
796 
801  bool has_unique_entries (const std::vector<std::string> &strings);
802 
803 
804 
826  double weighted_p_norm_average (const std::vector<double> &weights,
827  const std::vector<double> &values,
828  const double p);
829 
830 
856  template <typename T>
857  T derivative_of_weighted_p_norm_average (const double averaged_parameter,
858  const std::vector<double> &weights,
859  const std::vector<double> &values,
860  const std::vector<T> &derivatives,
861  const double p);
893  template <int dim>
894  double compute_spd_factor(const double eta,
895  const SymmetricTensor<2,dim> &strain_rate,
896  const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
897  const double SPD_safety_factor);
898 
902  template <int dim>
903  Point<dim> convert_array_to_point(const std::array<double,dim> &array);
904 
908  template <int dim>
909  std::array<double,dim> convert_point_to_array(const Point<dim> &point);
910 
918  class Operator
919  {
920  public:
925  {
931  replace_if_valid
932  };
933 
938  Operator();
939 
943  Operator(const operation op);
944 
949  double operator() (const double x, const double y) const;
950 
955  bool operator== (const operation op) const;
956 
957  private:
962  };
963 
968  std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
969 
973  const std::string get_model_operator_options();
974 
980  template <int dim>
981  SymmetricTensor<2,dim> nth_basis_for_symmetric_tensors (const unsigned int k);
982 
986  template <int dim>
988  {
989  public:
993  NaturalCoordinate(Point<dim> &position,
994  const GeometryModel::Interface<dim> &geometry_model);
995 
1000  NaturalCoordinate(const std::array<double, dim> &coord,
1001  const Utilities::Coordinates::CoordinateSystem &coord_system);
1002 
1007  std::array<double,dim> &get_coordinates();
1008 
1013  const std::array<double,dim> &get_coordinates() const;
1014 
1019  std::array<double,dim-1> get_surface_coordinates() const;
1020 
1025  double get_depth_coordinate() const;
1026 
1027  private:
1033 
1037  std::array<double,dim> coordinates;
1038  };
1039 
1040 
1059  template <int dim, typename VectorType>
1060  void
1061  project_cellwise(const Mapping<dim> &mapping,
1062  const DoFHandler<dim> &dof_handler,
1063  const unsigned int component_index,
1064  const Quadrature<dim> &quadrature,
1065  const std::function<void(
1066  const typename DoFHandler<dim>::active_cell_iterator &,
1067  const std::vector<Point<dim>> &,
1068  std::vector<double> &)> &function,
1069  VectorType &vec_result);
1070 
1094  void throw_linear_solver_failure_exception(const std::string &solver_name,
1095  const std::string &function_name,
1096  const std::vector<SolverControl> &solver_controls,
1097  const std::exception &exc,
1098  const MPI_Comm mpi_communicator,
1099  const std::string &output_filename = "");
1100 
1111  template <int dim>
1113  {
1114  public:
1124  VectorFunctionFromVelocityFunctionObject (const unsigned int n_components,
1125  const std::function<Tensor<1,dim> (const Point<dim> &)> &function_object);
1126  };
1127 
1134  template <typename T>
1135  inline
1136  std::vector<std::size_t>
1137  compute_sorting_permutation(const std::vector<T> &vector);
1138 
1146  template <typename T>
1147  inline
1148  std::vector<T>
1150  const std::vector<T> &vector,
1151  const std::vector<std::size_t> &permutation_vector);
1152 
1165  std::vector<Tensor<2,3>>
1166  rotation_matrices_random_draw_volume_weighting(const std::vector<double> &volume_fractions,
1167  const std::vector<Tensor<2,3>> &rotation_matrices,
1168  const unsigned int n_output_matrices,
1169  std::mt19937 &random_number_generator);
1170 
1174  double wrap_angle(const double angle);
1175 
1180  std::array<double,3> zxz_euler_angles_from_rotation_matrix(const Tensor<2,3> &rotation_matrix);
1181 
1186  Tensor<2,3> zxz_euler_angles_to_rotation_matrix(const double phi1,
1187  const double theta,
1188  const double phi2);
1189 
1190  }
1191 
1192 
1193 // inline implementations:
1194 #ifndef DOXYGEN
1195  namespace Utilities
1196  {
1197 
1198  template <typename T>
1199  inline
1200  std::vector<T>
1201  possibly_extend_from_1_to_N (const std::vector<T> &values,
1202  const unsigned int N,
1203  const std::string &id_text)
1204  {
1205  if (values.size() == 1)
1206  {
1207  return std::vector<T> (N, values[0]);
1208  }
1209  else if (values.size() == N)
1210  {
1211  return values;
1212  }
1213  else
1214  {
1215  // Non-specified behavior
1216  AssertThrow(false,
1217  ExcMessage("Length of " + id_text + " list must be " +
1218  "either one or " + Utilities::to_string(N) +
1219  ". Currently it is " + Utilities::to_string(values.size()) + "."));
1220  }
1221 
1222  // This should never happen, but return an empty vector so the compiler
1223  // will be happy
1224  return std::vector<T> ();
1225  }
1226 
1227  inline
1228  void
1229  extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
1230  const unsigned int q,
1231  std::vector<double> &composition_values_at_q_point)
1232  {
1233  Assert(q<composition_values.size(), ExcInternalError());
1234  Assert(composition_values_at_q_point.size() > 0,
1235  ExcInternalError());
1236 
1237  for (unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
1238  {
1239  Assert(composition_values[k].size() == composition_values_at_q_point.size(),
1240  ExcInternalError());
1241  composition_values_at_q_point[k] = composition_values[k][q];
1242  }
1243  }
1244 
1245  template <typename T>
1246  inline
1247  std::vector<std::size_t>
1248  compute_sorting_permutation(const std::vector<T> &vector)
1249  {
1250  std::vector<std::size_t> p(vector.size());
1251  std::iota(p.begin(), p.end(), 0);
1252  std::sort(p.begin(), p.end(),
1253  [&](std::size_t i, std::size_t j)
1254  {
1255  return vector[i] < vector[j];
1256  });
1257  return p;
1258  }
1259 
1260  template <typename T>
1261  inline
1262  std::vector<T>
1264  const std::vector<T> &vector,
1265  const std::vector<std::size_t> &permutation_vector)
1266  {
1267  std::vector<T> sorted_vec(vector.size());
1268  std::transform(permutation_vector.begin(), permutation_vector.end(), sorted_vec.begin(),
1269  [&](std::size_t i)
1270  {
1271  return vector[i];
1272  });
1273  return sorted_vec;
1274  }
1275 
1279  namespace Tensors
1280  {
1286  template <int dim, class Iterator>
1287  inline
1288  SymmetricTensor<2,dim>
1289  to_symmetric_tensor(const Iterator begin,
1290  const Iterator end)
1291  {
1292  AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
1293  (void) end;
1294 
1295  SymmetricTensor<2,dim> output;
1296 
1297  Iterator next = begin;
1298  for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
1299  output[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)] = *next;
1300 
1301  return output;
1302  }
1303 
1308  template <int dim, class Iterator>
1309  inline
1310  void
1311  unroll_symmetric_tensor_into_array(const SymmetricTensor<2,dim> &tensor,
1312  const Iterator begin,
1313  const Iterator end)
1314  {
1315  AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
1316  (void) end;
1317 
1318  Iterator next = begin;
1319  for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
1320  *next = tensor[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)];
1321  }
1322 
1326  SymmetricTensor<4,3>
1327  rotate_full_stiffness_tensor(const Tensor<2,3> &rotation_tensor, const SymmetricTensor<4,3> &input_tensor);
1328 
1333  SymmetricTensor<2,6>
1334  rotate_voigt_stiffness_matrix(const Tensor<2,3> &rotation_tensor, const SymmetricTensor<2,6> &input_tensor);
1335 
1340  SymmetricTensor<2,6>
1341  to_voigt_stiffness_matrix(const SymmetricTensor<4,3> &input_tensor);
1342 
1347  SymmetricTensor<4,3>
1348  to_full_stiffness_tensor(const SymmetricTensor<2,6> &input_tensor);
1349 
1354  Tensor<1,21>
1355  to_voigt_stiffness_vector(const SymmetricTensor<2,6> &input_tensor);
1356 
1361  SymmetricTensor<2,6>
1362  to_voigt_stiffness_matrix(const Tensor<1,21> &input_tensor);
1363 
1368  Tensor<1,21>
1369  to_voigt_stiffness_vector(const SymmetricTensor<4,3> &input);
1370 
1371 
1376  template <int dim>
1377  const Tensor<dim,dim> &levi_civita();
1378 
1379  // Declare the existence of a specialization:
1380  template <>
1381  const Tensor<3,3> &levi_civita<3>();
1382  }
1383 
1384  }
1385 #endif
1386 }
1387 
1388 #endif
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)
std::string read_and_distribute_file_content(const std::string &filename, const MPI_Comm comm)
Point< dim > convert_array_to_point(const std::array< double, dim > &array)
std::vector< double > m_y
Definition: utilities.h:743
Tensor< 1, dim > spherical_to_cartesian_vector(const Tensor< 1, dim > &spherical_vector, const ::Point< dim > &position)
Tensor< 2, 3 > compute_rotation_matrix_for_slice(const Tensor< 1, 3 > &point_one, const Tensor< 1, 3 > &point_two)
std::vector< std::string > list_of_required_keys
Definition: utilities.h:126
double weighted_p_norm_average(const std::vector< double > &weights, const std::vector< double > &values, const double p)
double distance_to_line(const std::array<::Point< 2 >, 2 > &point_list, const ::Point< 2 > &point)
std::vector< std::string > expand_dimensional_variable_names(const std::vector< std::string > &var_declarations)
std::array< double, dim > coordinates
Definition: utilities.h:1037
void collect_and_write_file_content(const std::string &filename, const std::string &file_content, const MPI_Comm comm)
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)
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)
std::pair< double, double > real_spherical_harmonic(unsigned int l, unsigned int m, double theta, double phi)
std::vector< double > m_x
Definition: utilities.h:735
std::string expand_ASPECT_SOURCE_DIR(const std::string &location)
Tensor< 2, 3 > zxz_euler_angles_to_rotation_matrix(const double phi1, const double theta, const double phi2)
std::array< Tensor< 1, dim >, dim-1 > orthogonal_vectors(const Tensor< 1, dim > &v)
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)
std::string to_string(const Newton::Parameters::Stabilization preconditioner_stabilization)
std::array< double, dim > WGS84_coordinates(const ::Point< dim > &position)
::Point< 3 > ellipsoidal_to_cartesian_coordinates(const std::array< double, 3 > &phi_theta_d, const double semi_major_axis_a, const double eccentricity)
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)
bool point_is_in_triangulation(const Mapping< dim > &mapping, const parallel::distributed::Triangulation< dim > &triangulation, const Point< dim > &point, const MPI_Comm mpi_communicator)
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="")
bool polygon_contains_point(const std::vector< Point< 2 >> &point_list, const ::Point< 2 > &point)
int mkdirp(std::string pathname, const mode_t mode=0755)
Utilities::Coordinates::CoordinateSystem coordinate_system
Definition: utilities.h:1032
void create_directory(const std::string &pathname, const MPI_Comm comm, const bool silent)
::Point< dim > spherical_to_cartesian_coordinates(const std::array< double, dim > &scoord)
std::string do_grouping() const override
Definition: utilities.h:586
std::vector< T > possibly_extend_from_1_to_N(const std::vector< T > &values, const unsigned int N, const std::string &id_text)
std::array< double, 3 > cartesian_to_ellipsoidal_coordinates(const ::Point< 3 > &position, const double semi_major_axis_a, const double eccentricity)
Options(const std::vector< std::string > &list_of_required_keys, const std::string &property_name)
Definition: utilities.h:197
std::array< double, dim > convert_point_to_array(const Point< dim > &point)
std::vector< bool > string_to_bool(const std::vector< std::string > &s)
std::vector< Point< dim > > get_unit_support_points(const SimulatorAccess< dim > &simulator_access)
std::vector< T > apply_permutation(const std::vector< T > &vector, const std::vector< std::size_t > &permutation_vector)
IndexSet extract_locally_active_dofs_with_component(const DoFHandler< dim > &dof_handler, const ComponentMask &component_mask)
std::vector< unsigned int > n_values_per_key
Definition: utilities.h:179
std::vector< unsigned int > string_to_unsigned_int(const std::vector< std::string > &s)
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)
double wrap_angle(const double angle)
bool fexists(const std::string &filename, const MPI_Comm comm)
std::vector< Operator > create_model_operator_list(const std::vector< std::string > &operator_names)
double signed_distance_to_polygon(const std::vector< Point< 2 >> &point_list, const ::Point< 2 > &point)
std::string parenthesize_if_nonempty(const std::string &s)
std::array< double, 3 > zxz_euler_angles_from_rotation_matrix(const Tensor< 2, 3 > &rotation_matrix)
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)
const std::string get_model_operator_options()
char do_thousands_sep() const override
Definition: utilities.h:581
std::vector< std::string > list_of_allowed_keys
Definition: utilities.h:117
bool filename_is_url(const std::string &filename)
SymmetricTensor< 2, dim > nth_basis_for_symmetric_tensors(const unsigned int k)
std::vector< std::size_t > compute_sorting_permutation(const std::vector< T > &vector)
CoordinateSystem string_to_coordinate_system(const std::string &)
std::array< double, dim > cartesian_to_spherical_coordinates(const ::Point< dim > &position)
bool has_unique_entries(const std::vector< std::string > &strings)
Tensor< 2, 3 > rotation_matrix_from_axis(const Tensor< 1, 3 > &rotation_axis, const double rotation_angle)