ASPECT
Namespaces | Classes | Functions
aspect::Utilities Namespace Reference

Namespaces

 Coordinates
 
 tk
 

Classes

class  AsciiDataBase
 
class  AsciiDataBoundary
 
class  AsciiDataInitial
 
class  AsciiDataLookup
 
class  AsciiDataProfile
 
class  NaturalCoordinate
 
class  Operator
 
struct  ThousandSep
 

Functions

template<typename T >
std::vector< Tpossibly_extend_from_1_to_N (const std::vector< T > &values, const unsigned int N, const std::string &id_text)
 
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)
 
template<typename T >
Table< 2, Tparse_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)
 
void split_by_block (const std::vector< types::global_dof_index > &dofs_per_block, const IndexSet &whole_set, std::vector< IndexSet > &partitioned)
 
template<int dim>
IndexSet extract_locally_active_dofs_with_component (const DoFHandler< dim > &dof_handler, const ComponentMask &component_mask)
 
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)
 
std::pair< double, double > real_spherical_harmonic (unsigned int l, unsigned int m, double theta, double phi)
 
bool fexists (const std::string &filename)
 
std::string read_and_distribute_file_content (const std::string &filename, 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, bool silent)
 
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 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)
 
std::vector< Operatorcreate_model_operator_list (const std::vector< std::string > &operator_names)
 
template<int dim>
SymmetricTensor< 2, dim > nth_basis_for_symmetric_tensors (const unsigned int k)
 

Detailed Description

A namespace for utility functions that might be used in many different places to prevent code duplication.

Function Documentation

§ possibly_extend_from_1_to_N()

template<typename T >
std::vector< T > aspect::Utilities::possibly_extend_from_1_to_N ( const std::vector< T > &  values,
const unsigned int  N,
const std::string &  id_text 
)
inline

Given an array values, consider three cases:

  • If it has size N, return the original array.
  • If it has size one, return an array of size N where all elements are equal to the one element of value.
  • If it has any other size, throw an exception that uses id_text as an identifying string.

This function is typically used for parameter lists that can either contain different values for each of a set of objects (e.g., for each compositional field), or contain a single value that is then used for each object.

Definition at line 482 of file utilities.h.

References AssertThrow, ExcMessage(), expand_ASPECT_SOURCE_DIR(), has_unique_entries(), parenthesize_if_nonempty(), and Utilities::to_string().

§ parse_map_to_double_array()

std::vector<double> aspect::Utilities::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 
)

This function takes a string argument that is interpreted as a map of the form "key1 : value1, key2 : value2, etc", and then parses it to return a vector of these values where the values are ordered in the same order as a given set of keys.

This function also considers a number of special cases:

  • If the input string consists of only a comma separated set of values "value1, value2, value3, ..." (i.e., without the "keyx :" part), then the input string is interpreted as if it had had the form "key1 : value1, key2 : value2, ..." where "key1", "key2", ... are exactly the keys provided by the list_of_keys in the same order as provided. In this situation, if a background field is required, the background value is assigned to the first element of the output vector.
  • Whether or not a background field is required depends on the parameter being parsed. Requiring a background field increases the size of the output vector by 1. For example, some Material models require background fields, but input files may not.
  • Three special keys are recognized: all –> Assign the associated value to all fields. Only one value is allowed in this case. background –> Assign associated value to the background.
Parameters
[in]key_value_mapThe string representation of the map to be parsed.
[in]list_of_keysA list of valid key names that are allowed to appear in the map. The order of these keys determines the order of values that are returned by this function.
[in]has_background_fieldIf true, expect N+1 values and allow setting of the background using the key "background".
[in]property_nameA name that identifies the type of property that is being parsed by this function and that is used in generating error messages if the map does not conform to the expected format.
Returns
A vector of values that are parsed from the map, provided in the order in which the keys appear in the list_of_keys argument.

§ parse_input_table()

template<typename T >
Table<2,T> aspect::Utilities::parse_input_table ( const std::string &  input_string,
const unsigned int  n_rows,
const unsigned int  n_columns,
const std::string &  property_name 
)

This function takes a string argument that is assumed to represent an input table, in which each row is separated by semicolons, and each column separated by commas. The function returns the parsed entries as a table. In addition this function utilizes the possibly_extend_from_1_to_N() function to accept inputs with only a single column/row, which will be extended to n_rows or n_columns respectively, to allow abbreviating the input string (e.g. you can provide a single value instead of n_rows by n_columns identical values). This function can for example be used by material models to read densities for different compositions and different phases for each composition.

