Namespaces
	Coordinates

	MapParsing

	tk

Classes
class	AsciiDataBase

class	AsciiDataBoundary

class	AsciiDataInitial

class	AsciiDataLayered

class	AsciiDataProfile

class	NaturalCoordinate

class	Operator

class	StructuredDataLookup

struct	ThousandSep

class	VectorFunctionFromVelocityFunctionObject

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

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)

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

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

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

template<int dim>
std::vector< Point< dim > >	get_unit_support_points (const SimulatorAccess< dim > &simulator_access)

template<int dim>
bool	point_is_in_triangulation (const Mapping< dim > &mapping, const parallel::distributed::Triangulation< dim > &triangulation, const Point< dim > &point, const MPI_Comm mpi_communicator)

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)

Tensor< 2, 3 >	rotation_matrix_from_axis (const Tensor< 1, 3 > &rotation_axis, const double rotation_angle)

Tensor< 2, 3 >	compute_rotation_matrix_for_slice (const Tensor< 1, 3 > &point_one, const Tensor< 1, 3 > &point_two)

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

bool	fexists (const std::string &filename)

bool	fexists (const std::string &filename, const MPI_Comm comm)

bool	filename_is_url (const std::string &filename)

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

void	collect_and_write_file_content (const std::string &filename, const std::string &file_content, 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, const 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	string_to_bool (const std::string &s)

std::vector< bool >	string_to_bool (const std::vector< std::string > &s)

unsigned int	string_to_unsigned_int (const std::string &s)

std::vector< unsigned int >	string_to_unsigned_int (const std::vector< 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< Operator >	create_model_operator_list (const std::vector< std::string > &operator_names)

const std::string	get_model_operator_options ()

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

template<int dim, typename VectorType >
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)

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="")

template<typename T >
std::vector< std::size_t >	compute_sorting_permutation (const std::vector< T > &vector)

template<typename T >
std::vector< T >	apply_permutation (const std::vector< T > &vector, const std::vector< std::size_t > &permutation_vector)

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)

double	wrap_angle (const double angle)

std::array< double, 3 >	zxz_euler_angles_from_rotation_matrix (const Tensor< 2, 3 > &rotation_matrix)

Tensor< 2, 3 >	zxz_euler_angles_to_rotation_matrix (const double phi1, const double theta, const double phi2)

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
	)

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.

§ parse_map_to_double_array()

