ASPECT
utilities.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2014 - 2019 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 
37 
38 
39 namespace aspect
40 {
41  template <int dim> class SimulatorAccess;
42 
43  namespace GeometryModel
44  {
45  template <int dim> class Interface;
46  }
47 
52  namespace Utilities
53  {
54  using namespace dealii;
55  using namespace ::Utilities;
56 
70  template <typename T>
71  std::vector<T>
72  possibly_extend_from_1_to_N (const std::vector<T> &values,
73  const unsigned int N,
74  const std::string &id_text);
75 
76 
77 
117  std::vector<double>
118  parse_map_to_double_array (const std::string &key_value_map,
119  const std::vector<std::string> &list_of_keys,
120  const bool has_background_field,
121  const std::string &property_name);
122 
136  template <typename T>
137  Table<2,T>
138  parse_input_table (const std::string &input_string,
139  const unsigned int n_rows,
140  const unsigned int n_columns,
141  const std::string &property_name);
142 
155  template <int dim>
156  std::vector<std::string>
157  expand_dimensional_variable_names (const std::vector<std::string> &var_declarations);
158 
159 
167  void split_by_block (const std::vector<types::global_dof_index> &dofs_per_block,
168  const IndexSet &whole_set,
169  std::vector<IndexSet> &partitioned);
170 
171 
178  template <int dim>
180  const ComponentMask &component_mask);
181 
182 
183 
184  namespace Coordinates
185  {
186 
192  template <int dim>
193  std::array<double,dim>
194  WGS84_coordinates(const Point<dim> &position);
195 
202  template <int dim>
203  std::array<double,dim>
205 
211  template <int dim>
212  Point<dim>
213  spherical_to_cartesian_coordinates(const std::array<double,dim> &scoord);
214 
220  template <int dim>
222  spherical_to_cartesian_vector(const Tensor<1,dim> &spherical_vector,
223  const Point<dim> &position);
224 
225 
231  template <int dim>
232  std::array<double,3>
234  const double semi_major_axis_a,
235  const double eccentricity);
236 
241  template <int dim>
242  Point<3>
243  ellipsoidal_to_cartesian_coordinates(const std::array<double,3> &phi_theta_d,
244  const double semi_major_axis_a,
245  const double eccentricity);
246 
247 
254  string_to_coordinate_system (const std::string &);
255  }
256 
257 
262  template <int dim>
263  bool
264  polygon_contains_point(const std::vector<Point<2> > &point_list,
265  const ::Point<2> &point);
266 
272  template <int dim>
273  double
274  signed_distance_to_polygon(const std::vector<Point<2> > &point_list,
275  const ::Point<2> &point);
276 
277 
284  double
285  distance_to_line(const std::array<::Point<2>,2 > &point_list,
286  const ::Point<2> &point);
287 
295  template <int dim>
296  std::array<Tensor<1,dim>,dim-1>
298 
335  std::pair<double,double> real_spherical_harmonic( unsigned int l, // degree
336  unsigned int m, // order
337  double theta, // colatitude (radians)
338  double phi ); // longitude (radians)
339 
343  struct ThousandSep : std::numpunct<char>
344  {
345  protected:
346  virtual char do_thousands_sep() const
347  {
348  return ',';
349  }
350  virtual std::string do_grouping() const
351  {
352  return "\003"; // groups of 3 digits (this string is in octal format)
353  }
354 
355  };
356 
362  bool fexists(const std::string &filename);
363 
374  std::string
375  read_and_distribute_file_content(const std::string &filename,
376  const MPI_Comm &comm);
377 
391  int
392  mkdirp(std::string pathname, const mode_t mode = 0755);
393 
404  void create_directory(const std::string &pathname,
405  const MPI_Comm &comm,
406  bool silent);
407 
412  namespace tk
413  {
417  class spline
418  {
419  public:
429  void set_points(const std::vector<double> &x,
430  const std::vector<double> &y,
431  const bool cubic_spline = true,
432  const bool monotone_spline = false);
436  double operator() (double x) const;
437 
438  private:
442  std::vector<double> m_x;
443 
450  std::vector<double> m_a, m_b, m_c, m_y;
451  };
452  }
453 
461  inline
462  void
463  extract_composition_values_at_q_point (const std::vector<std::vector<double> > &composition_values,
464  const unsigned int q,
465  std::vector<double> &composition_values_at_q_point)
466  {
467  Assert(q<composition_values.size(), ExcInternalError());
468  Assert(composition_values_at_q_point.size() > 0,
469  ExcInternalError());
470 
471  for (unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
472  {
473  Assert(composition_values[k].size() == composition_values_at_q_point.size(),
474  ExcInternalError());
475  composition_values_at_q_point[k] = composition_values[k][q];
476  }
477  }
478 
479  template <typename T>
480  inline
481  std::vector<T>
482  possibly_extend_from_1_to_N (const std::vector<T> &values,
483  const unsigned int N,
484  const std::string &id_text)
485  {
486  if (values.size() == 1)
487  {
488  return std::vector<T> (N, values[0]);
489  }
490  else if (values.size() == N)
491  {
492  return values;
493  }
494  else
495  {
496  // Non-specified behavior
497  AssertThrow(false,
498  ExcMessage("Length of " + id_text + " list must be " +
499  "either one or " + Utilities::to_string(N) +
500  ". Currently it is " + Utilities::to_string(values.size()) + "."));
501  }
502 
503  // This should never happen, but return an empty vector so the compiler
504  // will be happy
505  return std::vector<T> ();
506  }
507 
512  std::string
513  expand_ASPECT_SOURCE_DIR (const std::string &location);
514 
519  std::string parenthesize_if_nonempty (const std::string &s);
520 
525  bool has_unique_entries (const std::vector<std::string> &strings);
526 
543  template <int dim>
545  {
546  public:
556  AsciiDataLookup(const unsigned int components,
557  const double scale_factor);
558 
565  explicit AsciiDataLookup(const double scale_factor);
566 
572  void
573  load_file(const std::string &filename,
574  const MPI_Comm &communicator);
575 
586  double
587  get_data(const Point<dim> &position,
588  const unsigned int component) const;
589 
596  std::vector<std::string>
597  get_column_names() const;
598 
604  bool
605  has_equidistant_coordinates() const;
606 
612  unsigned int
613  get_column_index_from_name(const std::string &column_name) const;
614 
620  std::string
621  get_column_name_from_index(const unsigned int column_index) const;
622 
626  double get_maximum_component_value(const unsigned int component) const;
627 
628  private:
632  unsigned int components;
633 
639  std::vector<std::string> data_component_names;
640 
646  std::vector<std::unique_ptr<Function<dim>>> data;
647 
651  std::array<std::vector<double>,dim> coordinate_values;
652 
656  std::vector<double> maximum_component_value;
657 
661  std::array<std::pair<double,double>,dim> grid_extent;
662 
667 
672  const double scale_factor;
673 
679 
685  compute_table_indices(const unsigned int i) const;
686 
687  };
688 
693  template <int dim>
695  {
696  public:
700  AsciiDataBase();
701 
705  static
706  void
708  const std::string &default_directory,
709  const std::string &default_filename,
710  const std::string &subsection_name = "Ascii data model");
711 
715  void
716  parse_parameters (ParameterHandler &prm,
717  const std::string &subsection_name = "Ascii data model");
718 
722  std::string data_directory;
723 
729  std::string data_file_name;
730 
737  double scale_factor;
738  };
739 
744  template <int dim>
746  {
747  public:
752 
757  virtual
758  void
759  initialize (const std::set<types::boundary_id> &boundary_ids,
760  const unsigned int components);
761 
768  void
769  update();
770 
774  double
775  get_data_component (const types::boundary_id boundary_indicator,
776  const Point<dim> &position,
777  const unsigned int component) const;
778 
782  double
783  get_maximum_component_value (const types::boundary_id boundary_indicator,
784  const unsigned int component) const;
785 
789  static
790  void
792  const std::string &default_directory,
793  const std::string &default_filename,
794  const std::string &subsection_name = "Ascii data model");
795 
799  void
800  parse_parameters (ParameterHandler &prm,
801  const std::string &subsection_name = "Ascii data model");
802 
803  protected:
804 
818  std::array<unsigned int,dim-1>
819  get_boundary_dimensions (const types::boundary_id boundary_id) const;
820 
826 
835 
841 
849 
855 
860  double time_weight;
861 
868 
873  std::map<types::boundary_id,
874  std::unique_ptr<aspect::Utilities::AsciiDataLookup<dim-1> > > lookups;
875 
879  std::map<types::boundary_id,
881 
885  void
886  update_data (const types::boundary_id boundary_id,
887  const bool reload_both_files);
888 
893  void
894  end_time_dependence ();
895 
899  std::string
900  create_filename (const int filenumber,
901  const types::boundary_id boundary_id) const;
902  };
903 
904 
905 
910  template <int dim>
912  {
913  public:
918 
923  virtual
924  void
925  initialize (const unsigned int components);
926 
927 
931  double
932  get_data_component (const Point<dim> &position,
933  const unsigned int component) const;
934 
935  protected:
940  std::unique_ptr<aspect::Utilities::AsciiDataLookup<dim> > lookup;
941  };
942 
943 
948  template <int dim>
950  {
951  public:
956 
961  virtual
962  void
963  initialize (const unsigned int components);
964 
965 
969  double
970  get_data_component (const Point<dim> &position,
971  const unsigned int component) const;
972 
973 
977  static
978  void
980  const std::string &default_directory,
981  const std::string &default_filename,
982  const std::string &subsection_name = "Ascii data model");
983 
987  void
988  parse_parameters (ParameterHandler &prm,
989  const std::string &subsection_name = "Ascii data model");
990 
991  protected:
996  std::vector<std::unique_ptr<aspect::Utilities::AsciiDataLookup<dim-1> >> lookups;
997 
998  private:
999 
1003  std::string data_directory;
1004 
1008  std::vector<std::string> data_file_names;
1009 
1014 
1019 
1020 
1021  };
1022 
1023 
1027  template <int dim>
1029  {
1030  public:
1034  AsciiDataProfile();
1035 
1040  virtual
1041  void
1042  initialize (const MPI_Comm &communicator);
1043 
1044 
1048  double
1049  get_data_component (const Point<1> &position,
1050  const unsigned int component) const;
1051 
1058  std::vector<std::string>
1059  get_column_names() const;
1060 
1066  unsigned int
1067  get_column_index_from_name(const std::string &column_name) const;
1068 
1074  unsigned int
1075  maybe_get_column_index_from_name(const std::string &column_name) const;
1076 
1082  std::string
1083  get_column_name_from_index(const unsigned int column_index) const;
1084  protected:
1089  std::unique_ptr<aspect::Utilities::AsciiDataLookup<1> > lookup;
1090  };
1091 
1092 
1114  double weighted_p_norm_average (const std::vector<double> &weights,
1115  const std::vector<double> &values,
1116  const double p);
1117 
1118 
1144  template <typename T>
1145  T derivative_of_weighted_p_norm_average (const double averaged_parameter,
1146  const std::vector<double> &weights,
1147  const std::vector<double> &values,
1148  const std::vector<T> &derivatives,
1149  const double p);
1168  template <int dim>
1169  double compute_spd_factor(const double eta,
1171  const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
1172  const double SPD_safety_factor);
1173 
1177  template <int dim>
1178  Point<dim> convert_array_to_point(const std::array<double,dim> &array);
1179 
1183  template <int dim>
1184  std::array<double,dim> convert_point_to_array(const Point<dim> &point);
1185 
1193  class Operator
1194  {
1195  public:
1200  {
1206  replace_if_valid
1207  };
1208 
1213  Operator();
1214 
1218  Operator(const operation op);
1219 
1224  double operator() (const double x, const double y) const;
1225 
1230  bool operator== (const operation op) const;
1231 
1232  private:
1237  };
1238 
1243  std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
1244 
1248  const std::string get_model_operator_options();
1249 
1255  template <int dim>
1257 
1258  /*
1259  * A class that represents a point in a chosen coordinate system.
1260  */
1261  template <int dim>
1263  {
1264  public:
1268  NaturalCoordinate(Point<dim> &position,
1269  const GeometryModel::Interface<dim> &geometry_model);
1270 
1275  NaturalCoordinate(const std::array<double, dim> &coord,
1276  const Utilities::Coordinates::CoordinateSystem &coord_system);
1277 
1282  std::array<double,dim> &get_coordinates();
1283 
1288  std::array<double,dim-1> get_surface_coordinates() const;
1289 
1294  double get_depth_coordinate() const;
1295 
1296  private:
1302 
1306  std::array<double,dim> coordinates;
1307  };
1308  }
1309 }
1310 
1311 #endif
std::vector< double > parse_map_to_double_array(const std::string &key_value_map, const std::vector< std::string > &list_of_keys, const bool has_background_field, const std::string &property_name)
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:450
Tensor< 1, dim > spherical_to_cartesian_vector(const Tensor< 1, dim > &spherical_vector, const Point< dim > &position)
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:482
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)
virtual std::string do_grouping() const
Definition: utilities.h:350
Point< 3 > ellipsoidal_to_cartesian_coordinates(const std::array< double, 3 > &phi_theta_d, const double semi_major_axis_a, const double eccentricity)
std::array< double, dim > coordinates
Definition: utilities.h:1306
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:442
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::array< double, 3 > cartesian_to_ellipsoidal_coordinates(const Point< 3 > &position, const double semi_major_axis_a, const double eccentricity)
std::unique_ptr< aspect::Utilities::AsciiDataLookup< 1 > > lookup
Definition: utilities.h:1089
int mkdirp(std::string pathname, const mode_t mode=0755)
Utilities::Coordinates::CoordinateSystem coordinate_system
Definition: utilities.h:1301
virtual char do_thousands_sep() const
Definition: utilities.h:346
#define AssertThrow(cond, exc)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
std::array< double, dim > convert_point_to_array(const Point< dim > &point)
std::unique_ptr< aspect::Utilities::AsciiDataLookup< dim > > lookup
Definition: utilities.h:940
void declare_parameters(ParameterHandler &prm)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
std::array< double, dim > WGS84_coordinates(const Point< dim > &position)
IndexSet extract_locally_active_dofs_with_component(const DoFHandler< dim > &dof_handler, const ComponentMask &component_mask)
std::vector< std::string > data_file_names
Definition: utilities.h:1008
bool fexists(const std::string &filename)
std::vector< std::unique_ptr< aspect::Utilities::AsciiDataLookup< dim-1 > > > lookups
Definition: utilities.h:996
std::array< std::pair< double, double >, dim > grid_extent
Definition: utilities.h:661
double signed_distance_to_polygon(const std::vector< Point< 2 > > &point_list, const ::Point< 2 > &point)
std::array< std::vector< double >, dim > coordinate_values
Definition: utilities.h:651
std::vector< Operator > create_model_operator_list(const std::vector< std::string > &operator_names)
std::string read_and_distribute_file_content(const std::string &filename, const MPI_Comm &comm)
void create_directory(const std::string &pathname, const MPI_Comm &comm, bool silent)
std::string parenthesize_if_nonempty(const std::string &s)
Definition: compat.h:37
const std::string get_model_operator_options()
std::array< double, dim > cartesian_to_spherical_coordinates(const Point< dim > &position)
bool polygon_contains_point(const std::vector< Point< 2 > > &point_list, const ::Point< 2 > &point)
TableIndices< dim > table_points
Definition: utilities.h:666
std::vector< std::string > data_component_names
Definition: utilities.h:639
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:463
SymmetricTensor< 2, dim > nth_basis_for_symmetric_tensors(const unsigned int k)
void split_by_block(const std::vector< types::global_dof_index > &dofs_per_block, const IndexSet &whole_set, std::vector< IndexSet > &partitioned)
std::map< types::boundary_id, std::unique_ptr< aspect::Utilities::AsciiDataLookup< dim-1 > > > lookups
Definition: utilities.h:874
std::map< types::boundary_id, std::unique_ptr< aspect::Utilities::AsciiDataLookup< dim-1 > > > old_lookups
Definition: utilities.h:880
CoordinateSystem string_to_coordinate_system(const std::string &)
std::vector< std::unique_ptr< Function< dim > > > data
Definition: utilities.h:646
unsigned int boundary_id
bool has_unique_entries(const std::vector< std::string > &strings)
Point< dim > spherical_to_cartesian_coordinates(const std::array< double, dim > &scoord)
std::vector< double > maximum_component_value
Definition: utilities.h:656
static ::ExceptionBase & ExcInternalError()