diff --git a/include/aspect/particle/property/elastic_tensor_decomposition.h b/include/aspect/particle/property/elastic_tensor_decomposition.h new file mode 100644 index 00000000000..99204148e2f --- /dev/null +++ b/include/aspect/particle/property/elastic_tensor_decomposition.h @@ -0,0 +1,270 @@ +/* + Copyright (C) 2023 by the authors of the ASPECT code. + This file is part of ASPECT. + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + +#ifndef _aspect_particle_property_elastic_tensor_decomposition_h +#define _aspect_particle_property_elastic_tensor_decomposition_h + +#include +#include + +namespace aspect +{ + namespace Particle + { + namespace Property + { + namespace Utilities + { + /** + * Return an even permutation based on an index. This is an internal + * utilities function, also used by the unit tester. + */ + std::array + indexed_even_permutation(const unsigned int index); + + /** + * Computes the Voigt stiffness tensor from the elastic tensor. + * The Voigt stiffness tensor (see Browaeys and Chevrot, 2004) + * defines the stress needed to cause an isotropic strain in the + * material. + */ + SymmetricTensor<2,3> + compute_voigt_stiffness_tensor(const SymmetricTensor<2,6> &elastic_tensor); + + /** + * Computes the dilatation stiffness tensor from the elastic tensor + * The dilatational stiffness tensor (see Browaeys and Chevrot, 2004) + * defines the stress to cause isotropic dilatation in the material. + */ + SymmetricTensor<2,3> + compute_dilatation_stiffness_tensor(const SymmetricTensor<2,6> &elastic_tensor); + + /** + * Computes the Symmetry Cartesian Coordinate System (SCCS). + * + * This is based on Browaeys and Chevrot (2004), GJI (doi: 10.1111/j.1365-246X.2004.024115.x), + * which states at the end of paragraph 3.3 that "... an important property of an orthogonal projection + * is that the distance between a vector $X$ and its orthogonal projection $X_H = p(X)$ on a given + * subspace is minimum. These two features ensure that the decomposition is optimal once a 3-D Cartesian + * coordinate system is chosen.". The other property they talk about is that "The space of elastic + * vectors has a finite dimension [...], i.e. using a different norm from eq. 2.3 will change distances + * but not the resulting decomposition.". + * + * With the three SCCS directions, the elastic tensor can be decomposed into the different + * symmetries in those three SCCS direction, that is, triclinic, monoclinic, orthorhombic, + * tetragonal, hexagonal, and isotropic (Browaeys & Chevrot, 2004). + * + * The dilatation_stiffness_tensor defines the stress to cause isotropic dilatation in the material. + * The voigt_stiffness_tensor defines the stress needed to cause an isotropic strain in the material + */ + Tensor<2,3> compute_unpermutated_SCCS(const SymmetricTensor<2,3> &dilatation_stiffness_tensor, + const SymmetricTensor<2,3> &voigt_stiffness_tensor); + + /** + * Uses the Symmetry Cartesian Coordinate System (SCCS) to try the different permutations to + * determine what is the best projection. + * This is based on Browaeys and Chevrot (2004), GJI (doi: 10.1111/j.1365-246X.2004.024115.x), + * which states at the end of paragraph 3.3 that "... an important property of an orthogonal projection + * is that the distance between a vector $X$ and its orthogonal projection $X_H = p(X)$ on a given + * subspace is minimum. These two features ensure that the decomposition is optimal once a 3-D Cartesian + * coordinate system is chosen.". The other property they talk about is that "The space of elastic + * vectors has a finite dimension [...], i.e. using a different norm from eq. 2.3 will change distances + * but not the resulting decomposition.". + */ + std::array,7> + compute_elastic_tensor_SCCS_decompositions( + const Tensor<2,3> &unpermutated_SCCS, + const SymmetricTensor<2,6> &elastic_matrix); + + + /** + * The tensors below can be used to project matrices to different symmetries. + * See Browaeys and Chevrot, 2004. + */ + static const SymmetricTensor<2,21> projection_matrix_triclinic_to_monoclinic( + Tensor<2,21>( + { + {1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0}, + {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1}, + }) + ); + static const SymmetricTensor<2,9> projection_matrix_monoclinic_to_orthorhombic( + Tensor<2,9>( + { + {1,0,0,0,0,0,0,0,0}, + {0,1,0,0,0,0,0,0,0}, + {0,0,1,0,0,0,0,0,0}, + {0,0,0,1,0,0,0,0,0}, + {0,0,0,0,1,0,0,0,0}, + {0,0,0,0,0,1,0,0,0}, + {0,0,0,0,0,0,1,0,0}, + {0,0,0,0,0,0,0,1,0}, + {0,0,0,0,0,0,0,0,1} + }) + ); + static const SymmetricTensor<2,9> projection_matrix_orthorhombic_to_tetragonal( + Tensor<2,9>( + { + {0.5,0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0}, + {0.5,0.5,0.0,0.0,0.0,0.0,0.0,0.0,0.0}, + {0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0}, + {0.0,0.0,0.0,0.5,0.5,0.0,0.0,0.0,0.0}, + {0.0,0.0,0.0,0.5,0.5,0.0,0.0,0.0,0.0}, + {0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0}, + {0.0,0.0,0.0,0.0,0.0,0.0,0.5,0.5,0.0}, + {0.0,0.0,0.0,0.0,0.0,0.0,0.5,0.5,0.0}, + {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0} + }) + ); + static const SymmetricTensor<2,9> projection_matrix_tetragonal_to_hexagonal( + Tensor<2,9>( + { + {3./8. , 3./8. , 0.0, 0.0, 0.0, 1./(4.*sqrt(2.)) , 0.0, 0.0, 1./4. }, + {3./8. , 3./8. , 0.0, 0.0, 0.0, 1./(4.*sqrt(2.)) , 0.0, 0.0, 1./4. }, + {0.0 , 0.0 , 1.0, 0.0, 0.0, 0.0 , 0.0, 0.0, 0.0 }, + {0.0 , 0.0 , 0.0, 0.5, 0.5, 0.0 , 0.0, 0.0, 0.0 }, + {0.0 , 0.0 , 0.0, 0.5, 0.5, 0.0 , 0.0, 0.0, 0.0 }, + {1./(4.*sqrt(2.)), 1./(4.*sqrt(2.)), 0.0, 0.0, 0.0, 3./4. , 0.0, 0.0, -1./(2.*sqrt(2.))}, + {0.0 , 0.0 , 0.0, 0.0, 0.0, 0.0 , 0.5, 0.5, 0.0 }, + {0.0 , 0.0 , 0.0, 0.0, 0.0, 0.0 , 0.5, 0.5, 0.0 }, + {1./4. , 1./4. , 0.0, 0.0, 0.0, -1./(2.*sqrt(2.)) , 0.0, 0.0, 0.5 } + }) + ); + static const SymmetricTensor<2,9> projection_matrix_hexagonal_to_isotropic( + Tensor<2,9>( + { + {3./15. , 3./15. , 3./15. , sqrt(2.)/15. , sqrt(2.)/15. , sqrt(2.)/15. , 2./15. , 2./15. , 2./15. }, + {3./15. , 3./15. , 3./15. , sqrt(2.)/15. , sqrt(2.)/15. , sqrt(2.)/15. , 2./15. , 2./15. , 2./15. }, + {3./15. , 3./15. , 3./15. , sqrt(2.)/15. , sqrt(2.)/15. , sqrt(2.)/15. , 2./15. , 2./15. , 2./15. }, + {sqrt(2.)/15., sqrt(2.)/15., sqrt(2.)/15., 4./15. , 4./15. , 4./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15. }, + {sqrt(2.)/15., sqrt(2.)/15., sqrt(2.)/15., 4./15. , 4./15. , 4./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15. }, + {sqrt(2.)/15., sqrt(2.)/15., sqrt(2.)/15., 4./15. , 4./15. , 4./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15. }, + {2./15. , 2./15. , 2./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15., 1./5. , 1./5. , 1./5. }, + {2./15. , 2./15. , 2./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15., 1./5. , 1./5. , 1./5. }, + {2./15. , 2./15. , 2./15. , -sqrt(2.)/15., -sqrt(2.)/15., -sqrt(2.)/15., 1./5. , 1./5. , 1./5. } + }) + ); + } + + /** + * Computes several properties of a elastic tensor stored on a particle. + * These include the eigenvectors and conversions between different forms of 4th order tensors. + * + * @ingroup ParticleProperties + */ + template + class ElasticTensorDecomposition : public Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Constructor + */ + ElasticTensorDecomposition(); + + /** + * Initialization function. This function is called once at the + * beginning of the program after parse_parameters is run. + */ + void + initialize () override; + + /** + * Initialization function. This function is called once at the + * creation of every particle for every property to initialize its + * value. + * + * @param [in] position The current particle position. + * @param [in,out] particle_properties The properties of the particle + * that is initialized within the call of this function. The purpose + * of this function should be to extend this vector by a number of + * properties. + */ + void + initialize_one_particle_property (const Point &position, + std::vector &particle_properties) const override; + + /** + * Update function. This function is called every time an update is + * requested by need_update() for every particle for every property. + * + * @param [in] data_position. An unsigned integer that denotes which + * component of the particle property vector is associated with the + * current property. For properties that own several components it + * denotes the first component of this property, all other components + * fill consecutive entries in the @p particle_properties vector. + * + * @param [in] position. The current particle position. + * + * @param [in] solution. The values of the solution variables at the + * current particle position. + * + * @param [in] gradients. The gradients of the solution variables at + * the current particle position. + * + * @param [in,out] particle_properties. The properties of the particle + * that is updated within the call of this function. + */ + void + update_one_particle_property (const unsigned int data_position, + const Point &position, + const Vector &solution, + const std::vector> &gradients, + const ArrayView &particle_properties) const override; + + /** + * This function tells the particle manager that + * we need to update particle properties. + */ + UpdateTimeFlags + need_update () const override; + + /** + * Set up the information about the names and number of components + * this property requires. + * + * @return A vector that contains pairs of the property names and the + * number of components this property plugin defines. + */ + std::vector> + get_property_information() const override; + + private: + unsigned int cpo_elastic_tensor_data_position; + }; + } + } +} + +#endif diff --git a/source/particle/property/elastic_tensor_decomposition.cc b/source/particle/property/elastic_tensor_decomposition.cc new file mode 100644 index 00000000000..3485412d56f --- /dev/null +++ b/source/particle/property/elastic_tensor_decomposition.cc @@ -0,0 +1,472 @@ +/* + Copyright (C) 2023 by the authors of the ASPECT code. + This file is part of ASPECT. + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + +#include +#include +#include +#include + +#include + +namespace aspect +{ + namespace Particle + { + namespace Property + { + namespace Utilities + { + std::array + indexed_even_permutation(const unsigned int index) + { + // there are 6 permutations, but only the odd or even are needed. We use the even + // permutation here. + switch (index) + { + case 0 : + return {{0,1,2}}; + case 1 : + return {{1,2,0}}; + case 2: + return {{2,0,1}}; + /*case 3: + return {0,2,1}; + case 4 : + return {1,0,2}; + case 5: + return {2,1,0};*/ + default: + AssertThrow(false,ExcMessage("Provided index larger then 2 (" + std::to_string(index)+ ").")); + return {{0,0,0}}; + } + + } + + + + SymmetricTensor<2,3> + compute_voigt_stiffness_tensor(const SymmetricTensor<2,6> &elastic_matrix) + { + /** + * the Voigt stiffness tensor (see Browaeys and Chevrot, 2004) + * It defines the stress needed to cause an isotropic strain in the + * material + */ + SymmetricTensor<2,3> voigt_stiffness_tensor; + voigt_stiffness_tensor[0][0]=elastic_matrix[0][0]+elastic_matrix[5][5]+elastic_matrix[4][4]; + voigt_stiffness_tensor[1][1]=elastic_matrix[5][5]+elastic_matrix[1][1]+elastic_matrix[3][3]; + voigt_stiffness_tensor[2][2]=elastic_matrix[4][4]+elastic_matrix[3][3]+elastic_matrix[2][2]; + voigt_stiffness_tensor[1][0]=elastic_matrix[0][5]+elastic_matrix[1][5]+elastic_matrix[3][4]; + voigt_stiffness_tensor[2][0]=elastic_matrix[0][4]+elastic_matrix[2][4]+elastic_matrix[3][5]; + voigt_stiffness_tensor[2][1]=elastic_matrix[1][3]+elastic_matrix[2][3]+elastic_matrix[4][5]; + + return voigt_stiffness_tensor; + } + + + + SymmetricTensor<2,3> + compute_dilatation_stiffness_tensor(const SymmetricTensor<2,6> &elastic_matrix) + { + /** + * The dilatational stiffness tensor (see Browaeys and Chevrot, 2004) + * It defines the stress to cause isotropic dilatation in the material. + */ + SymmetricTensor<2,3> dilatation_stiffness_tensor; + for (size_t i = 0; i < 3; i++) + { + dilatation_stiffness_tensor[0][0]=elastic_matrix[0][i]+dilatation_stiffness_tensor[0][0]; + dilatation_stiffness_tensor[1][1]=elastic_matrix[1][i]+dilatation_stiffness_tensor[1][1]; + dilatation_stiffness_tensor[2][2]=elastic_matrix[2][i]+dilatation_stiffness_tensor[2][2]; + dilatation_stiffness_tensor[1][0]=elastic_matrix[5][i]+dilatation_stiffness_tensor[1][0]; + dilatation_stiffness_tensor[2][0]=elastic_matrix[4][i]+dilatation_stiffness_tensor[2][0]; + dilatation_stiffness_tensor[2][1]=elastic_matrix[3][i]+dilatation_stiffness_tensor[2][1]; + } + return dilatation_stiffness_tensor; + } + + + + Tensor<2,3> + compute_unpermutated_SCCS(const SymmetricTensor<2,3> &dilatation_stiffness_tensor, + const SymmetricTensor<2,3> &voigt_stiffness_tensor) + { + // computing the eigenvector of the dilation and Voigt stiffness matrices and then averaging them by bysection. + const std::array>, 3> voigt_eigenvectors_a = eigenvectors(voigt_stiffness_tensor, SymmetricTensorEigenvectorMethod::jacobi); + const std::array>, 3> dilatation_eigenvectors_a = eigenvectors(dilatation_stiffness_tensor, SymmetricTensorEigenvectorMethod::jacobi); + + + std::array,3> unpermutated_SCCS; + // Averaging dilatation and voigt eigenvectors + // To do this we need to find the eigenvectors which are closest to each other and average those. + // The next function looks for the smallest angle and returns the corresponding vector index for + // that vector. + int vector_index_signed = 0.; + for (unsigned int i1 = 0; i1 < 3; i1++) + { + vector_index_signed = 0; + double smallest_angle = std::numeric_limits::max(); // angle D VeCtor + for (unsigned int i2 = 0; i2 < 3; i2++) + { + double dv_dot_product = dilatation_eigenvectors_a[i1].second*voigt_eigenvectors_a[i2].second; + // limit the dot product between 1 and -1 so we can use the arccos function safely. + if (std::abs(dv_dot_product) >= 1.0) + dv_dot_product = std::copysign(1.0,dv_dot_product); + // Compute the angle between the vectors and account for that vector in the opposite + // direction are the same (0 == 180 degrees). So limit the direction of the vectors between + // 0 and 90 degrees such that it represents the minimum angle between the two lines. + const double angle = dv_dot_product < 0.0 ? std::acos(-1.0)-std::acos(dv_dot_product) : std::acos(dv_dot_product); + // store this if the angle is smaller + if (angle < smallest_angle) + { + vector_index_signed = std::copysign(1.0, dv_dot_product)*i2; + smallest_angle = angle; + } + } + + // Adds to the dilatation eigenvectors to the Voigt eigenvectors with the smallest angle + // Note that the voigt eigenvector is multiplied with a value which can be negative, which means it would be a subtraction. + // Lastly we normalize the unpermutated_SCCS. + unpermutated_SCCS[i1] = 0.5*(dilatation_eigenvectors_a[i1].second + static_cast(vector_index_signed)*voigt_eigenvectors_a[std::abs(vector_index_signed)].second); + unpermutated_SCCS[i1] /= unpermutated_SCCS[i1].norm(); + } + + return Tensor<2,3>( + { + {unpermutated_SCCS[0][0],unpermutated_SCCS[0][1],unpermutated_SCCS[0][2]}, + {unpermutated_SCCS[1][0],unpermutated_SCCS[1][1],unpermutated_SCCS[1][2]}, + {unpermutated_SCCS[2][0],unpermutated_SCCS[2][1],unpermutated_SCCS[2][2]} + }); + } + + + + std::array,7> + compute_elastic_tensor_SCCS_decompositions( + const Tensor<2,3> &unpermutated_SCCS, + const SymmetricTensor<2,6> &elastic_matrix) + { + /** + * Try the different permutations to determine what is the best hexagonal projection. + * This is based on Browaeys and Chevrot (2004), GJI (doi: 10.1111/j.1365-246X.2004.024115.x), + * which states at the end of paragraph 3.3 that "... an important property of an orthogonal projection + * is that the distance between a vector $X$ and its orthogonal projection $X_H = p(X)$ on a given + * subspace is minimum. These two features ensure that the decomposition is optimal once a 3-D Cartesian + * coordinate system is chosen.". The other property they talk about is that "The space of elastic + * vectors has a finite dimension [...], i.e. using a different norm from eq. (2.3 will change distances + * but not the resulting decomposition.". + */ + std::array,3> permutated; + std::array,3> rotated_elastic_matrix; + // The norms variable contains the square norms of the different symmetry approximations of the elastic tensor. + std::array,7> squared_norms_to_projections; + + for (unsigned int permutation_i = 0; permutation_i < 3; permutation_i++) + { + std::array permutation = indexed_even_permutation(permutation_i); + + for (size_t j = 0; j < 3; j++) + { + permutated[permutation_i][j] = unpermutated_SCCS[permutation[j]]; + } + + rotated_elastic_matrix[permutation_i] = aspect::Utilities::Tensors::rotate_voigt_stiffness_matrix((permutated[permutation_i]),elastic_matrix); + + const Tensor<1,21> full_elastic_vector_rotated = aspect::Utilities::Tensors::to_voigt_stiffness_vector(rotated_elastic_matrix[permutation_i]); + + + const double full_norm_square = full_elastic_vector_rotated.norm_square(); + squared_norms_to_projections[6][permutation_i] = full_norm_square; + + // Get the monoclinic and higher symmetry axes, which can be comptued by taking specific elements from the full vector. + // This means that this vector contains all symmetry axes, but the isotropic part is removed. + // The following line would do the same as the lines below, but is is very slow. It has therefore been + // replaced by the lines below. + //auto monoclinic_and_higher_vector = projection_matrix_triclinic_to_monoclinic*full_elastic_vector_rotated; + dealii::Tensor<1,21> monoclinic_and_higher_vector; + monoclinic_and_higher_vector[0] = full_elastic_vector_rotated[0]; + monoclinic_and_higher_vector[1] = full_elastic_vector_rotated[1]; + monoclinic_and_higher_vector[2] = full_elastic_vector_rotated[2]; + monoclinic_and_higher_vector[3] = full_elastic_vector_rotated[3]; + monoclinic_and_higher_vector[4] = full_elastic_vector_rotated[4]; + monoclinic_and_higher_vector[5] = full_elastic_vector_rotated[5]; + monoclinic_and_higher_vector[6] = full_elastic_vector_rotated[6]; + monoclinic_and_higher_vector[7] = full_elastic_vector_rotated[7]; + monoclinic_and_higher_vector[8] = full_elastic_vector_rotated[8]; + monoclinic_and_higher_vector[11] = full_elastic_vector_rotated[11]; + monoclinic_and_higher_vector[14] = full_elastic_vector_rotated[14]; + monoclinic_and_higher_vector[17] = full_elastic_vector_rotated[17]; + monoclinic_and_higher_vector[20] = full_elastic_vector_rotated[20]; + + // The triclinic vector is the full elastic tensor minux the monoclinic and higher symmetry axes vector. + auto triclinic_vector = full_elastic_vector_rotated-monoclinic_and_higher_vector; + squared_norms_to_projections[0][permutation_i] = triclinic_vector.norm_square(); + + // Only the first 9 elements are now non-zero, so crop the vector to the first 9 elements. + // The following line would do the same as the lines below, but it is slow. It has therefore been + // replaced by the lines below. + //auto orthorhombic_and_higher_vector = projection_matrix_monoclinic_to_orthorhombic*monoclinic_and_higher_vector; + dealii::Tensor<1,9> monoclinic_and_higher_vector_cropped; + monoclinic_and_higher_vector_cropped[0] = monoclinic_and_higher_vector[0]; + monoclinic_and_higher_vector_cropped[1] = monoclinic_and_higher_vector[1]; + monoclinic_and_higher_vector_cropped[2] = monoclinic_and_higher_vector[2]; + monoclinic_and_higher_vector_cropped[3] = monoclinic_and_higher_vector[3]; + monoclinic_and_higher_vector_cropped[4] = monoclinic_and_higher_vector[4]; + monoclinic_and_higher_vector_cropped[5] = monoclinic_and_higher_vector[5]; + monoclinic_and_higher_vector_cropped[6] = monoclinic_and_higher_vector[6]; + monoclinic_and_higher_vector_cropped[7] = monoclinic_and_higher_vector[7]; + monoclinic_and_higher_vector_cropped[8] = monoclinic_and_higher_vector[8]; + dealii::Tensor<1,9> orthorhombic_and_higher_vector; + orthorhombic_and_higher_vector[0] = monoclinic_and_higher_vector[0]; + orthorhombic_and_higher_vector[1] = monoclinic_and_higher_vector[1]; + orthorhombic_and_higher_vector[2] = monoclinic_and_higher_vector[2]; + orthorhombic_and_higher_vector[3] = monoclinic_and_higher_vector[3]; + orthorhombic_and_higher_vector[4] = monoclinic_and_higher_vector[4]; + orthorhombic_and_higher_vector[5] = monoclinic_and_higher_vector[5]; + orthorhombic_and_higher_vector[6] = monoclinic_and_higher_vector[6]; + orthorhombic_and_higher_vector[7] = monoclinic_and_higher_vector[7]; + orthorhombic_and_higher_vector[8] = monoclinic_and_higher_vector[8]; + + // The monoclinic vector is the monoclinic and higher symmetry axes vector minux the orthoclinic and higher symmetry axes vector. + auto monoclinic_vector = monoclinic_and_higher_vector_cropped-orthorhombic_and_higher_vector; + squared_norms_to_projections[1][permutation_i] = monoclinic_vector.norm_square(); + + + auto tetragonal_and_higher_vector = projection_matrix_orthorhombic_to_tetragonal*orthorhombic_and_higher_vector; + auto orthorhombic_vector = orthorhombic_and_higher_vector-tetragonal_and_higher_vector; + squared_norms_to_projections[2][permutation_i] = orthorhombic_vector.norm_square(); + + auto hexagonal_and_higher_vector = projection_matrix_tetragonal_to_hexagonal*tetragonal_and_higher_vector; + auto tetragonal_vector = tetragonal_and_higher_vector-hexagonal_and_higher_vector; + squared_norms_to_projections[3][permutation_i] = tetragonal_vector.norm_square(); + + auto isotropic_vector = projection_matrix_hexagonal_to_isotropic*hexagonal_and_higher_vector; + auto hexagonal_vector = hexagonal_and_higher_vector-isotropic_vector; + squared_norms_to_projections[4][permutation_i] = hexagonal_vector.norm_square(); + squared_norms_to_projections[5][permutation_i] = isotropic_vector.norm_square(); + + } + return squared_norms_to_projections; + + } + } + + template + ElasticTensorDecomposition::ElasticTensorDecomposition () + {} + + + + template + void + ElasticTensorDecomposition::initialize () + { + const Particle::Property::Manager &manager = this->get_particle_world().get_property_manager(); + AssertThrow(manager.plugin_name_exists("crystal preferred orientation"), + ExcMessage("No cpo property plugin found.")); + AssertThrow(manager.plugin_name_exists("cpo elastic tensor"), + ExcMessage("No cpo elastic tensor property plugin found.")); + + AssertThrow(manager.check_plugin_order("crystal preferred orientation","elastic tensor decomposition"), + ExcMessage("To use the elastic tensor decomposition plugin, the cpo plugin needs to be defined before this plugin.")); + + AssertThrow(manager.check_plugin_order("cpo elastic tensor","elastic tensor decomposition"), + ExcMessage("To use the elastic tensor decomposition plugin, the cpo elastic tensor plugin needs to be defined before this plugin.")); + + cpo_elastic_tensor_data_position = manager.get_data_info().get_position_by_plugin_index(manager.get_plugin_index_by_name("cpo elastic tensor")); + } + + + + template + void + ElasticTensorDecomposition::initialize_one_particle_property(const Point &, + std::vector &data) const + { + const SymmetricTensor<2,6> elastic_matrix = Particle::Property::CpoElasticTensor::get_elastic_tensor(cpo_elastic_tensor_data_position, + data); + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = Property::Utilities::compute_dilatation_stiffness_tensor(elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = Property::Utilities::compute_voigt_stiffness_tensor(elastic_matrix); + const Tensor<2,3> SCCS_full = Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + + const std::array,7 > norms = Property::Utilities::compute_elastic_tensor_SCCS_decompositions(SCCS_full, elastic_matrix); + + // get max hexagonal element index, which is the same as the permutation index + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + // reorder the SCCS be the SCCS permutation which yields the largest hexagonal vector (percentage of anisotropy) + Tensor<2,3> hexagonal_permutated_SCCS; + for (size_t index = 0; index < 3; ++index) + { + hexagonal_permutated_SCCS[index] = SCCS_full[permutation[index]]; + } + + data.push_back(SCCS_full[0][0]); + data.push_back(SCCS_full[0][1]); + data.push_back(SCCS_full[0][2]); + data.push_back(SCCS_full[1][0]); + data.push_back(SCCS_full[1][1]); + data.push_back(SCCS_full[1][2]); + data.push_back(SCCS_full[2][0]); + data.push_back(SCCS_full[2][1]); + data.push_back(SCCS_full[2][2]); + data.push_back(hexagonal_permutated_SCCS[2][0]); + data.push_back(hexagonal_permutated_SCCS[2][1]); + data.push_back(hexagonal_permutated_SCCS[2][2]); + data.push_back(norms[6][0]); + data.push_back(norms[0][0]); // triclinic + data.push_back(norms[0][1]); // triclinic + data.push_back(norms[0][2]); // triclinic + data.push_back(norms[1][0]); // monoclinic + data.push_back(norms[1][1]); // monoclinic + data.push_back(norms[1][2]); // monoclinic + data.push_back(norms[2][0]); // orthorhomic + data.push_back(norms[2][1]); // orthorhomic + data.push_back(norms[2][2]); // orthorhomic + data.push_back(norms[3][0]); // tetragonal + data.push_back(norms[3][1]); // tetragonal + data.push_back(norms[3][2]); // tetragonal + data.push_back(norms[4][0]); // hexagonal + data.push_back(norms[4][1]); // hexagonal + data.push_back(norms[4][2]); // hexagonal + data.push_back(norms[5][0]); // isotropic + + } + + + + template + void + ElasticTensorDecomposition::update_one_particle_property(const unsigned int data_position, + const Point &, + const Vector &, + const std::vector> &, + const ArrayView &data) const + { + const SymmetricTensor<2,6> elastic_matrix = Particle::Property::CpoElasticTensor::get_elastic_tensor(cpo_elastic_tensor_data_position, + data); + + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = Utilities::compute_dilatation_stiffness_tensor(elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = Utilities::compute_voigt_stiffness_tensor(elastic_matrix); + const Tensor<2,3> SCCS_full = Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + + const std::array,7 > norms = Utilities::compute_elastic_tensor_SCCS_decompositions(SCCS_full, elastic_matrix); + + // get max hexagonal element index, which is the same as the permutation index + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = Utilities::indexed_even_permutation(max_hexagonal_element_index); + // reorder the SCCS by the SCCS permutation which yields the largest hexagonal vector (percentage of anisotropy) + Tensor<2,3> hexagonal_permutated_SCCS; + for (size_t index = 0; index < 3; ++index) + { + hexagonal_permutated_SCCS[index] = SCCS_full[permutation[index]]; + } + + data[data_position] = SCCS_full[0][0]; + data[data_position+1] = SCCS_full[0][1]; + data[data_position+2] = SCCS_full[0][2]; + data[data_position+3] = SCCS_full[1][0]; + data[data_position+4] = SCCS_full[1][1]; + data[data_position+5] = SCCS_full[1][2]; + data[data_position+6] = SCCS_full[2][0]; + data[data_position+7] = SCCS_full[2][1]; + data[data_position+8] = SCCS_full[2][2]; + data[data_position+9] = hexagonal_permutated_SCCS[2][0]; + data[data_position+10] = hexagonal_permutated_SCCS[2][1]; + data[data_position+11] = hexagonal_permutated_SCCS[2][2]; + data[data_position+12] = norms[6][0]; + data[data_position+13] = norms[0][0]; // triclinic + data[data_position+14] = norms[0][1]; // triclinic + data[data_position+15] = norms[0][2]; // triclinic + data[data_position+16] = norms[1][0]; // monoclinic + data[data_position+17] = norms[1][1]; // monoclinic + data[data_position+18] = norms[1][2]; // monoclinic + data[data_position+19] = norms[2][0]; // orthorhomic + data[data_position+20] = norms[2][1]; // orthorhomic + data[data_position+21] = norms[2][2]; // orthorhomic + data[data_position+22] = norms[3][0]; // tetragonal + data[data_position+23] = norms[3][1]; // tetragonal + data[data_position+24] = norms[3][2]; // tetragonal + data[data_position+25] = norms[4][0]; // hexagonal + data[data_position+26] = norms[4][1]; // hexagonal + data[data_position+27] = norms[4][2]; // hexagonal + data[data_position+28] = norms[5][0]; // isotropic + } + + + + template + UpdateTimeFlags + ElasticTensorDecomposition::need_update() const + { + return update_output_step; + } + + + + template + std::vector> + ElasticTensorDecomposition::get_property_information() const + { + std::vector> property_information; + + property_information.push_back(std::make_pair("cpo elastic axis e1",3)); + property_information.push_back(std::make_pair("cpo elastic axis e2",3)); + property_information.push_back(std::make_pair("cpo elastic axis e3",3)); + property_information.push_back(std::make_pair("cpo elastic hexagonal axis",3)); + property_information.push_back(std::make_pair("cpo elastic vector norm square",1)); + property_information.push_back(std::make_pair("cpo elastic triclinic vector norm square p1",1)); + property_information.push_back(std::make_pair("cpo elastic triclinic vector norm square p2",1)); + property_information.push_back(std::make_pair("cpo elastic triclinic vector norm square p3",1)); + property_information.push_back(std::make_pair("cpo elastic monoclinic vector norm square p1",1)); + property_information.push_back(std::make_pair("cpo elastic monoclinic vector norm square p2",1)); + property_information.push_back(std::make_pair("cpo elastic monoclinic vector norm square p3",1)); + property_information.push_back(std::make_pair("cpo elastic orthorhombic vector norm square p1",1)); + property_information.push_back(std::make_pair("cpo elastic orthorhombic vector norm square p2",1)); + property_information.push_back(std::make_pair("cpo elastic orthorhombic vector norm square p3",1)); + property_information.push_back(std::make_pair("cpo elastic tetragonal vector norm square p1",1)); + property_information.push_back(std::make_pair("cpo elastic tetragonal vector norm square p2",1)); + property_information.push_back(std::make_pair("cpo elastic tetragonal vector norm square p3",1)); + property_information.push_back(std::make_pair("cpo elastic hexagonal vector norm square p1",1)); + property_information.push_back(std::make_pair("cpo elastic hexagonal vector norm square p2",1)); + property_information.push_back(std::make_pair("cpo elastic hexagonal vector norm square p3",1)); + property_information.push_back(std::make_pair("cpo elastic isotropic vector norm square",1)); + + return property_information; + } + } + } +} + + + +// explicit instantiations +namespace aspect +{ + namespace Particle + { + namespace Property + { + ASPECT_REGISTER_PARTICLE_PROPERTY(ElasticTensorDecomposition, + "elastic tensor decomposition", + "A plugin which decomposes the elastic tensor into different approximations " + "(Isotropic, Hexagonal, Tetragonal, Orthorhombic, Monoclinic and Triclinic) " + "and provides the eigenvectors of the tensor.") + } + } +} diff --git a/unit_tests/crystal_preferred_orientation.cc b/unit_tests/crystal_preferred_orientation.cc index b47d88d7254..72e8f5c6b94 100644 --- a/unit_tests/crystal_preferred_orientation.cc +++ b/unit_tests/crystal_preferred_orientation.cc @@ -19,9 +19,11 @@ */ #include "common.h" +#include #include #include #include +#include /** @@ -1447,3 +1449,856 @@ TEST_CASE("CPO elastic tensor") CHECK(array_plus_100[cpo_data_position + i] == tensor[dealii::SymmetricTensor<2,6>::unrolled_to_component_indices(i)]); } } + + + + +TEST_CASE("LPO elastic tensor decomposition") +{ + using namespace dealii; + using namespace aspect::Particle::Property; + using namespace aspect::Particle::Property::Utilities; + { + auto elastic_composition_class = ElasticTensorDecomposition<3>(); + + // the matrix from Browaeys and Chevrot, 2004, Geophys. J. Int. (doi: 10.1111/j.1365-246X.2004.02415.x) + // Note that the computed norms are not in the papers, but are an extra check. + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 192; + full_elastic_matrix[0][1] = 66; + full_elastic_matrix[0][2] = 60; + full_elastic_matrix[1][1] = 160; + full_elastic_matrix[1][2] = 56; + full_elastic_matrix[2][2] = 272; + full_elastic_matrix[3][3] = 60; + full_elastic_matrix[4][4] = 62; + full_elastic_matrix[5][5] = 49; + + SymmetricTensor<2,6> reference_isotropic_matrix; + reference_isotropic_matrix[0][0] = 194.7; + reference_isotropic_matrix[0][1] = 67.3; + reference_isotropic_matrix[0][2] = 67.3; + reference_isotropic_matrix[1][1] = 194.7; + reference_isotropic_matrix[1][2] = 67.3; + reference_isotropic_matrix[2][2] = 194.7; + reference_isotropic_matrix[3][3] = 63.7; + reference_isotropic_matrix[4][4] = 63.7; + reference_isotropic_matrix[5][5] = 63.7; + + SymmetricTensor<2,6> reference_hex_matrix; + reference_hex_matrix[0][0] = -21.7; + reference_hex_matrix[1][1] = -21.7; + reference_hex_matrix[2][2] = 77.3; + reference_hex_matrix[1][2] = -9.3; + reference_hex_matrix[0][2] = -9.3; + reference_hex_matrix[0][1] = 1.7; + reference_hex_matrix[3][3] = -2.7; + reference_hex_matrix[4][4] = -2.7; + reference_hex_matrix[5][5] = -11.7; + + SymmetricTensor<2,6> reference_T_matrix; + reference_T_matrix[0][0] = 3; + reference_T_matrix[1][1] = 3; + reference_T_matrix[0][1] = -3; + reference_T_matrix[5][5] = -3; + + SymmetricTensor<2,6> reference_O_matrix; + reference_O_matrix[0][0] = 16; + reference_O_matrix[1][1] = -16; + reference_O_matrix[1][2] = -2; + reference_O_matrix[0][2] = 2; + reference_O_matrix[3][3] = -1; + reference_O_matrix[4][4] = 1; + + const Tensor<1,21> full_elastic_vector = aspect::Utilities::Tensors::to_voigt_stiffness_vector(full_elastic_matrix); + const Tensor<1,21> reference_isotropic_vector = aspect::Utilities::Tensors::to_voigt_stiffness_vector(reference_isotropic_matrix); + const Tensor<1,21> reference_hex_vector = aspect::Utilities::Tensors::to_voigt_stiffness_vector(reference_hex_matrix); + const Tensor<1,21> reference_T_vector = aspect::Utilities::Tensors::to_voigt_stiffness_vector(reference_T_matrix); + const Tensor<1,21> reference_O_vector = aspect::Utilities::Tensors::to_voigt_stiffness_vector(reference_O_matrix); + + // this matrix is already in the SCCSS, so no rotation needed, only the projection. + CHECK(full_elastic_vector.norm_square() == Approx(198012.0)); + const Tensor<1,21> mono_and_higher_vector = aspect::Particle::Property::Utilities::projection_matrix_triclinic_to_monoclinic * full_elastic_vector; + const Tensor<1,21> tric_vector = full_elastic_vector-mono_and_higher_vector; + for (size_t iii = 0; iii < 21; iii++) + { + CHECK(std::abs(tric_vector[iii]) < 1e-15); + } + CHECK(tric_vector.norm_square() == Approx(0.0)); + + dealii::Tensor<1,9> mono_and_higher_vector_croped; + mono_and_higher_vector_croped[0] = mono_and_higher_vector[0]; + mono_and_higher_vector_croped[1] = mono_and_higher_vector[1]; + mono_and_higher_vector_croped[2] = mono_and_higher_vector[2]; + mono_and_higher_vector_croped[3] = mono_and_higher_vector[3]; + mono_and_higher_vector_croped[4] = mono_and_higher_vector[4]; + mono_and_higher_vector_croped[5] = mono_and_higher_vector[5]; + mono_and_higher_vector_croped[6] = mono_and_higher_vector[6]; + mono_and_higher_vector_croped[7] = mono_and_higher_vector[7]; + mono_and_higher_vector_croped[8] = mono_and_higher_vector[8]; + const Tensor<1,9> ortho_and_higher_vector = aspect::Particle::Property::Utilities::projection_matrix_monoclinic_to_orthorhombic*mono_and_higher_vector_croped; + const Tensor<1,9> mono_vector = mono_and_higher_vector_croped-ortho_and_higher_vector; + for (size_t iii = 0; iii < 9; iii++) + { + CHECK(std::abs(mono_vector[iii]) < 1e-15); + } + CHECK(mono_vector.norm_square() == Approx(0.0)); + + const Tensor<1,9> tetra_and_higher_vector = aspect::Particle::Property::Utilities::projection_matrix_orthorhombic_to_tetragonal*ortho_and_higher_vector; + const Tensor<1,9> ortho_vector = ortho_and_higher_vector-tetra_and_higher_vector; + for (size_t iii = 0; iii < 9; iii++) + { + CHECK(ortho_vector[iii] == Approx(reference_O_vector[iii])); + } + CHECK(ortho_vector.norm_square() == Approx(536.0)); + + const Tensor<1,9> hexa_and_higher_vector = aspect::Particle::Property::Utilities::projection_matrix_tetragonal_to_hexagonal*tetra_and_higher_vector; + const Tensor<1,9> tetra_vector = tetra_and_higher_vector-hexa_and_higher_vector; + for (size_t iii = 0; iii < 9; iii++) + { + CHECK(tetra_vector[iii] == Approx(reference_T_vector[iii])); + } + CHECK(tetra_vector.norm_square() == Approx(72.0)); + + const Tensor<1,9> iso_vector = aspect::Particle::Property::Utilities::projection_matrix_hexagonal_to_isotropic*hexa_and_higher_vector; + for (size_t iii = 0; iii < 9; iii++) + { + // something weird is going on with the approx, somehow without a very large + // epsilon it thinks for example that "194.6666666667 == Approx( 194.7 )" is false. + CHECK(iso_vector[iii] == Approx(reference_isotropic_vector[iii]).epsilon(1e-3)); + } + CHECK(iso_vector.norm_square() == Approx(189529.3333333334)); + + const Tensor<1,9> hexa_vector = hexa_and_higher_vector-iso_vector; + for (size_t iii = 0; iii < 9; iii++) + { + // same as above. + CHECK(hexa_vector[iii] == Approx(reference_hex_vector[iii]).epsilon(2.5e-2)); + } + CHECK(hexa_vector.norm_square() == Approx(7874.6666666667)); + + } + + { + // now compute the distribution in a elastic tensor which is not in SCCS. + // We will use the result (with the results of the steps in between) from + // Cowin and Mehrabadi, 1985, Mech. appl. Math. + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 1769.50; + full_elastic_matrix[0][1] = 873.50; + full_elastic_matrix[0][2] = 838.22; + full_elastic_matrix[0][3] = -17.68; + full_elastic_matrix[0][4] = -110.32; + full_elastic_matrix[0][5] = 144.92; + full_elastic_matrix[1][0] = 873.50; + full_elastic_matrix[1][1] = 1846.64; + full_elastic_matrix[1][2] = 836.66; + full_elastic_matrix[1][3] = -37.60; + full_elastic_matrix[1][4] = -32.32; + full_elastic_matrix[1][5] = 153.80; + full_elastic_matrix[2][0] = 838.22; + full_elastic_matrix[2][1] = 836.66; + full_elastic_matrix[2][2] = 1603.28; + full_elastic_matrix[2][3] = -29.68; + full_elastic_matrix[2][4] = -93.52; + full_elastic_matrix[2][5] = 22.40; + full_elastic_matrix[3][0] = -17.68; + full_elastic_matrix[3][1] = -37.60; + full_elastic_matrix[3][2] = -29.68; + full_elastic_matrix[3][3] = 438.77; + full_elastic_matrix[3][4] = 57.68; + full_elastic_matrix[3][5] = -50.50; + full_elastic_matrix[4][0] = -110.32; + full_elastic_matrix[4][1] = -32.32; + full_elastic_matrix[4][2] = -93.52; + full_elastic_matrix[4][3] = 57.68; + full_elastic_matrix[4][4] = 439.79; + full_elastic_matrix[4][5] = -34.78; + full_elastic_matrix[5][0] = 144.92; + full_elastic_matrix[5][1] = 153.80; + full_elastic_matrix[5][2] = 22.40; + full_elastic_matrix[5][3] = -50.50; + full_elastic_matrix[5][4] = -34.78; + full_elastic_matrix[5][5] = 501.80; + full_elastic_matrix /= 81.; + + SymmetricTensor<2,3> reference_dilatation_matrix; + reference_dilatation_matrix[0][0] = 386.80; + reference_dilatation_matrix[0][1] = 35.68; + reference_dilatation_matrix[0][2] = -26.24; + //reference_dilatation_matrix[1][0] = ; + reference_dilatation_matrix[1][1] = 395.20; + reference_dilatation_matrix[1][2] = -9.44; + //reference_dilatation_matrix[2][0] = ; + //reference_dilatation_matrix[2][1] = ; + reference_dilatation_matrix[2][2] = 364.24; + reference_dilatation_matrix /= 9.; + + SymmetricTensor<2,3> reference_voigt_matrix; + reference_voigt_matrix[0][0] = 33.47; + reference_voigt_matrix[0][1] = 4.4; + reference_voigt_matrix[0][2] = -3.14; + //reference_voigt_matrix[1][0] = ; + reference_voigt_matrix[1][1] = 34.41; + reference_voigt_matrix[1][2] = -1.26; + //reference_voigt_matrix[2][0] = ; + //reference_voigt_matrix[2][1] = ; + reference_voigt_matrix[2][2] = 30.64; + + // Q in equation 2.17, but reordered from + // largest to smallest eigenvalue. + Tensor<2,3> reference_SCCS; + reference_SCCS[2][0] = 2./3.; + reference_SCCS[2][1] = -1./3.; + reference_SCCS[2][2] = 2./3.; + reference_SCCS[1][0] = 1./3.; + reference_SCCS[1][1] = -2./3.; + reference_SCCS[1][2] = -2./3.; + reference_SCCS[0][0] = 2./3.; + reference_SCCS[0][1] = 2./3.; + reference_SCCS[0][2] = -1./3.; + + SymmetricTensor<2,6> reference_rotated_elastic_matrix; + reference_rotated_elastic_matrix[0][0] = 18; + reference_rotated_elastic_matrix[0][1] = 9.98; + reference_rotated_elastic_matrix[0][2] = 10.1; + //reference_rotated_elastic_matrix[1][0] = ; + reference_rotated_elastic_matrix[1][1] = 20.2; + reference_rotated_elastic_matrix[1][2] = 10.07; + //reference_rotated_elastic_matrix[2][0] = ; + //reference_rotated_elastic_matrix[2][1] = ; + reference_rotated_elastic_matrix[2][2] = 27.6; + reference_rotated_elastic_matrix[3][3] = 6.23; + reference_rotated_elastic_matrix[4][4] = 5.61; + reference_rotated_elastic_matrix[5][5] = 4.52; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + for (size_t iii = 0; iii < 3; iii++) + { + for (size_t jjj = 0; jjj < 3; jjj++) + { + CHECK(reference_dilatation_matrix[iii][jjj] == Approx(dilatation_stiffness_tensor_full[iii][jjj])); + } + } + + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + for (size_t iii = 0; iii < 3; iii++) + { + for (size_t jjj = 0; jjj < 3; jjj++) + { + CHECK(reference_voigt_matrix[iii][jjj] == Approx(voigt_stiffness_tensor_full[iii][jjj])); + } + } + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + // chose the permutation which fits the solution given in the paper. + // This is not a even permutation. + std::array permutation = {{2,1,0}}; + Tensor<2,3> permutated_SCCS; + for (size_t j = 0; j < 3; j++) + { + permutated_SCCS[j] = unpermutated_SCCS[permutation[j]]; + } + SymmetricTensor<2,6> rotated_elastic_matrix = aspect::Utilities::Tensors::rotate_voigt_stiffness_matrix(permutated_SCCS,full_elastic_matrix); + + for (size_t iii = 0; iii < 6; iii++) + { + for (size_t jjj = 0; jjj < 6; jjj++) + { + CHECK(reference_rotated_elastic_matrix[iii][jjj] == Approx(rotated_elastic_matrix[iii][jjj]).margin(1e-4).epsilon(7.5e-2)); + } + } + } + { + // This is testing some other elastic matrices of which the properties can + // be resonaly predicted (for example, mostly hexagonal) and some coverage + // testing: Case 1. + + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 192; + full_elastic_matrix[0][1] = 66; + full_elastic_matrix[0][2] = 56; + full_elastic_matrix[1][1] = 192; + full_elastic_matrix[1][2] = 56; + full_elastic_matrix[2][2] = 272; + full_elastic_matrix[3][3] = 60; + full_elastic_matrix[4][4] = 60; + full_elastic_matrix[5][5] = 0.5*(full_elastic_matrix[0][0] + full_elastic_matrix[0][1]); + + std::array,7 > reference_norms; + reference_norms[0] = {{0,0,0}}; // no triclinic + reference_norms[1] = {{0,0,0}}; // no monoclinic + reference_norms[2] = {{12822.0,0,12822.0}}; // orthorhomic dependent on direction (remember, this is basically just one grain) + reference_norms[3] = {{1568.0,8712.0,1568.0}}; // tetragonal could be present + reference_norms[4] = {{2759.3333333333,8437.3333333333,2759.3333333333}}; // hexagonal is expected + reference_norms[5] = {{247182.6666666667,247182.6666666667,247182.6666666667}}; // isotropic should be the same for every permutation + reference_norms[6] = {{264332.0,264332.0,264332.0}}; // total should be the same for every direction + + // It is already in a SCCS frame, so they should be + // unit vectors. The order is irrelevant. + Tensor<2,3> reference_unpermutated_SCCS; + reference_unpermutated_SCCS[0][0] = 0; + reference_unpermutated_SCCS[0][1] = 0; + reference_unpermutated_SCCS[0][2] = 1; + reference_unpermutated_SCCS[1][0] = 1; + reference_unpermutated_SCCS[1][1] = 0; + reference_unpermutated_SCCS[1][2] = 0; + reference_unpermutated_SCCS[2][0] = 0; + reference_unpermutated_SCCS[2][1] = 1; + reference_unpermutated_SCCS[2][2] = 0; + + // Permutation 2 has the highest norm, so reference_unpermutated_SCCS + // is permutated by {1,2,0}; + Tensor<2,3> reference_hexa_permutated_SCCS; + reference_hexa_permutated_SCCS[0][0] = 1; + reference_hexa_permutated_SCCS[0][1] = 0; + reference_hexa_permutated_SCCS[0][2] = 0; + reference_hexa_permutated_SCCS[1][0] = 0; + reference_hexa_permutated_SCCS[1][1] = 1; + reference_hexa_permutated_SCCS[1][2] = 0; + reference_hexa_permutated_SCCS[2][0] = 0; + reference_hexa_permutated_SCCS[2][1] = 0; + reference_hexa_permutated_SCCS[2][2] = 1; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("unpermutated_SCCS[" << iii << "] = " << unpermutated_SCCS[iii][0] << ":" << unpermutated_SCCS[iii][1] << ":" << unpermutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_unpermutated_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_unpermutated_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + std::array,7 > norms = aspect::Particle::Property::Utilities::compute_elastic_tensor_SCCS_decompositions(unpermutated_SCCS, full_elastic_matrix); + std::array totals = {{0,0,0}}; + for (size_t iii = 0; iii < 7; iii++) + { + // also check that the total is equal to the full norm. + + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("i:j = " << iii << ":" << jjj ); + CHECK(reference_norms[iii][jjj] == Approx(norms[iii][jjj])); + if (iii < 6) + { + totals[jjj] += norms[iii][jjj]; + } + } + } + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("j = " << jjj ); + CHECK(totals[jjj] == Approx(norms[6][jjj])); + } + + // also check that all the isotropic and totals are the same + for (size_t iii = 1; iii < 3; iii++) + { + CHECK(norms[5][0] == Approx(norms[5][iii])); + CHECK(norms[6][0] == Approx(norms[6][iii])); + } + + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = aspect::Particle::Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + Tensor<2,3> hexa_permutated_SCCS; + for (size_t index = 0; index < 3; index++) + { + hexa_permutated_SCCS[index] = unpermutated_SCCS[permutation[index]]; + } + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("hexa_permutated_SCCS[" << iii << "] = " << hexa_permutated_SCCS[iii][0] << ":" << hexa_permutated_SCCS[iii][1] << ":" << hexa_permutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_hexa_permutated_SCCS[iii][0] == Approx(hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(hexa_permutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_hexa_permutated_SCCS[iii][0] == Approx(-hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(-hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(-hexa_permutated_SCCS[iii][2]).epsilon(epsilon)))); + } + } + { + // This is testing some other elastic matrices of which the properties can + // be resonaly predicted (for example, mostly hexagonal) and some coverage + // testing: Case 2: a more complicated matrix. + + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 225; + full_elastic_matrix[0][1] = 54; + full_elastic_matrix[0][2] = 72; + full_elastic_matrix[1][1] = 214; + full_elastic_matrix[1][2] = 53; + full_elastic_matrix[2][2] = 178; + full_elastic_matrix[3][3] = 78; + full_elastic_matrix[4][4] = 82; + full_elastic_matrix[5][5] = 76; + + std::array,7 > reference_norms; + reference_norms[0] = {{0,0,0}}; // no triclinic + reference_norms[1] = {{0,0,0}}; // no monoclinic + reference_norms[2] = {{453.5,1044.0,1113.5}}; // orthorhomic dependent on direction + reference_norms[3] = {{91.125,84.5,595.125}}; // tetragonal could be present + reference_norms[4] = {{1350.175,766.3,186.175}}; // hexagonal is expected + reference_norms[5] = {{222364.2,222364.2,222364.2}}; // isotropic should be the same for every permutation + reference_norms[6] = {{224259.0,224259.0,224259.0}}; // total should be the same for every direction + + // It is already in a SCCS frame, so they should be + // unit vectors. The order is irrelevant. + Tensor<2,3> reference_unpermutated_SCCS; + reference_unpermutated_SCCS[0][0] = 1; + reference_unpermutated_SCCS[0][1] = 0; + reference_unpermutated_SCCS[0][2] = 0; + reference_unpermutated_SCCS[1][0] = 0; + reference_unpermutated_SCCS[1][1] = 1; + reference_unpermutated_SCCS[1][2] = 0; + reference_unpermutated_SCCS[2][0] = 0; + reference_unpermutated_SCCS[2][1] = 0; + reference_unpermutated_SCCS[2][2] = 1; + + // Permutation 1 has the highest norm, so reference_unpermutated_SCCS + // is permutated by {0,1,2}; + Tensor<2,3> reference_hexa_permutated_SCCS; + reference_hexa_permutated_SCCS[0][0] = 1; + reference_hexa_permutated_SCCS[0][1] = 0; + reference_hexa_permutated_SCCS[0][2] = 0; + reference_hexa_permutated_SCCS[1][0] = 0; + reference_hexa_permutated_SCCS[1][1] = 1; + reference_hexa_permutated_SCCS[1][2] = 0; + reference_hexa_permutated_SCCS[2][0] = 0; + reference_hexa_permutated_SCCS[2][1] = 0; + reference_hexa_permutated_SCCS[2][2] = 1; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("unpermutated_SCCS[" << iii << "] = " << unpermutated_SCCS[iii][0] << ":" << unpermutated_SCCS[iii][1] << ":" << unpermutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_unpermutated_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_unpermutated_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + std::array,7 > norms = aspect::Particle::Property::Utilities::compute_elastic_tensor_SCCS_decompositions(unpermutated_SCCS, full_elastic_matrix); + std::array totals = {{0,0,0}}; + for (size_t iii = 0; iii < 7; iii++) + { + // also check that the total is equal to the full norm. + + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("i:j = " << iii << ":" << jjj ); + CHECK(reference_norms[iii][jjj] == Approx(norms[iii][jjj])); + if (iii < 6) + { + totals[jjj] += norms[iii][jjj]; + } + } + } + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("j = " << jjj ); + CHECK(totals[jjj] == Approx(norms[6][jjj])); + } + + // also check that all the isotropic and totals are the same + for (size_t iii = 1; iii < 3; iii++) + { + CHECK(norms[5][0] == Approx(norms[5][iii])); + CHECK(norms[6][0] == Approx(norms[6][iii])); + } + + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = aspect::Particle::Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + Tensor<2,3> hexa_permutated_SCCS; + for (size_t index = 0; index < 3; index++) + { + hexa_permutated_SCCS[index] = unpermutated_SCCS[permutation[index]]; + } + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("hexa_permutated_SCCS[" << iii << "] = " << hexa_permutated_SCCS[iii][0] << ":" << hexa_permutated_SCCS[iii][1] << ":" << hexa_permutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_hexa_permutated_SCCS[iii][0] == Approx(hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(hexa_permutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_hexa_permutated_SCCS[iii][0] == Approx(-hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(-hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(-hexa_permutated_SCCS[iii][2]).epsilon(epsilon)))); + } + } + + { + // This is testing some other elastic matrices of which the properties can + // be resonaly predicted (for example, mostly hexagonal) and some coverage + // testing: Case 3: olivine from experiments as used in D-Rex. + + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 320.71; + full_elastic_matrix[0][1] = 69.84; + full_elastic_matrix[0][2] = 71.22; + full_elastic_matrix[1][1] = 197.25; + full_elastic_matrix[1][2] = 74.8; + full_elastic_matrix[2][2] = 234.32; + full_elastic_matrix[3][3] = 63.77; + full_elastic_matrix[4][4] = 77.67; + full_elastic_matrix[5][5] = 78.36; + + std::array,7 > reference_norms; + reference_norms[0] = {{0,0,0}}; // no triclinic + reference_norms[1] = {{0,0,0}}; // no monoclinic + reference_norms[2] = {{4181.95385,689.94905,8020.4222}}; // orthorhomic dependent on direction + reference_norms[3] = {{1298.2060125,90.3840125,525.5282}}; // tetragonal could be present + reference_norms[4] = {{4364.6051908333,9064.4319908333,1298.8146533333}}; // hexagonal is expected + reference_norms[5] = {{282871.5975466666,282871.5975466666,282871.5975466666}}; // isotropic should be the same for every permutation + reference_norms[6] = {{292716.3626,292716.3626,292716.3626}}; // total should be the same for every direction + + // It is already in a SCCS frame, so they should be + // unit vectors. The order is irrelevant. + Tensor<2,3> reference_unpermutated_SCCS; + reference_unpermutated_SCCS[0][0] = 1; + reference_unpermutated_SCCS[0][1] = 0; + reference_unpermutated_SCCS[0][2] = 0; + reference_unpermutated_SCCS[1][0] = 0; + reference_unpermutated_SCCS[1][1] = 0; + reference_unpermutated_SCCS[1][2] = 1; + reference_unpermutated_SCCS[2][0] = 0; + reference_unpermutated_SCCS[2][1] = 1; + reference_unpermutated_SCCS[2][2] = 0; + + // Permutation 2 has the highest norm, so reference_unpermutated_SCCS + // is permutated by {1,2,0}; + Tensor<2,3> reference_hexa_permutated_SCCS; + reference_hexa_permutated_SCCS[0][0] = 0; + reference_hexa_permutated_SCCS[0][1] = 0; + reference_hexa_permutated_SCCS[0][2] = 1; + reference_hexa_permutated_SCCS[1][0] = 0; + reference_hexa_permutated_SCCS[1][1] = 1; + reference_hexa_permutated_SCCS[1][2] = 0; + reference_hexa_permutated_SCCS[2][0] = 1; + reference_hexa_permutated_SCCS[2][1] = 0; + reference_hexa_permutated_SCCS[2][2] = 0; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("unpermutated_SCCS[" << iii << "] = " << unpermutated_SCCS[iii][0] << ":" << unpermutated_SCCS[iii][1] << ":" << unpermutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_unpermutated_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_unpermutated_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + std::array,7 > norms = aspect::Particle::Property::Utilities::compute_elastic_tensor_SCCS_decompositions(unpermutated_SCCS, full_elastic_matrix); + std::array totals = {{0,0,0}}; + for (size_t iii = 0; iii < 7; iii++) + { + // also check that the total is equal to the full norm. + + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("i:j = " << iii << ":" << jjj ); + CHECK(reference_norms[iii][jjj] == Approx(norms[iii][jjj])); + if (iii < 6) + { + totals[jjj] += norms[iii][jjj]; + } + } + } + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("j = " << jjj ); + CHECK(totals[jjj] == Approx(norms[6][jjj])); + } + + // also check that all the isotropic and totals are the same + for (size_t iii = 1; iii < 3; iii++) + { + CHECK(norms[5][0] == Approx(norms[5][iii])); + CHECK(norms[6][0] == Approx(norms[6][iii])); + } + + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = aspect::Particle::Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + Tensor<2,3> hexa_permutated_SCCS; + for (size_t index = 0; index < 3; index++) + { + hexa_permutated_SCCS[index] = unpermutated_SCCS[permutation[index]]; + } + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("hexa_permutated_SCCS[" << iii << "] = " << hexa_permutated_SCCS[iii][0] << ":" << hexa_permutated_SCCS[iii][1] << ":" << hexa_permutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_hexa_permutated_SCCS[iii][0] == Approx(hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(hexa_permutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_hexa_permutated_SCCS[iii][0] == Approx(-hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(-hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(-hexa_permutated_SCCS[iii][2]).epsilon(epsilon)))); + } + } + + { + // This is testing some other elastic matrices of which the properties can + // be resonaly predicted (for example, mostly hexagonal) and some coverage + // testing: Case 4: enstatite from experiments as used in D-Rex. + + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 236.9; + full_elastic_matrix[0][1] = 79.6; + full_elastic_matrix[0][2] = 63.2; + full_elastic_matrix[1][1] = 180.5; + full_elastic_matrix[1][2] = 56.8; + full_elastic_matrix[2][2] = 230.4; + full_elastic_matrix[3][3] = 84.3; + full_elastic_matrix[4][4] = 79.4; + full_elastic_matrix[5][5] = 80.1; + + std::array,7 > reference_norms; + reference_norms[0] = {{0,0,0}}; // no triclinic + reference_norms[1] = {{0,0,0}}; // no monoclinic + reference_norms[2] = {{576.245,1514.945,1679.46}}; // orthorhomic dependent on direction + reference_norms[3] = {{67.86125,199.00125,483.605}}; // tetragonal could be present + reference_norms[4] = {{2076.64175,1006.80175,557.683}}; // hexagonal is expected + reference_norms[5] = {{245485.992,245485.992,245485.992}}; // isotropic should be the same for every permutation + reference_norms[6] = {{248206.74,248206.74,248206.74}}; // total should be the same for every direction + + // It is already in a SCCS frame, so they should be + // unit vectors. The order is irrelevant. + Tensor<2,3> reference_unpermutated_SCCS; + reference_unpermutated_SCCS[0][0] = 1; + reference_unpermutated_SCCS[0][1] = 0; + reference_unpermutated_SCCS[0][2] = 0; + reference_unpermutated_SCCS[1][0] = 0; + reference_unpermutated_SCCS[1][1] = 0; + reference_unpermutated_SCCS[1][2] = 1; + reference_unpermutated_SCCS[2][0] = 0; + reference_unpermutated_SCCS[2][1] = 1; + reference_unpermutated_SCCS[2][2] = 0; + + // Permutation 1 has the highest norm, so reference_unpermutated_SCCS + // is permutated by {0,1,2}; + Tensor<2,3> reference_hexa_permutated_SCCS; + reference_hexa_permutated_SCCS[0][0] = 1; + reference_hexa_permutated_SCCS[0][1] = 0; + reference_hexa_permutated_SCCS[0][2] = 0; + reference_hexa_permutated_SCCS[1][0] = 0; + reference_hexa_permutated_SCCS[1][1] = 0; + reference_hexa_permutated_SCCS[1][2] = 1; + reference_hexa_permutated_SCCS[2][0] = 0; + reference_hexa_permutated_SCCS[2][1] = 1; + reference_hexa_permutated_SCCS[2][2] = 0; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("unpermutated_SCCS[" << iii << "] = " << unpermutated_SCCS[iii][0] << ":" << unpermutated_SCCS[iii][1] << ":" << unpermutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_unpermutated_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_unpermutated_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + std::array,7 > norms = aspect::Particle::Property::Utilities::compute_elastic_tensor_SCCS_decompositions(unpermutated_SCCS, full_elastic_matrix); + std::array totals = {{0,0,0}}; + for (size_t iii = 0; iii < 7; iii++) + { + // also check that the total is equal to the full norm. + + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("i:j = " << iii << ":" << jjj ); + CHECK(reference_norms[iii][jjj] == Approx(norms[iii][jjj])); + if (iii < 6) + { + totals[jjj] += norms[iii][jjj]; + } + } + } + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("j = " << jjj ); + CHECK(totals[jjj] == Approx(norms[6][jjj])); + } + + // also check that all the isotropic and totals are the same + for (size_t iii = 1; iii < 3; iii++) + { + CHECK(norms[5][0] == Approx(norms[5][iii])); + CHECK(norms[6][0] == Approx(norms[6][iii])); + } + + + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = aspect::Particle::Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + Tensor<2,3> hexa_permutated_SCCS; + for (size_t index = 0; index < 3; index++) + { + hexa_permutated_SCCS[index] = unpermutated_SCCS[permutation[index]]; + } + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("hexa_permutated_SCCS[" << iii << "] = " << hexa_permutated_SCCS[iii][0] << ":" << hexa_permutated_SCCS[iii][1] << ":" << hexa_permutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_hexa_permutated_SCCS[iii][0] == Approx(hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(hexa_permutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_hexa_permutated_SCCS[iii][0] == Approx(-hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(-hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(-hexa_permutated_SCCS[iii][2]).epsilon(epsilon)))); + } + } + + { + // This is testing some other elastic matrices of which the properties can + // be resonaly predicted (for example, mostly hexagonal) and some coverage + // testing: Case 4: from http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php (Be). + // This one should have at least in one direction almost all the the anistropy be + // hexagonal anisotropy. + + SymmetricTensor<2,6> full_elastic_matrix; + full_elastic_matrix[0][0] = 192.3; + full_elastic_matrix[0][1] = 26.7; + full_elastic_matrix[0][2] = 14; + full_elastic_matrix[1][1] = full_elastic_matrix[0][0]; + full_elastic_matrix[1][2] = full_elastic_matrix[0][2]; + full_elastic_matrix[2][2] = 336.4; + full_elastic_matrix[3][3] = 162.5; + full_elastic_matrix[4][4] = full_elastic_matrix[3][3]; + full_elastic_matrix[5][5] = 0.5*(full_elastic_matrix[0][0] + full_elastic_matrix[0][1]); + + std::array,7 > reference_norms; + reference_norms[0] = {{0,0,0}}; // no triclinic + reference_norms[1] = {{0,0,0}}; // no monoclinic + reference_norms[2] = {{16161.695,0.0,16161.695}}; // orthorhomic dependent on direction + reference_norms[3] = {{2786.31125,1425.78,2786.31125}}; // tetragonal could be present + reference_norms[4] = {{8079.22575,25601.452,8079.22575}}; // hexagonal is expected + reference_norms[5] = {{421517.088,421517.088,421517.088}}; // isotropic should be the same for every permutation + reference_norms[6] = {{448544.32,448544.32,448544.32}}; // total should be the same for every direction + + // It is already in a SCCS frame, so they should be + // unit vectors. The order is irrelevant. + Tensor<2,3> reference_unpermutated_SCCS; + reference_unpermutated_SCCS[0][0] = 0; + reference_unpermutated_SCCS[0][1] = 0; + reference_unpermutated_SCCS[0][2] = 1; + reference_unpermutated_SCCS[1][0] = 1; + reference_unpermutated_SCCS[1][1] = 0; + reference_unpermutated_SCCS[1][2] = 0; + reference_unpermutated_SCCS[2][0] = 0; + reference_unpermutated_SCCS[2][1] = 1; + reference_unpermutated_SCCS[2][2] = 0; + + // Permutation 2 has the highest norm, so reference_unpermutated_SCCS + // is permutated by {1,2,0}; + Tensor<2,3> reference_hexa_permutated_SCCS; + reference_hexa_permutated_SCCS[0][0] = 1; + reference_hexa_permutated_SCCS[0][1] = 0; + reference_hexa_permutated_SCCS[0][2] = 0; + reference_hexa_permutated_SCCS[1][0] = 0; + reference_hexa_permutated_SCCS[1][1] = 1; + reference_hexa_permutated_SCCS[1][2] = 0; + reference_hexa_permutated_SCCS[2][0] = 0; + reference_hexa_permutated_SCCS[2][1] = 0; + reference_hexa_permutated_SCCS[2][2] = 1; + + const SymmetricTensor<2,3> dilatation_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_dilatation_stiffness_tensor(full_elastic_matrix); + const SymmetricTensor<2,3> voigt_stiffness_tensor_full = aspect::Particle::Property::Utilities::compute_voigt_stiffness_tensor(full_elastic_matrix); + Tensor<2,3> unpermutated_SCCS = aspect::Particle::Property::Utilities::compute_unpermutated_SCCS(dilatation_stiffness_tensor_full, voigt_stiffness_tensor_full); + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("unpermutated_SCCS[" << iii << "] = " << unpermutated_SCCS[iii][0] << ":" << unpermutated_SCCS[iii][1] << ":" << unpermutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_unpermutated_SCCS[iii][0] == Approx(unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(unpermutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_unpermutated_SCCS[iii][0] == Approx(-unpermutated_SCCS[iii][0]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][1] == Approx(-unpermutated_SCCS[iii][1]).epsilon(epsilon) + && reference_unpermutated_SCCS[iii][2] == Approx(-unpermutated_SCCS[iii][2]).epsilon(epsilon)))); + } + std::array,7 > norms = aspect::Particle::Property::Utilities::compute_elastic_tensor_SCCS_decompositions(unpermutated_SCCS, full_elastic_matrix); + std::array totals = {{0,0,0}}; + for (size_t iii = 0; iii < 7; iii++) + { + // also check that the total is equal to the full norm. + + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("i:j = " << iii << ":" << jjj ); + CHECK(reference_norms[iii][jjj] == Approx(norms[iii][jjj])); + if (iii < 6) + { + totals[jjj] += norms[iii][jjj]; + } + } + } + for (size_t jjj = 0; jjj < 3; jjj++) + { + INFO("j = " << jjj ); + CHECK(totals[jjj] == Approx(norms[6][jjj])); + } + + // also check that all the isotropic and totals are the same + for (size_t iii = 1; iii < 3; iii++) + { + CHECK(norms[5][0] == Approx(norms[5][iii])); + CHECK(norms[6][0] == Approx(norms[6][iii])); + } + + + const size_t max_hexagonal_element_index = std::max_element(norms[4].begin(),norms[4].end())-norms[4].begin(); + std::array permutation = aspect::Particle::Property::Utilities::indexed_even_permutation(max_hexagonal_element_index); + Tensor<2,3> hexa_permutated_SCCS; + for (size_t index = 0; index < 3; index++) + { + hexa_permutated_SCCS[index] = unpermutated_SCCS[permutation[index]]; + } + for (size_t iii = 0; iii < 3; iii++) + { + const double epsilon = 1e-3; + INFO("hexa_permutated_SCCS[" << iii << "] = " << hexa_permutated_SCCS[iii][0] << ":" << hexa_permutated_SCCS[iii][1] << ":" << hexa_permutated_SCCS[iii][2]); + // The direction of the eigenvectors is not important, but they do need to be consistent + CHECK((reference_hexa_permutated_SCCS[iii][0] == Approx(hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(hexa_permutated_SCCS[iii][2]).epsilon(epsilon) + || (reference_hexa_permutated_SCCS[iii][0] == Approx(-hexa_permutated_SCCS[iii][0]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][1] == Approx(-hexa_permutated_SCCS[iii][1]).epsilon(epsilon) + && reference_hexa_permutated_SCCS[iii][2] == Approx(-hexa_permutated_SCCS[iii][2]).epsilon(epsilon)))); + } + } +}