Merge pull request #5640 from gassmoeller/add_unroll_function_symmetr…

…ic_tensor

Simplify creating and unrolling symmetric tensors

include/aspect/utilities.h

@@ -1245,6 +1245,47 @@ namespace aspect
  */
  namespace Tensors
  {
+ /**
+ * Convert a series of doubles to a SymmetricTensor of the same size.
+ * The input is expected to be ordered according to the
+ * SymmetricTensor::unrolled_to_component_indices() function.
+ */
+ template <int dim, class Iterator>
+ inline
+ SymmetricTensor<2,dim>
+ to_symmetric_tensor(const Iterator begin,
+ const Iterator end)
+ {
+ AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
+ (void) end;
+
+ SymmetricTensor<2,dim> output;
+
+ Iterator next = begin;
+ for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
+ output[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)] = *next;
+
+ return output;
+ }
+
+ /**
+ * Unroll a SymmetricTensor into a series of doubles. The output is ordered
+ * according to the SymmetricTensor::unrolled_to_component_indices() function.
+ */
+ template <int dim, class Iterator>
+ inline
+ void
+ unroll_symmetric_tensor_into_array(const SymmetricTensor<2,dim> &tensor,
+ const Iterator begin,
+ const Iterator end)
+ {
+ AssertDimension(std::distance(begin, end), (SymmetricTensor<2,dim>::n_independent_components));
+ (void) end;
+
+ Iterator next = begin;
+ for (unsigned int i=0; i < SymmetricTensor<2,dim>::n_independent_components; ++i, ++next)
+ *next = tensor[SymmetricTensor<2,dim>::unrolled_to_component_indices(i)];
+ }
 
  /**
  * Rotate a 3D 4th order tensor representing the full stiffnexx matrix using a 3D 2nd order rotation tensor

source/material_model/rheology/elasticity.cc

@@ -314,9 +314,8 @@ namespace aspect
  for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
  {
  // Get old stresses from compositional fields
- SymmetricTensor<2,dim> stress_old;
- for (unsigned int j=0; j < SymmetricTensor<2,dim>::n_independent_components; ++j)
- stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j];
+ const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][0],
+ &in.composition[i][0]+SymmetricTensor<2,dim>::n_independent_components));
 
  elastic_out->elastic_force[i] = -effective_creep_viscosities[i] / calculate_elastic_viscosity(average_elastic_shear_moduli[i]) * stress_old;
  // The viscoelastic strain rate is needed only when the Newton method is selected.
@@ -393,9 +392,8 @@ namespace aspect
  for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
  {
  // Get old stresses from compositional fields
- SymmetricTensor<2,dim> stress_old;
- for (unsigned int j=0; j < SymmetricTensor<2,dim>::n_independent_components; ++j)
- stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)] = in.composition[i][j];
+ const SymmetricTensor<2,dim> stress_old(Utilities::Tensors::to_symmetric_tensor<dim>(&in.composition[i][0],
+ &in.composition[i][0]+SymmetricTensor<2,dim>::n_independent_components));
 
  // Calculate the rotated stresses
  // Rotation (vorticity) tensor (equation 25 in Moresi et al., 2003, J. Comp. Phys.)
@@ -426,11 +424,12 @@ namespace aspect
  // timestep, then no averaging occurs as dt/dte = 1.
  stress_new = ( ( 1. - ( dt / dte ) ) * stress_old ) + ( ( dt / dte ) * stress_new ) ;
 
- // Fill reaction terms
- for (unsigned int j = 0; j < SymmetricTensor<2,dim>::n_independent_components ; ++j)
- out.reaction_terms[i][j] = -stress_old[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)]
- + stress_new[SymmetricTensor<2,dim>::unrolled_to_component_indices(j)];
+ const SymmetricTensor<2,dim> stress_update = stress_new - stress_old;
 
+ // Fill reaction terms
+ Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_update,
+ &out.reaction_terms[i][0],
+ &out.reaction_terms[i][0]+SymmetricTensor<2,dim>::n_independent_components);
  }
  }
  }

source/particle/property/cpo_elastic_tensor.cc

@@ -174,10 +174,8 @@ namespace aspect
  CpoElasticTensor<dim>::get_elastic_tensor(unsigned int cpo_data_position,
  const ArrayView<double> &data)
  {
- SymmetricTensor<2,6> elastic_tensor;
- for (unsigned int i = 0; i < SymmetricTensor<2,6>::n_independent_components ; ++i)
- elastic_tensor[SymmetricTensor<2,6>::unrolled_to_component_indices(i)] = data[cpo_data_position + i];
- return elastic_tensor;
+ return Utilities::Tensors::to_symmetric_tensor<6>(&data[cpo_data_position],
+ &data[cpo_data_position]+SymmetricTensor<2,6>::n_independent_components);
  }
 
 
@@ -188,8 +186,9 @@ namespace aspect
  const ArrayView<double> &data,
  const SymmetricTensor<2,6> &elastic_tensor)
  {
- for (unsigned int i = 0; i < SymmetricTensor<2,6>::n_independent_components ; ++i)
- data[cpo_data_position + i] = elastic_tensor[SymmetricTensor<2,6>::unrolled_to_component_indices(i)];
+ Utilities::Tensors::unroll_symmetric_tensor_into_array(elastic_tensor,
+ &data[cpo_data_position],
+ &data[cpo_data_position]+SymmetricTensor<2,6>::n_independent_components);
  }