DEAL_II_DEPRECATED 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	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`
	)

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. If allow_multiple_values_per_key is set to 'true' it also allows entries of the form "key1: value1|value2|value3, etc", in which case the returned vector will have more entries than the provided list_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. This form only allows a single value per key.
Whether or not a background field is required depends on expects_background_field. Requiring a background field inserts "background" into the list_of_keys as the first entry in the list. Some calling functions (like some material models) require values for background fields, while others may not need background values.
Two special keys are recognized: all –> Assign the associated value to all fields. Only one key is allowed in this case. background –> Assign associated value to the background.

Parameters

[in]	key_value_map	The string representation of the map to be parsed.
[in]	list_of_keys	A list of N 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]	expects_background_field	If true, expect N+1 values and allow setting of the background using the key "background".
[in]	property_name	A 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.
[in]	allow_multiple_values_per_key	If true allow having multiple values for each key. If false only allow a single value per key. In either case each key is only allowed to appear once.
[in,out]	n_values_per_key	A pointer to a vector of unsigned integers. If no pointer is handed over nothing happens. If a pointer to an empty vector is handed over, the vector is resized with as many components as keys (+1 if there is a background field). Each value then stores how many values were found for this key. The sum over all entries is the length of the return value of this function. If a pointer to an existing vector with one or more entries is handed over (e.g. a n_values_per_key vector created by a previous call to this function) then this vector is used as expected structure of the input parameter, and it is checked that `key_value_map` fulfills this structure. This can be used to assert that several input parameters prescribe the same number of values to each key.
[in]	allow_missing_keys	Whether to allow that some keys are not set to any values, i.e. they do not appear at all in the `key_value_map`. This also allows a completely empty map.

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. If multiple values per key are allowed, the vector contains first all values for key 1, then all values for key 2 and so forth. Using the n_values_per_key vector allows the caller to associate entries in the returned vector with specific keys.

Deprecated:: : This function is deprecated in favor of the more general Utilities::MapParsing::parse_map_to_double_array() function. Please use the other function instead.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ get_unit_support_points()

template<int dim>

std::vector<Point<dim> > aspect::Utilities::get_unit_support_points ( const SimulatorAccess< dim > & simulator_access )

This function retrieves the unit support points (in the unit cell) for the current element. The DGP element used when 'set Use locally conservative discretization = true' does not have support points. If these elements are in use, a fictitious support point at the cell center is returned for each shape function that corresponds to the pressure variable, whereas the support points for the velocity are correct. The fictitious points don't matter because we only use this function when interpolating the velocity variable, and ignore the evaluation at the pressure support points.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ point_is_in_triangulation()

template<int dim>

bool aspect::Utilities::point_is_in_triangulation	(	const Mapping< dim > &	mapping,
		const parallel::distributed::Triangulation< dim > &	triangulation,
		const Point< dim > &	point,
		const MPI_Comm	mpi_communicator
	)

Given a point point, find out if any of the MPI processes own the cell in which this point lies. If not, the point lies outside the triangulation.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ rotation_matrix_from_axis()

Tensor<2,3> aspect::Utilities::rotation_matrix_from_axis	(	const Tensor< 1, 3 > &	rotation_axis,
		const double	rotation_angle
	)

A function that returns the corresponding euler angles for a rotation described by rotation axis and angle.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ compute_rotation_matrix_for_slice()

Tensor<2,3> aspect::Utilities::compute_rotation_matrix_for_slice	(	const Tensor< 1, 3 > &	point_one,
		const Tensor< 1, 3 > &	point_two
	)

Compute the 3d rotation matrix that describes the rotation of a plane defined by the two points point_one and point_two onto the x-y-plane in a way that the vector from the origin to point_one points into the (0,1,0) direction after the rotation.

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ 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

Referenced by aspect::Utilities::MapParsing::Options::Options().

§ fexists() [1/2]

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

Checks whether a file named filename exists and is readable.

Note: This function performs file access on all MPI ranks that call it. Only use this function if the file access is performed on all MPI ranks and on a reasonably fast file system.

Returns: True if the file exists and is readable on all MPI ranks that call this function.

Parameters

filename File to check existence

§ fexists() [2/2]

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

Check whether a file named filename exists and is readable on MPI rank 0. Then, broadcast the result to all MPI ranks.

Note that in contrast to the other fexists() function, this function only checks for the existence of the file on a single process. This is useful to avoid overloading the file system and is sufficient if the subsequent file access is also only performed on MPI rank 0.

Returns: True if the file exists and is readable on MPI rank 0.

Parameters

filename	File to check existence
comm	MPI communicator to use.

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

§ filename_is_url()

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

Checks to see if the user is trying to use data from a url.

Parameters

filename File to check

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. The purpose of this function is to reduce parallel file access to a single file access on process 0.

Parameters

[in]	filename	The name of the ascii file to load. If the file name ends in `.gz`, then the function assumes that the file has been compressed using gzip; it then reads and uncompresses the file before distributing it. If the file name is a URL (starting with either `http://`, `https://`, or `file://`), and if ASPECT has been configured with libDAP, then the file is read from that location via libDAP and the returned string is an ASCII data representation of what was obtained this way.
[in]	comm	The MPI communicator in which the content is distributed.

Returns: A string which contains the data in filename.

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

§ collect_and_write_file_content()

void aspect::Utilities::collect_and_write_file_content	(	const std::string &	filename,
		const std::string &	file_content,
		const MPI_Comm	comm
	)

Collect the content of file_content using MPI_Gather to process 0. Then write the content to the file filename on process 0. The purpose of this function is to reduce parallel file access to process 0. Note that this function assumes that the content from all processes fits into the memory of process 0.

Parameters

[in]	filename	The name of the file to write.
[in]	file_content	The content that should be written to file.
[in]	comm	The MPI communicator from which the content is collected.

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

pathname	String that contains the path to create. '/' is used as directory separator.
mode	Permissions (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,
		const bool	silent
	)

Create directory pathname, optionally printing a message.

Parameters

