add symmetric_independent_component_matrix function to utilities.

Author: MFraters

include/aspect/utilities.h

@@ -1060,6 +1060,13 @@ namespace aspect
      * entry in the list must match one of the allowed operations.
      */
     std::vector<Operator> create_model_operator_list(const std::vector<std::string> &operator_names);
+
+    /**
+     * Create matrix with unit at independent indices
+     */
+    template <int dim>
+    SymmetricTensor<2,dim> symmetric_independent_component_matrix (const unsigned int k);
+
   }
 }
 

source/utilities.cc

@@ -2532,6 +2532,21 @@ namespace aspect
     }
 
 
+    template <int dim>
+    SymmetricTensor<2,dim>
+    symmetric_independent_component_matrix (const unsigned int k)
+    {
+      // somehow I need this, because otherwise it complains that I am passing 3 arguments and it only takes 2...
+      const bool boolian = k < SymmetricTensor<2,dim>::n_independent_components;
+      Assert(boolian, ExcMessage("The component is larger then the amount of independent components in the matrix.") );
+
+      const TableIndices<2> indices_ij = SymmetricTensor<2,dim>::unrolled_to_component_indices (k);
+
+      Tensor<2,dim> result;
+      result[indices_ij] = 1;
+
+      return symmetrize(result);
+    }
 
 
 // Explicit instantiations
@@ -2618,5 +2633,8 @@ namespace aspect
     template std_cxx11::array<double,2> convert_point_to_array<2>(const Point<2> &point);
     template std_cxx11::array<double,3> convert_point_to_array<3>(const Point<3> &point);
 
+    template SymmetricTensor<2,2> symmetric_independent_component_matrix (const unsigned int k);
+    template SymmetricTensor<2,3> symmetric_independent_component_matrix (const unsigned int k);
+
   }
 }