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);
       }