Use constructor and unroll functions for Tensors

include/aspect/utilities.h

@@ -1251,7 +1251,7 @@ namespace aspect
       template <int dim>
       inline
       SymmetricTensor<2,dim>
-      to_symmetric_tensor(const ArrayView<double> &values)
+      to_symmetric_tensor(const ArrayView<const double> &values)
       {
         // Assemble stress tensor if elastic behavior is enabled
         SymmetricTensor<2,dim> output;

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];
+                ArrayView<const double> stress_values(&in.composition[i][0], SymmetricTensor<2,dim>::n_independent_components);
+                const SymmetricTensor<2,dim> stress_old = Utilities::Tensors::to_symmetric_tensor<dim>(stress_values);
 
                 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 ArrayView<const double> stress_values(&in.composition[i][0], SymmetricTensor<2,dim>::n_independent_components);
+                const SymmetricTensor<2,dim> stress_old(Utilities::Tensors::to_symmetric_tensor<dim>(stress_values));
 
                 // Calculate the rotated stresses
                 // Rotation (vorticity) tensor (equation 25 in Moresi et al., 2003, J. Comp. Phys.)
@@ -426,11 +424,11 @@ 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
+                ArrayView<double> stress_update_values(&out.reaction_terms[i][0], SymmetricTensor<2,dim>::n_independent_components);
+                Utilities::Tensors::unroll(stress_update, stress_update_values);
               }
           }
       }

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;
+        ArrayView<const double> elastic_values(&data[cpo_data_position],SymmetricTensor<2,6>::n_independent_components);
+        return Utilities::Tensors::to_symmetric_tensor<6>(elastic_values);
       }
 
 
@@ -188,8 +186,8 @@ 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)];
+        const ArrayView<double> elastic_values(&data[cpo_data_position],SymmetricTensor<2,6>::n_independent_components);
+        Utilities::Tensors::unroll(elastic_tensor, data);
       }
 
 

source/particle/property/integrated_strain.cc

@@ -43,11 +43,8 @@ namespace aspect
                                                       const std::vector<Tensor<1,dim>> &gradients,
                                                       typename ParticleHandler<dim>::particle_iterator &particle) const
       {
-        const auto data = particle->get_properties();
-
-        Tensor<2,dim> old_strain;
-        for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
-          old_strain[Tensor<2,dim>::unrolled_to_component_indices(i)] = data[data_position + i];
+        ArrayView<double> strain_values(&particle->get_properties()[data_position], Tensor<2,dim>::n_independent_components);
+        Tensor<2,dim> old_strain(strain_values);
 
         Tensor<2,dim> grad_u;
         for (unsigned int d=0; d<dim; ++d)
@@ -75,8 +72,7 @@ namespace aspect
         // strain of the current time step
         new_strain = old_strain + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
 
-        for (unsigned int i = 0; i < Tensor<2,dim>::n_independent_components ; ++i)
-          data[data_position + i] = new_strain[Tensor<2,dim>::unrolled_to_component_indices(i)];
+        new_strain.unroll(strain_values.begin(),strain_values.end());
       }
 
       template <int dim>