ASPECT
utilities.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2014 - 2022 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 <deal.II/base/point.h>
29 #include <deal.II/base/conditional_ostream.h>
30 #include <deal.II/base/table_indices.h>
31 #include <deal.II/base/function_lib.h>
32 #include <deal.II/dofs/dof_handler.h>
33 #include <deal.II/fe/component_mask.h>
34 
36 #include <aspect/structured_data.h>
37 
38 
39 
40 namespace aspect
41 {
42  template <int dim> class SimulatorAccess;
43 
44  namespace GeometryModel
45  {
46  template <int dim> class Interface;
47  }
48 
53  namespace Utilities
54  {
55  using namespace dealii;
56  using namespace ::Utilities;
57 
71  template <typename T>
72  std::vector<T>
73  possibly_extend_from_1_to_N (const std::vector<T> &values,
74  const unsigned int N,
75  const std::string &id_text);
76 
146  std::vector<double>
147  parse_map_to_double_array (const std::string &key_value_map,
148  const std::vector<std::string> &list_of_keys,
149  const bool expects_background_field,
150  const std::string &property_name,
151  const bool allow_multiple_values_per_key = false,
152  const std::unique_ptr<std::vector<unsigned int>> &n_values_per_key = nullptr,
153  const bool allow_missing_keys = false);
154 
168  template <typename T>
169  Table<2,T>
170  parse_input_table (const std::string &input_string,
171  const unsigned int n_rows,
172  const unsigned int n_columns,
173  const std::string &property_name);
174 
187  template <int dim>
188  std::vector<std::string>
189  expand_dimensional_variable_names (const std::vector<std::string> &var_declarations);
190 
197  template <int dim>
198  IndexSet extract_locally_active_dofs_with_component(const DoFHandler<dim> &dof_handler,
199  const ComponentMask &component_mask);
200 
210  template <int dim>
211  std::vector<Point<dim>> get_unit_support_points(const SimulatorAccess<dim> &simulator_access);
212 
213 
214 
215  namespace Coordinates
216  {
217 
223  template <int dim>
224  std::array<double,dim>
225  WGS84_coordinates(const ::Point<dim> &position);
226 
233  template <int dim>
234  std::array<double,dim>
235  cartesian_to_spherical_coordinates(const ::Point<dim> &position);
236 
242  template <int dim>
243  ::Point<dim>
244  spherical_to_cartesian_coordinates(const std::array<double,dim> &scoord);
245 
251  template <int dim>
252  Tensor<1,dim>
253  spherical_to_cartesian_vector(const Tensor<1,dim> &spherical_vector,
254  const ::Point<dim> &position);
255 
256 
262  template <int dim>
263  std::array<double,3>
264  cartesian_to_ellipsoidal_coordinates(const ::Point<3> &position,
265  const double semi_major_axis_a,
266  const double eccentricity);
267 
272  template <int dim>
273  ::Point<3>
274  ellipsoidal_to_cartesian_coordinates(const std::array<double,3> &phi_theta_d,
275  const double semi_major_axis_a,
276  const double eccentricity);
277 
278 
285  string_to_coordinate_system (const std::string &);
286  }
287 
288 
293  template <int dim>
294  bool
295  polygon_contains_point(const std::vector<Point<2>> &point_list,
296  const ::Point<2> &point);
297 
303  template <int dim>
304  double
305  signed_distance_to_polygon(const std::vector<Point<2>> &point_list,
306  const ::Point<2> &point);
307 
308 
315  double
316  distance_to_line(const std::array<::Point<2>,2> &point_list,
317  const ::Point<2> &point);
318 
326  template <int dim>
327  std::array<Tensor<1,dim>,dim-1>
328  orthogonal_vectors (const Tensor<1,dim> &v);
329 
334  Tensor<2,3>
335  rotation_matrix_from_axis (const Tensor<1,3> &rotation_axis,
336  const double rotation_angle);
337 
344  Tensor<2,3>
345  compute_rotation_matrix_for_slice (const Tensor<1,3> &point_one,
346  const Tensor<1,3> &point_two);
347 
384  std::pair<double,double> real_spherical_harmonic( unsigned int l, // degree
385  unsigned int m, // order
386  double theta, // colatitude (radians)
387  double phi ); // longitude (radians)
388 
392  struct ThousandSep : std::numpunct<char>
393  {
394  protected:
395  char do_thousands_sep() const override
396  {
397  return ',';
398  }
399 
400  std::string do_grouping() const override
401  {
402  return "\003"; // groups of 3 digits (this string is in octal format)
403  }
404 
405  };
406 
412  bool fexists(const std::string &filename);
413 
419  bool filename_is_url(const std::string &filename);
420 
438  std::string
439  read_and_distribute_file_content(const std::string &filename,
440  const MPI_Comm &comm);
441 
455  int
456  mkdirp(std::string pathname, const mode_t mode = 0755);
457 
468  void create_directory(const std::string &pathname,
469  const MPI_Comm &comm,
470  bool silent);
471 
476  namespace tk
477  {
481  class spline
482  {
483  public:
493  void set_points(const std::vector<double> &x,
494  const std::vector<double> &y,
495  const bool cubic_spline = true,
496  const bool monotone_spline = false);
500  double operator() (double x) const;
501 
502  private:
506  std::vector<double> m_x;
507 
514  std::vector<double> m_a, m_b, m_c, m_y;
515  };
516  }
517 
525  inline
526  void
527  extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
528  const unsigned int q,
529  std::vector<double> &composition_values_at_q_point)
530  {
531  Assert(q<composition_values.size(), ExcInternalError());
532  Assert(composition_values_at_q_point.size() > 0,
533  ExcInternalError());
534 
535  for (unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
536  {
537  Assert(composition_values[k].size() == composition_values_at_q_point.size(),
538  ExcInternalError());
539  composition_values_at_q_point[k] = composition_values[k][q];
540  }
541  }
542 
543  template <typename T>
544  inline
545  std::vector<T>
546  possibly_extend_from_1_to_N (const std::vector<T> &values,
547  const unsigned int N,
548  const std::string &id_text)
549  {
550  if (values.size() == 1)
551  {
552  return std::vector<T> (N, values[0]);
553  }
554  else if (values.size() == N)
555  {
556  return values;
557  }
558  else
559  {
560  // Non-specified behavior
561  AssertThrow(false,
562  ExcMessage("Length of " + id_text + " list must be " +
563  "either one or " + Utilities::to_string(N) +
564  ". Currently it is " + Utilities::to_string(values.size()) + "."));
565  }
566 
567  // This should never happen, but return an empty vector so the compiler
568  // will be happy
569  return std::vector<T> ();
570  }
571 
576  std::string
577  expand_ASPECT_SOURCE_DIR (const std::string &location);
578 
583  std::string parenthesize_if_nonempty (const std::string &s);
584 
589  bool has_unique_entries (const std::vector<std::string> &strings);
590 
591 
592 
614  double weighted_p_norm_average (const std::vector<double> &weights,
615  const std::vector<double> &values,
616  const double p);
617 
618 
644  template <typename T>
645  T derivative_of_weighted_p_norm_average (const double averaged_parameter,
646  const std::vector<double> &weights,
647  const std::vector<double> &values,
648  const std::vector<T> &derivatives,
649  const double p);
668  template <int dim>
669  double compute_spd_factor(const double eta,
670  const SymmetricTensor<2,dim> &strain_rate,
671  const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
672  const double SPD_safety_factor);
673 
677  template <int dim>
678  Point<dim> convert_array_to_point(const std::array<double,dim> &array);
679 
683  template <int dim>
684  std::array<double,dim> convert_point_to_array(const Point<dim> &point);
685 
693  class Operator
694  {
695  public:
700  {
706  replace_if_valid
707  };
708 
713  Operator();
714 
718  Operator(const operation op);
719 
724  double operator() (const double x, const double y) const;
725 
730  bool operator== (const operation op) const;
731 
732  private:
737  };
738 
743  std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
744 
748  const std::string get_model_operator_options();
749 
755  template <int dim>
756  SymmetricTensor<2,dim> nth_basis_for_symmetric_tensors (const unsigned int k);
757 
761  template <int dim>
763  {
764  public:
768  NaturalCoordinate(Point<dim> &position,
769  const GeometryModel::Interface<dim> &geometry_model);
770 
775  NaturalCoordinate(const std::array<double, dim> &coord,
776  const Utilities::Coordinates::CoordinateSystem &coord_system);
777 
782  std::array<double,dim> &get_coordinates();
783 
788  const std::array<double,dim> &get_coordinates() const;
789 
794  std::array<double,dim-1> get_surface_coordinates() const;
795 
800  double get_depth_coordinate() const;
801 
802  private:
808 
812  std::array<double,dim> coordinates;
813  };
814 
815 
834  template <int dim, typename VectorType>
835  void
836  project_cellwise(const Mapping<dim> &mapping,
837  const DoFHandler<dim> &dof_handler,
838  const unsigned int component_index,
839  const Quadrature<dim> &quadrature,
840  const std::function<void(
841  const typename DoFHandler<dim>::active_cell_iterator &,
842  const std::vector<Point<dim>> &,
843  std::vector<double> &)> &function,
844  VectorType &vec_result);
845 
865  void linear_solver_failed(const std::string &solver_name,
866  const std::string &function_name,
867  const std::vector<SolverControl> &solver_controls,
868  const std::exception &exc,
869  const MPI_Comm &mpi_communicator,
870  const std::string &output_filename = "");
871 
879  template <int dim>
881  {
882  public:
892  VectorFunctionFromVelocityFunctionObject (const unsigned int n_components,
893  const std::function<Tensor<1,dim> (const Point<dim> &)> &function_object);
894 
902  double value (const Point<dim> &p,
903  const unsigned int component = 0) const override;
904 
914  void vector_value (const Point<dim> &p,
915  Vector<double> &values) const override;
916 
917  private:
922  const std::function<Tensor<1,dim> (const Point<dim> &)> function_object;
923  };
924  }
925 }
926 
927 #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)
Point< dim > convert_array_to_point(const std::array< double, dim > &array)
std::vector< double > m_y
Definition: utilities.h:514
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)
double weighted_p_norm_average(const std::vector< double > &weights, const std::vector< double > &values, const double p)
std::vector< T > possibly_extend_from_1_to_N(const std::vector< T > &values, const unsigned int N, const std::string &id_text)
Definition: utilities.h:546
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:812
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:506
std::string expand_ASPECT_SOURCE_DIR(const std::string &location)
std::array< Tensor< 1, dim >, dim-1 > orthogonal_vectors(const Tensor< 1, dim > &v)
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 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:807
::Point< dim > spherical_to_cartesian_coordinates(const std::array< double, dim > &scoord)
std::string do_grouping() const override
Definition: utilities.h:400
std::array< double, 3 > cartesian_to_ellipsoidal_coordinates(const ::Point< 3 > &position, const double semi_major_axis_a, const double eccentricity)
std::array< double, dim > convert_point_to_array(const Point< dim > &point)
std::vector< Point< dim > > get_unit_support_points(const SimulatorAccess< dim > &simulator_access)
IndexSet extract_locally_active_dofs_with_component(const DoFHandler< dim > &dof_handler, const ComponentMask &component_mask)
bool fexists(const std::string &filename)
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 read_and_distribute_file_content(const std::string &filename, const MPI_Comm &comm)
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)
void create_directory(const std::string &pathname, const MPI_Comm &comm, bool silent)
std::string parenthesize_if_nonempty(const std::string &s)
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)
Definition: utilities.h:527
const std::string get_model_operator_options()
char do_thousands_sep() const override
Definition: utilities.h:395
bool filename_is_url(const std::string &filename)
SymmetricTensor< 2, dim > nth_basis_for_symmetric_tensors(const unsigned int k)
void linear_solver_failed(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="")
CoordinateSystem string_to_coordinate_system(const std::string &)
Definition: compat.h:88
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)