pathname	String that contains path to create. '/' is used as directory separator.
comm	MPI communicator, used to limit creation of directory to processor 0.
silent	Print a nicely formatted message on processor 0 if set to false and the directory does not exist yet and is therefore created in this function.

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.

§ 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.

§ 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.

§ string_to_bool() [1/2]

bool aspect::Utilities::string_to_bool ( const std::string & s )

Given a string s, convert it to a boolean value.

§ string_to_bool() [2/2]

std::vector<bool> aspect::Utilities::string_to_bool ( const std::vector< std::string > & s )

Given a vector of strings s, convert it to a vector of boolean values.

§ string_to_unsigned_int() [1/2]

unsigned int aspect::Utilities::string_to_unsigned_int ( const std::string & s )

Given a string s, convert it to an unsigned int.

§ string_to_unsigned_int() [2/2]

std::vector<unsigned int> aspect::Utilities::string_to_unsigned_int ( const std::vector< std::string > & s )

Given a vector of strings s, convert it to a vector of unsigned int values.

§ 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.

§ 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 rank-4 tensor (i.e., a positive definite operator mapping rank-2 tensors to rank-2 tensors). By definition, the whole operator is symmetric. In the definition above, $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 operator. 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.

One can do some algebra to determine what the optimal factor is. We did this in the Newton paper (Fraters et al., Geophysical Journal International, 2019) where we derived a factor of $\frac{2\eta(\varepsilon(\mathbf u))}{\left[1-\frac{b:a}{\|a\| \|b\|} \right]^2\|a\|\|b\|}$, which we reset to a maximum of one, and if it is smaller then one, a safety_factor scales the value to make sure that 1-alpha won't get to close to zero. However, as later pointed out by Yimin Jin, the computation is wrong, see https://github.com/geodynamics/aspect/issues/5555. Instead, the function now computes the factor as $(2 \eta) / (a:b + b:a)$, again capped at a maximal value of 1, and using a safety factor from below.

In practice, $a$ and $b$ are almost always parallel to each other, and $a:b + b:a = 2a:b$, in which case one can drop the factor of $2$ everywhere in the computations.

§ 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.

§ get_model_operator_options()

const std::string aspect::Utilities::get_model_operator_options ( )

Create a string of model operators for use in declare_parameters

§ 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.

§ project_cellwise()

template<int dim, typename VectorType >

void aspect::Utilities::project_cellwise	(	const Mapping< dim > &	mapping,
		const DoFHandler< dim > &	dof_handler,
		const unsigned int	component_index,
		const Quadrature< dim > &	quadrature
	)

Compute the cellwise projection of component component_index of function to the finite element space described by dof_handler.

Parameters

[in]	mapping	The mapping object to use.
[in]	dof_handler	The DoFHandler the describes the finite element space to project into and that corresponds to `vec_result`.
[in]	component_index	The component index of the `dof_handler` for which the projection is being performed. This component should be described by a DG finite element.
[in]	quadrature	The quadrature formula to use to evaluate `function` on each cell.
[in]	function	The function to project into the finite element space. This function should store the value of the function at the points described by the 2nd argument into the 3rd argument as an std::vector<double> of size quadrature.size().
[out]	vec_result	The output vector where the projected function will be stored in.

§ throw_linear_solver_failure_exception()

void aspect::Utilities::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 = `""`
	)

Throw an exception that a linear solver failed with some helpful information for the user. This function is needed because we have multiple solvers that all require similar treatment and we would like to keep the output consistent. If the optional parameter output_filename is given, the solver history is additionally written to this file.

solver_name A name that identifies the solver and appears in the error message. function_name The name of the function that used the solver (to identify where in the code a solver failed). solver_controls One or more solver controls that describe the history of the solver(s) that failed. The reason the function takes multiple controls is we sometimes use multi-stage solvers, e.g. we try a cheap solver first, and use an expensive solver if the cheap solver fails. exc The exception that was thrown by the solver when it failed, containing additional information about what happened. mpi_communicator The MPI Communicator of the problem. output_filename An optional file name into which (if present) the solver history will be written.

Note: This function never returns normally. It always exits via an exception, either of type ExcMessage (on rank 0 of the parallel computation) or QuietException (on all other ranks).