§ expand_dimensional_variable_names()

template<int dim>
std::vector<std::string> aspect::Utilities::expand_dimensional_variable_names ( const std::vector< std::string > &  var_declarations)

Given a vector var_declarations expand any entries of the form vector(str) or tensor(str) to sublists with component names of the form str_x, str_y, str_z or str_xx, str_xy... for the correct dimension value.

This function is to be used for expanding lists of variable names where one or more such variable is actually intended to be a list of components.

Returns the generated list of variable names

§ split_by_block()

void aspect::Utilities::split_by_block ( const std::vector< types::global_dof_index > &  dofs_per_block,
const IndexSet whole_set,
std::vector< IndexSet > &  partitioned 
)

Split the set of DoFs (typically locally owned or relevant) in whole_set into blocks given by the dofs_per_block structure.

The numbers of dofs per block need to add up to the size of the index space described by whole_set.

§ extract_locally_active_dofs_with_component()

template<int dim>
IndexSet aspect::Utilities::extract_locally_active_dofs_with_component ( const DoFHandler< dim > &  dof_handler,
const ComponentMask component_mask 
)

Returns an IndexSet that contains all locally active DoFs that belong to the given component_mask.

This function should be moved into deal.II at some point.

§ polygon_contains_point()

template<int dim>
bool aspect::Utilities::polygon_contains_point ( const std::vector< Point< 2 > > &  point_list,
const ::Point< 2 > &  point 
)

Given a 2d point and a list of points which form a polygon, computes if the point falls within the polygon.

§ signed_distance_to_polygon()

template<int dim>
double aspect::Utilities::signed_distance_to_polygon ( const std::vector< Point< 2 > > &  point_list,
const ::Point< 2 > &  point 
)

Given a 2d point and a list of points which form a polygon, compute the smallest distance of the point to the polygon. The sign is negative for points outside of the polygon and positive for points inside the polygon.

§ distance_to_line()

double aspect::Utilities::distance_to_line ( const std::array<::Point< 2 >, 2 > &  point_list,
const ::Point< 2 > &  point 
)

Given a 2d point and a list of two points that define a line, compute the smallest distance of the point to the line segment. When the point's perpendicular base does not lie on the line segment, the smallest distance to the segment's end points is calculated.

§ orthogonal_vectors()

template<int dim>
std::array<Tensor<1,dim>,dim-1> aspect::Utilities::orthogonal_vectors ( const Tensor< 1, dim > &  v)

Given a vector v in dim dimensional space, return a set of (dim-1) vectors that are orthogonal to v and to each other. The length of each of these vectors equals that of the original vector v to ensure that the resulting set of vectors represents a well-conditioned basis.

§ real_spherical_harmonic()

std::pair<double,double> aspect::Utilities::real_spherical_harmonic ( unsigned int  l,
unsigned int  m,
double  theta,
double  phi 
)

A function for evaluating real spherical harmonics. It takes the degree (l) and the order (m) of the spherical harmonic, where $l \geq 0$ and $0 \leq m \leq l$. It also takes the colatitude (theta) and longitude (phi), which are in radians.

There are an unfortunate number of normalization conventions in existence for spherical harmonics. Here we use fully normalized spherical harmonics including the Condon-Shortley phase. This corresponds to the definitions given in equations B.72 and B.99-B.102 in Dahlen and Tromp (1998, ISBN: 9780691001241). The functional form of the real spherical harmonic is given by

\[ Y_{lm}(\theta, \phi) = \sqrt{2} X_{l \left| m \right| }(\theta) \cos m \phi \qquad \mathrm{if} \qquad -l \le m < 0 \]

\[ Y_{lm}(\theta, \phi) = X_{l 0 }(\theta) \qquad \mathrm{if} \qquad m = 0 \]

\[ Y_{lm}(\theta, \phi) = \sqrt{2} X_{lm}(\theta) \sin m \phi \qquad \mathrm{if} \qquad 0< m \le m \]

where $X_{lm}( \theta )$ is an associated Legendre function.

In practice it is often convenient to compute the sine ( $-l \le m < 0$) and cosine ( $0 < m \le l$) variants of the real spherical harmonic at the same time. That is the approach taken here, where we return a pair of numbers, the first corresponding the cosine part and the second corresponding to the sine part. Given this, it is no longer necessary to distinguish between positive and negative $m$, so this function only accepts $ m \ge 0$. For $m = 0$, there is only one part, which is stored in the first entry of the pair.

