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  AsciiDataLookup(const double scale_factor);
566 
572  void
573  load_file(const std::string &filename,
574  const MPI_Comm &communicator);
575 
584  double
585  get_data(const Point<dim> &position,
586  const unsigned int component) const;
587 
594  std::vector<std::string>
595  get_column_names() const;
596 
602  bool
603  has_equidistant_coordinates() const;
604 
610  unsigned int
611  get_column_index_from_name(const std::string &column_name) const;
612 
618  std::string
619  get_column_name_from_index(const unsigned int column_index) const;
620 
624  double get_maximum_component_value(const unsigned int component) const;
625 
626  private:
630  unsigned int components;
631 
637  std::vector<std::string> data_component_names;
638 
644  std::vector<std::unique_ptr<Function<dim>>> data;
645 
649  std::array<std::vector<double>,dim> coordinate_values;
650 
654  std::vector<double> maximum_component_value;
655 
659  std::array<std::pair<double,double>,dim> grid_extent;
660 
665 
670  const double scale_factor;
671 
677 
683  compute_table_indices(const unsigned int i) const;
684 
685  };
686 
691  template <int dim>
693  {
694  public:
698  AsciiDataBase();
699 
703  static
704  void
706  const std::string &default_directory,
707  const std::string &default_filename,
708  const std::string &subsection_name = "Ascii data model");
709 
713  void
714  parse_parameters (ParameterHandler &prm,
715  const std::string &subsection_name = "Ascii data model");
716 
720  std::string data_directory;
721 
727  std::string data_file_name;
728 
735  double scale_factor;
736  };
737 
742  template <int dim>
744  {
745  public:
750 
755  virtual
756  void
757  initialize (const std::set<types::boundary_id> &boundary_ids,
758  const unsigned int components);
759 
766  void
767  update();
768 
772  double
773  get_data_component (const types::boundary_id boundary_indicator,
774  const Point<dim> &position,
775  const unsigned int component) const;
776 
780  double
781  get_maximum_component_value (const types::boundary_id boundary_indicator,
782  const unsigned int component) const;
783 
787  static
788  void
790  const std::string &default_directory,
791  const std::string &default_filename,
792  const std::string &subsection_name = "Ascii data model");
793 
797  void
798  parse_parameters (ParameterHandler &prm,
799  const std::string &subsection_name = "Ascii data model");
800 
801  protected:
802 
816  std::array<unsigned int,dim-1>
817  get_boundary_dimensions (const types::boundary_id boundary_id) const;
818 
824 
833 
839 
847 
853 
858  double time_weight;
859 
866 
871  std::map<types::boundary_id,
872  std::unique_ptr<aspect::Utilities::AsciiDataLookup<dim-1> > > lookups;
873 
877  std::map<types::boundary_id,
879 
883  void
884  update_data (const types::boundary_id boundary_id,
885  const bool reload_both_files);
886 
891  void
892  end_time_dependence ();
893 
897  std::string
898  create_filename (const int timestep,
899  const types::boundary_id boundary_id) const;
900  };
901 
902 
903 
908  template <int dim>
910  {
911  public:
916 
921  virtual
922  void
923  initialize (const unsigned int components);
924 
925 
929  double
930  get_data_component (const Point<dim> &position,
931  const unsigned int component) const;
932 
933  protected:
938  std::unique_ptr<aspect::Utilities::AsciiDataLookup<dim> > lookup;
939  };
940 
941 
942 
946  template <int dim>
948  {
949  public:
954 
959  virtual
960  void
961  initialize (const MPI_Comm &communicator);
962 
963 
967  double
968  get_data_component (const Point<1> &position,
969  const unsigned int component) const;
970 
977  std::vector<std::string>
978  get_column_names() const;
979 
985  unsigned int
986  get_column_index_from_name(const std::string &column_name) const;
987 
993  unsigned int
994  maybe_get_column_index_from_name(const std::string &column_name) const;
995 
1001  std::string
1002  get_column_name_from_index(const unsigned int column_index) const;
1003  protected:
1008  std::unique_ptr<aspect::Utilities::AsciiDataLookup<1> > lookup;
1009  };
1010 
1011 
1033  double weighted_p_norm_average (const std::vector<double> &weights,
1034  const std::vector<double> &values,
1035  const double p);
1036 
1037 
1063  template <typename T>
1064  T derivative_of_weighted_p_norm_average (const double averaged_parameter,
1065  const std::vector<double> &weights,
1066  const std::vector<double> &values,
1067  const std::vector<T> &derivatives,
1068  const double p);
1087  template <int dim>
1088  double compute_spd_factor(const double eta,
1090  const SymmetricTensor<2,dim> &dviscosities_dstrain_rate,
1091  const double SPD_safety_factor);
1092 
1096  template <int dim>
1097  Point<dim> convert_array_to_point(const std::array<double,dim> &array);
1098 
1102  template <int dim>
1103  std::array<double,dim> convert_point_to_array(const Point<dim> &point);
1104 
1112  class Operator
1113  {
1114  public:
1119  {
1124  maximum
1125  };
1126 
1131  Operator();
1132 
1136  Operator(const operation op);
1137 
1142  double operator() (const double x, const double y) const;
1143 
1148  bool operator== (const operation op) const;
1149 
1150  private:
1155  };
1156 
1161  std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
1162 
1168  template <int dim>
1170 
1171  /*
1172  * A class that represents a point in a chosen coordinate system.
1173  */
1174  template <int dim>
1176  {
1177  public:
1181  NaturalCoordinate(Point<dim> &position,
1182  const GeometryModel::Interface<dim> &geometry_model);
1183 
1188  NaturalCoordinate(const std::array<double, dim> &coord,
1189  const Utilities::Coordinates::CoordinateSystem &coord_system);
1190 
1195  std::array<double,dim> &get_coordinates();
1196 
1201  std::array<double,dim-1> get_surface_coordinates() const;
1202 
1207  double get_depth_coordinate() const;
1208 
1209  private:
1215 
1219  std::array<double,dim> coordinates;
1220  };
1221  }
1222 }
1223 
1224 #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:1219
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:1008
int mkdirp(std::string pathname, const mode_t mode=0755)
Utilities::Coordinates::CoordinateSystem coordinate_system
Definition: utilities.h:1214
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:938
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)
bool fexists(const std::string &filename)
std::array< std::pair< double, double >, dim > grid_extent
Definition: utilities.h:659
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:649
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:38
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:664
std::vector< std::string > data_component_names
Definition: utilities.h:637
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:872
std::map< types::boundary_id, std::unique_ptr< aspect::Utilities::AsciiDataLookup< dim-1 > > > old_lookups
Definition: utilities.h:878
CoordinateSystem string_to_coordinate_system(const std::string &)
std::vector< std::unique_ptr< Function< dim > > > data
Definition: utilities.h:644
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:654
static ::ExceptionBase & ExcInternalError()