Referenced by aspect::internal::BlockSchurGMGPreconditioner< StokesMatrixType, ABlockMatrixType, BTBlockOperatorType, SchurComplementMatrixType, ABlockPreconditionerType, SchurComplementPreconditionerType >::vmult().

§ compute_sorting_permutation()

template<typename T >

std::vector<std::size_t> aspect::Utilities::compute_sorting_permutation ( const std::vector< T > & vector )

inline

Create and return a permutation vector which can be used by the apply_permutation() function to put the vector in sorted order.

Parameters

vector vector to sort

§ apply_permutation()

template<typename T >

std::vector<T> aspect::Utilities::apply_permutation	(	const std::vector< T > &	vector,
		const std::vector< std::size_t > &	permutation_vector
	)

inline

Applies a permutation vector to another vector and return the resulting vector.

Parameters

vector	vector to sort
permutation_vector	The permutation vector used to sort the input vector.

Returns: The permuted input vector.

§ rotation_matrices_random_draw_volume_weighting()

std::vector<Tensor<2,3> > aspect::Utilities::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
	)

Get volume weighted rotation matrices, using random draws to convert to a discrete number of orientations, weighted by volume. The input is a vector of volume fractions and a vector of rotation matrices. The vectors need to have the same length.

Parameters

volume_fractions	a vector of doubles representing the volume fraction of each grain
rotation_matrices	a vector of 2nd order 3D tensors representing the rotation matrix of each grain
n_output_matrices	The number of rotation matrices which are output by this function. This can be different from the number of entries in the volume fraction and rotation matrices vectors.
random_number_generator	a reference to a mt19937 random number generator.

§ wrap_angle()

double aspect::Utilities::wrap_angle ( const double angle )

Wraps angle between 0 and 360 degrees.

§ zxz_euler_angles_from_rotation_matrix()

std::array<double,3> aspect::Utilities::zxz_euler_angles_from_rotation_matrix ( const Tensor< 2, 3 > & rotation_matrix )

Compute Z-X-Z Euler angles (https://en.wikipedia.org/wiki/Euler_angles) from rotation matrix. The Z-X-Z indicates the order of axis rotations to generate the Euler angles.

§ zxz_euler_angles_to_rotation_matrix()

Tensor<2,3> aspect::Utilities::zxz_euler_angles_to_rotation_matrix	(	const double	phi1,
		const double	theta,
		const double	phi2
	)

Compute rotation matrix from Z-X-Z Euler angles (https://en.wikipedia.org/wiki/Euler_angles) The Z-X-Z indicates the order of axis axis rotations to generate the Euler angles.

Namespaces

Classes

Functions

Detailed Description

Function Documentation

§ possibly_extend_from_1_to_N()

§ parse_map_to_double_array()

§ parse_input_table()

§ expand_dimensional_variable_names()

§ extract_locally_active_dofs_with_component()

§ get_unit_support_points()

§ point_is_in_triangulation()

§ polygon_contains_point()

§ signed_distance_to_polygon()

§ distance_to_line()

§ orthogonal_vectors()

§ rotation_matrix_from_axis()

§ compute_rotation_matrix_for_slice()

§ real_spherical_harmonic()

§ fexists() [1/2]

§ fexists() [2/2]

§ filename_is_url()

§ read_and_distribute_file_content()

§ collect_and_write_file_content()

§ mkdirp()

§ create_directory()

§ extract_composition_values_at_q_point()

§ expand_ASPECT_SOURCE_DIR()

§ parenthesize_if_nonempty()

§ string_to_bool() [1/2]

§ string_to_bool() [2/2]

§ string_to_unsigned_int() [1/2]

§ string_to_unsigned_int() [2/2]

§ has_unique_entries()

§ weighted_p_norm_average()

§ derivative_of_weighted_p_norm_average()

§ compute_spd_factor()

§ convert_array_to_point()

§ convert_point_to_array()

§ create_model_operator_list()

§ get_model_operator_options()

§ nth_basis_for_symmetric_tensors()

§ project_cellwise()

§ throw_linear_solver_failure_exception()

§ compute_sorting_permutation()

§ apply_permutation()

§ rotation_matrices_random_draw_volume_weighting()

§ wrap_angle()

§ zxz_euler_angles_from_rotation_matrix()

§ zxz_euler_angles_to_rotation_matrix()