Note
This function uses the Boost spherical harmonics implementation internally, which is not designed for very high order (> 100) spherical harmonics computation. If you use spherical harmonics of a high order be sure to confirm the accuracy first. For more information, see: http://www.boost.org/doc/libs/1_49_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_poly/sph_harm.html

§ fexists()

bool aspect::Utilities::fexists ( const std::string &  filename)

Checks whether a file named filename exists and is readable.

Parameters
filenameFile to check existence

Referenced by aspect::Utilities::ThousandSep::do_grouping().

§ read_and_distribute_file_content()

std::string aspect::Utilities::read_and_distribute_file_content ( const std::string &  filename,
const MPI_Comm &  comm 
)

Reads the content of the ascii file filename on process 0 and distributes the content by MPI_Bcast to all processes. The function returns the content of the file on all processes.

Parameters
[in]filenameThe name of the ascii file to load.
[in]commThe MPI communicator in which the content is distributed.
Returns
A string which contains the data in filename.

Referenced by aspect::Utilities::ThousandSep::do_grouping().

§ mkdirp()

int aspect::Utilities::mkdirp ( std::string  pathname,
const mode_t  mode = 0755 
)

Creates a path as if created by the shell command "mkdir -p", therefore generating directories from the highest to the lowest level if they are not already existing.

Parameters
pathnameString that contains the path to create. '/' is used as directory separator.
modePermissions (mode bits) of the created directories. See the documentation of the chmod() command for more information.
Returns
The function returns the error value of the last mkdir call inside. It returns zero on success. See the man page of mkdir() for more information.

Referenced by aspect::Utilities::ThousandSep::do_grouping().

§ create_directory()

void aspect::Utilities::create_directory ( const std::string &  pathname,
const MPI_Comm &  comm,
bool  silent 
)

Create directory pathname, optionally printing a message.

Parameters
pathnameString that contains path to create. '/' is used as directory separator.
commMPI communicator, used to limit creation of directory to processor 0.
silentPrint a nicely formatted message on processor 0 if set to true.

Referenced by aspect::Utilities::ThousandSep::do_grouping().

§ extract_composition_values_at_q_point()

void aspect::Utilities::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 
)
inline

Extract the compositional values at a single quadrature point with index q from composition_values, which is indexed by compositional index and quadrature point, and write them into composition_values_at_q_point. In other words, this extracts composition_values[i][q] for all i.

Definition at line 463 of file utilities.h.

References Assert, and ExcInternalError().

§ expand_ASPECT_SOURCE_DIR()

std::string aspect::Utilities::expand_ASPECT_SOURCE_DIR ( const std::string &  location)

Replace the string $ASPECT_SOURCE_DIR in location by the current source directory of ASPECT and return the resulting string.

Referenced by possibly_extend_from_1_to_N().

§ parenthesize_if_nonempty()

std::string aspect::Utilities::parenthesize_if_nonempty ( const std::string &  s)

Given a string s, return it in the form ' ("s")' if nonempty. Otherwise just return the empty string itself.

Referenced by possibly_extend_from_1_to_N().

§ has_unique_entries()

bool aspect::Utilities::has_unique_entries ( const std::vector< std::string > &  strings)

Returns if a vector of strings strings only contains unique entries.

Referenced by possibly_extend_from_1_to_N().

§ weighted_p_norm_average()

double aspect::Utilities::weighted_p_norm_average ( const std::vector< double > &  weights,
const std::vector< double > &  values,
const double  p 
)

This function computes the weighted average $\bar y$ of $y_i$ for a weighted p norm. This leads for a general p to: $\bar y = \left(\frac{\sum_{i=1}^k w_i y_i^p}{\sum_{i=1}^k w_i}\right)^{\frac{1}{p}}$. When p = 0 we take the geometric average: $\bar y = \exp\left(\frac{\sum_{i=1}^k w_i \log\left(y_i\right)}{\sum_{i=1}^k w_i}\right)$, and when $p \le -1000$ or $p \ge 1000$ we take the minimum and maximum norm respectively. This means that the smallest and largest value is respectively taken taken.

This function has been set up to be very tolerant to strange values, such as negative weights. The only things we require in for the general p is that the sum of the weights and the sum of the weights times the values to the power p may not be smaller or equal to zero. Furthermore, when a value is zero, the exponent is smaller then one and the correspondent weight is non-zero, this corresponds to no resistance in a parallel system. This means that this 'path' will be followed, and we return zero.

