Skip to content

Commit

Permalink
Further adress Bob's comments.
Browse files Browse the repository at this point in the history
  • Loading branch information
@MFraters
MFraters committed Aug 12, 2023
1 parent e799ff0 commit f8a4bcd
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 34 deletions.
32 changes: 23 additions & 9 deletions include/aspect/particle/property/elastic_tensor_decomposition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,38 +101,49 @@ namespace aspect
need_update () const override;

/**
* Return which data has to be provided to update the property.
* The integrated strains needs the gradients of the velocity.
*/
UpdateFlags
get_needed_update_flags () const override;

/**
* Return a specific permutation based on an index
* Return an even permutation based on an index. This is an internal
* utility function, also used by the unit tester.
*/
static
std::array<unsigned short int, 3>
indexed_even_permutation(const unsigned short int index);

/**
* Computes the voigt stiffness tensor from the elastic tensor
* 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
*/
static
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.
*/
static
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
*/
static
Tensor<2,3> compute_unpermutated_SCCS(const SymmetricTensor<2,3> &dilatation_stiffness_tensor,
Expand Down Expand Up @@ -179,6 +190,9 @@ namespace aspect
parse_parameters (ParameterHandler &prm) override;


/**
* The tenors below can be used to project matrices to different symmetries.
*/
static SymmetricTensor<2,21> projection_matrix_tric_to_mono;
static SymmetricTensor<2,9> projection_matrix_mono_to_ortho;
static SymmetricTensor<2,9> projection_matrix_ortho_to_tetra;
Expand Down
39 changes: 14 additions & 25 deletions source/particle/property/elastic_tensor_decomposition.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,14 +162,12 @@ namespace aspect
ExcMessage("No cpo property plugin found."));
AssertThrow(manager.plugin_name_exists("cpo elastic tensor"),
ExcMessage("No cpo elastic tensor property plugin found."));
Assert(manager.plugin_name_exists("decompose elastic matrix"),
ExcMessage("No hexagonal axes property plugin found."));

AssertThrow(manager.check_plugin_order("crystal preferred orientation","decompose elastic matrix"),
ExcMessage("To use the decompose elastic matrix plugin, the cpo plugin need to be defined before this plugin."));
AssertThrow(manager.check_plugin_order("crystal preferred orientation","elastic tensor decomposition"),
ExcMessage("To use the elastic tensor decomposition plugin, the cpo plugin need to be defined before this plugin."));

AssertThrow(manager.check_plugin_order("cpo elastic tensor","decompose elastic matrix"),
ExcMessage("To use the decompose elastic matrix plugin, the cpo elastic tensor plugin need 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 need 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"));
}
Expand Down Expand Up @@ -487,15 +485,15 @@ namespace aspect
// replaced by the lines below.
//auto ortho_and_higher_vector = projection_matrix_mono_to_ortho*mono_and_higher_vector;
dealii::Tensor<1,9> mono_and_higher_vector_cropped;
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];
mono_and_higher_vector_cropped[0] = mono_and_higher_vector[0];
mono_and_higher_vector_cropped[1] = mono_and_higher_vector[1];
mono_and_higher_vector_cropped[2] = mono_and_higher_vector[2];
mono_and_higher_vector_cropped[3] = mono_and_higher_vector[3];
mono_and_higher_vector_cropped[4] = mono_and_higher_vector[4];
mono_and_higher_vector_cropped[5] = mono_and_higher_vector[5];
mono_and_higher_vector_cropped[6] = mono_and_higher_vector[6];
mono_and_higher_vector_cropped[7] = mono_and_higher_vector[7];
mono_and_higher_vector_cropped[8] = mono_and_higher_vector[8];
dealii::Tensor<1,9> ortho_and_higher_vector;
ortho_and_higher_vector[0] = mono_and_higher_vector[0];
ortho_and_higher_vector[1] = mono_and_higher_vector[1];
Expand All @@ -506,7 +504,7 @@ namespace aspect
ortho_and_higher_vector[6] = mono_and_higher_vector[6];
ortho_and_higher_vector[7] = mono_and_higher_vector[7];
ortho_and_higher_vector[8] = mono_and_higher_vector[8];
auto mono_vector = mono_and_higher_vector_croped-ortho_and_higher_vector;
auto mono_vector = mono_and_higher_vector_cropped-ortho_and_higher_vector;
norms[1][permutation_i] = mono_vector.norm_square();


Expand Down Expand Up @@ -539,15 +537,6 @@ namespace aspect



template <int dim>
UpdateFlags
ElasticTensorDecomposition<dim>::get_needed_update_flags () const
{
return update_default;
}



template <int dim>
std::vector<std::pair<std::string, unsigned int>>
ElasticTensorDecomposition<dim>::get_property_information() const
Expand Down

0 comments on commit f8a4bcd

Please sign in to comment.