The implemented special cases (which are minimum (p <= -1000), harmonic average (p = -1), geometric average (p = 0), arithmetic average (p = 1), quadratic average (RMS) (p = 2), cubic average (p = 3) and maximum (p >= 1000) ) is, except for the harmonic and quadratic averages even more tolerant of negative values, because they only require the sum of weights to be non-zero.

§ derivative_of_weighted_p_norm_average()

template<typename T >
T aspect::Utilities::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 
)

This function computes the derivative ( $\frac{\partial\bar y}{\partial x}$) of an average of the values $y_i(x)$ with regard to $x$, using $\frac{\partial y_i(x)}{\partial x}$. This leads for a general p to: $\frac{\partial\bar y}{\partial x} = \frac{1}{p}\left(\frac{\sum_{i=1}^k w_i y_i^p}{\sum_{i=1}^k w_i}\right)^{\frac{1}{p}-1} \frac{\sum_{i=1}^k w_i p y_i^{p-1} y'_i}{\sum_{i=1}^k w_i}$. When p = 0 we take the geometric average as a reference, which results in: $\frac{\partial\bar y}{\partial x} = \exp\left(\frac{\sum_{i=1}^k w_i \log\left(y_i\right)}{\sum_{i=1}^k w_i}\right) \frac{\sum_{i=1}^k\frac{w_i y'_i}{y_i}}{\sum_{i=1}^k w_i}$ and when $p \le -1000$ or $p \ge 1000$ we take the min and max norm respectively. This means that the derivative is taken which has the min/max value.

This function has, like the function weighted_p_norm_average been set up to be very tolerant to strange values, such as negative weights. The only things we require in for the general p is that the sum of the weights and the sum of the weights times the values to the power p may not be smaller or equal to zero. Furthermore, when a value is zero, the exponent is smaller then one and the correspondent weight is non-zero, this corresponds to no resistance in a parallel system. This means that this 'path' will be followed, and we return the corresponding derivative.

The implemented special cases (which are minimum (p <= -1000), harmonic average (p = -1), geometric average (p = 0), arithmetic average (p = 1), and maximum (p >= 1000) ) is, except for the harmonic average even more tolerant of negative values, because they only require the sum of weights to be non-zero.

§ compute_spd_factor()

template<int dim>
double aspect::Utilities::compute_spd_factor ( const double  eta,
const SymmetricTensor< 2, dim > &  strain_rate,
const SymmetricTensor< 2, dim > &  dviscosities_dstrain_rate,
const double  SPD_safety_factor 
)

This function computes a factor which can be used to make sure that the Jacobian remains positive definite.

The goal of this function is to find a factor $\alpha$ so that $2\eta(\varepsilon(\mathbf u)) I \otimes I + \alpha\left[a \otimes b + b \otimes a\right]$ remains a positive definite matrix. Here, $a=\varepsilon(\mathbf u)$ is the strain_rate and $b=\frac{\partial\eta(\varepsilon(\mathbf u),p)}{\partial \varepsilon}$ is the derivative of the viscosity with respect to the strain rate and is given by dviscosities_dstrain_rate. Since the viscosity $\eta$ must be positive, there is always a value of $\alpha$ (possibly small) so that the result is a positive definite matrix. In the best case, we want to choose $\alpha=1$ because that corresponds to the full Newton step, and so the function never returns anything larger than one.

The factor is defined by: $\frac{2\eta(\varepsilon(\mathbf u))}{\left[1-\frac{b:a}{\|a\| \|b\|} \right]^2\|a\|\|b\|}$. Alpha is reset to a maximum of one, and if it is smaller then one, a safety_factor scales the alpha to make sure that the 1-alpha won't get to close to zero.

§ convert_array_to_point()

template<int dim>
Point<dim> aspect::Utilities::convert_array_to_point ( const std::array< double, dim > &  array)

Converts an array of size dim to a Point of size dim.

§ convert_point_to_array()

template<int dim>
std::array<double,dim> aspect::Utilities::convert_point_to_array ( const Point< dim > &  point)

Converts a Point of size dim to an array of size dim.

§ create_model_operator_list()

std::vector<Operator> aspect::Utilities::create_model_operator_list ( const std::vector< std::string > &  operator_names)

Create a vector of operator objects out of a list of strings. Each entry in the list must match one of the allowed operations.

§ nth_basis_for_symmetric_tensors()

template<int dim>
SymmetricTensor<2,dim> aspect::Utilities::nth_basis_for_symmetric_tensors ( const unsigned int  k)

A function that returns a SymmetricTensor, whose entries are zero, except for the k'th component, which is set to one. If k is not on the main diagonal the resulting tensor is symmetrized.