]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move closer to general symmetric tensors.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 29 Mar 2005 15:46:41 +0000 (15:46 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 29 Mar 2005 15:46:41 +0000 (15:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@10288 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/symmetric_tensor.h

index c52b1bd999b822874cdcac65e5593c6f7c6f0663..d8a3b40889c22d9d21490158c64abeb5b80a965e 100644 (file)
 
 
 #include <base/tensor.h>
+#include <base/table_indices.h>
 
 
 template <int rank, int dim> class SymmetricTensor;
-template <int dim> class SymmetricTensor<2,dim>;
 
 
 namespace internal
@@ -122,9 +122,7 @@ namespace internal
     template <int rank, int dim>
     struct AccessorTypes<rank, dim,true>
     {
-        typedef
-        const typename StorageType<rank,dim>::base_tensor_type
-        base_tensor_type;
+        typedef const ::SymmetricTensor<rank,dim> tensor_type;
 
         typedef double reference;
     };
@@ -141,86 +139,95 @@ namespace internal
     template <int rank, int dim>
     struct AccessorTypes<rank,dim,false>
     {
-        typedef
-        typename StorageType<rank,dim>::base_tensor_type
-        base_tensor_type;
+        typedef ::SymmetricTensor<rank,dim> tensor_type;
 
         typedef double &reference;
     };
 
+
+    template <int rank, int dim, bool constness>
+    class Accessor;
     
-    namespace Rank2Accessors
+                                     /**
+                                      * Accessor class to access the elements
+                                      * of individual rows in a symmetric
+                                      * tensor of rank 2. Since the elements
+                                      * of symmetric tensors are not stored as
+                                      * in a table, the accessors are a little
+                                      * more involved. However, for tensors of
+                                      * rank 2 they are still relatively
+                                      * simple in that an accessor is created
+                                      * by the SymmetricTensor class with the
+                                      * first access to <tt>operator[]</tt>;
+                                      * the accessor thereby points to a row
+                                      * of the tensor. Calling
+                                      * <tt>operator[]</tt> on the accessor
+                                      * then selects an entry of this
+                                      * row. Note that if this entry is not
+                                      * actually stored, then the transpose
+                                      * entry is chosen as that is guaranteed
+                                      * to be stored.
+                                      *
+                                      * @author Wolfgang Bangerth, 2005
+                                      */
+    template <int dim, bool constness>
+    class Accessor<2,dim,constness> 
     {
+      public:
+                                         /**
+                                          * Import which tensor we work on.
+                                          */
+        typedef
+        typename AccessorTypes<2,dim,constness>::tensor_type
+        tensor_type;
+
+                                         /**
+                                          * The type of a reference to an
+                                          * individual element of the
+                                          * symmetric tensor. If the tensor
+                                          * is constant, we can only return
+                                          * a value instead of a reference.
+                                          */
+        typedef typename AccessorTypes<2,dim,constness>::reference reference;
+
+                                         /**
+                                          * Constructor. Take the tensor to
+                                          * access as well as the row we
+                                          * point to as arguments.
+                                          */
+        Accessor (tensor_type        &tensor,
+                  const unsigned int  row);
 
-                                       /**
-                                        * Accessor class to access the
-                                        * elements of individual rows in a
-                                        * symmetric tensor. Since the elements
-                                        * of symmetric tensors are not stored
-                                        * as in a table, the accessors are a
-                                        * little more involved.
-                                        *
-                                        * @author Wolfgang Bangerth, 2005
-                                        */
-      template <int dim, bool constness>
-      class Accessor 
-      {
-        public:
-                                           /**
-                                            * Import which tensor we work on.
-                                            */
-          typedef
-          typename AccessorTypes<2,dim,constness>::base_tensor_type
-          base_tensor_type;
-
-                                           /**
-                                            * The type of a reference to an
-                                            * individual element of the
-                                            * symmetric tensor. If the tensor
-                                            * is constant, we can only return
-                                            * a value instead of a reference.
-                                            */
-          typedef typename AccessorTypes<2,dim,constness>::reference reference;
-
-                                           /**
-                                            * Constructor. Take the tensor to
-                                            * access as well as the row we
-                                            * point to as arguments.
-                                            */
-          Accessor (base_tensor_type  &tensor,
-                    const unsigned int  row);
-
-                                           /**
-                                            * Return a reference to an element
-                                            * of this row (if we point to a
-                                            * non-const tensor), or the value
-                                            * of the element (in case this is
-                                            * a constant tensor).
-                                            */
-          reference operator[] (const unsigned int column);
+                                         /**
+                                          * Return a reference to an element
+                                          * of this row (if we point to a
+                                          * non-const tensor), or the value
+                                          * of the element (in case this is
+                                          * a constant tensor).
+                                          */
+        reference operator[] (const unsigned int column);
           
-        private:
-                                           /**
-                                            * Reference to the tensor we
-                                            * access.
-                                            */
-          base_tensor_type &base_tensor;
-
-                                           /**
-                                            * Index of the row we access.
-                                            */
-          const unsigned int row;
-
-                                           /**
-                                            * Make the symmetric tensor
-                                            * classes a friend, since they are
-                                            * the only ones who can create
-                                            * objects like this.
-                                            */
-          template <int,int> class ::SymmetricTensor;
-      };
+      private:
+                                         /**
+                                          * Reference to the tensor we
+                                          * access.
+                                          */
+        tensor_type &tensor;
+
+                                         /**
+                                          * Index of the row we access.
+                                          */
+        const unsigned int row;
+
+                                         /**
+                                          * Make the symmetric tensor
+                                          * classes a friend, since they are
+                                          * the only ones who can create
+                                          * objects like this.
+                                          */
+        template <int,int> class ::SymmetricTensor;
+    };
 
-    }
   }
 }
 
@@ -228,19 +235,42 @@ namespace internal
 
 
 /**
- * Provide a class that stores symmetric tensors of rank 2 efficiently,
- * i.e. only store half of the off-diagonal elements of the full tensor.
+ * Provide a class that stores symmetric tensors of rank 2,4,... efficiently,
+ * i.e. only store those off-diagonal elements of the full tensor that are not
+ * redundant. For example, for symmetric 2x2 tensors, this would be the
+ * elements 11, 22, and 12, while the element 21 is equal to the 12 element.
+ *
+ * Using this class for symmetric tensors of rank 2 has advantages over
+ * matrices in many cases since the dimension is known to the compiler as well
+ * as the location of the data. It is therefore possible to produce far more
+ * efficient code than for matrices with runtime-dependent dimension. It is
+ * also more efficient than using the more general <tt>Tensor</tt> class,
+ * since less elements are stored, and the class automatically makes sure that
+ * the tensor represents a symmetric object.
+ *
+ * For tensors of higher rank, the savings in storage are even higher. For
+ * example for the 3x3x3x3 tensors of rank 4, only 36 instead of the full 81
+ * entries have to be stored.
+ *
+ * Tensors of rank 4 are considered symmetric if they are operators mapping
+ * symmetric rank-2 tensors onto symmetric rank-2 tensors. This entails
+ * certain symmetry properties on the elements in their 4-dimensional index
+ * space.
+ *
+ * Symmetric tensors are most often used in structural and fluid mechanics,
+ * where strains and stresses are usually symmetric tensors, and the
+ * stress-strain relationship is given by a symmetric rank-4 tensor.
  *
- * Using this tensor class for objects of rank 2 has advantages over
- * matrices in many cases since the dimension is known to the compiler
- * as well as the location of the data. It is therefore possible to
- * produce far more efficient code than for matrices with
- * runtime-dependent dimension.
+ * Note that symmetric tensors only exist with even numbers of indices. In
+ * other words, the only objects that you can use are
+ * <tt>SymmetricTensor<2,dim></tt>, <tt>SymmetricTensor<4,dim></tt>, etc, but
+ * <tt>SymmetricTensor<1,dim></tt> and <tt>SymmetricTensor<3,dim></tt> do not
+ * exist and their use will most likely lead to compiler errors.
  *
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
-class SymmetricTensor<2,dim>
+template <int rank, int dim>
+class SymmetricTensor
 {
   public:
                                     /**
@@ -261,12 +291,6 @@ class SymmetricTensor<2,dim>
                                      * data types.
                                      */
     static const unsigned int dimension = dim;
-
-                                    /**
-                                     * Publish the rank of this tensor to
-                                     * the outside world.
-                                     */
-    static const unsigned int rank      = 2;
     
                                      /**
                                       * Default constructor. Creates a zero
@@ -355,12 +379,33 @@ class SymmetricTensor<2,dim>
                                       */
     double operator * (const SymmetricTensor &s) const;
     
+                                     /**
+                                      * Return a read-write reference
+                                      * to the indicated element.
+                                      */
+    double & operator() (const TableIndices<rank> &indices);
+  
+                                     /**
+                                      * Return the value of the
+                                      * indicated element as a
+                                      * read-only reference.
+                                      *
+                                      * We return the requested value
+                                      * as a constant reference rather
+                                      * than by value since this
+                                      * object may hold data types
+                                      * that may be large, and we
+                                      * don't know here whether
+                                      * copying is expensive or not.
+                                      */
+    double operator() (const TableIndices<rank> &indices) const;
+
                                      /**
                                       * Access the elements of a row of this
                                       * symmetric tensor. This function is
                                       * called for constant tensors.
                                       */
-    internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true>
+    internal::SymmetricTensor::Accessor<rank,dim,true>
     operator [] (const unsigned int row) const;
 
                                      /**
@@ -368,7 +413,7 @@ class SymmetricTensor<2,dim>
                                       * symmetric tensor. This function is
                                       * called for non-constant tensors.
                                       */
-    internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false>
+    internal::SymmetricTensor::Accessor<rank,dim,false>
     operator [] (const unsigned int row);
     
                                      /**
@@ -424,74 +469,32 @@ namespace internal
 {
   namespace SymmetricTensor
   {
-    namespace Rank2Accessors
+    template <int dim, bool constness>
+    Accessor<2,dim,constness>::
+    Accessor (tensor_type       &tensor,
+              const unsigned int  row)
+                    :
+                    tensor (tensor),
+                    row (row)
+    {}
+
+
+
+    template <int dim, bool constness>
+    typename Accessor<2,dim,constness>::reference
+    Accessor<2,dim,constness>::
+    operator[] (const unsigned int column)
     {
-      template <int dim, bool constness>
-      Accessor<dim,constness>::
-      Accessor (base_tensor_type  &base_tensor,
-                   const unsigned int row)
-                      :
-                      base_tensor (base_tensor),
-                      row (row)
-      {
-        Assert (row < dim, ExcIndexRange (row, 0, dim));
-      }
-
-
-
-      template <int dim, bool constness>
-      typename Accessor<dim,constness>::reference
-      Accessor<dim,constness>::
-      operator[] (const unsigned int column)
-      {
-        Assert (column < dim, ExcIndexRange (column, 0, dim));
-
-                                         // first treat the main diagonal
-                                         // elements, which are stored
-                                         // consecutively at the beginning
-        if (row == column)
-          return base_tensor[row];
-
-                                         // the rest is messier and requires a
-                                         // few switches. if someone has a
-                                         // better idea, help is welcome
-        switch (dim)
-          {
-            case 2:
-                  Assert (((row==1) && (column==0)) || ((row==0) && (column==1)),
-                          ExcInternalError());
-                  return base_tensor[2];
-
-            case 3:
-                  if (((row==0) && (column==1)) ||
-                      ((row==1) && (column==0)))
-                    return base_tensor[3];
-                  else if (((row==0) && (column==2)) ||
-                           ((row==2) && (column==0)))
-                    return base_tensor[4];
-                  else if (((row==1) && (column==2)) ||
-                           ((row==2) && (column==1)))
-                    return base_tensor[5];
-                  else
-                    Assert (false, ExcInternalError());
-
-            default:
-                  Assert (false, ExcNotImplemented());
-          }
-
-        Assert (false, ExcInternalError());
-        static double dummy_but_referenceable = 0;
-        return dummy_but_referenceable;
-      }
+      return tensor(TableIndices<2> (row, column));
     }
   }
 }
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim>::SymmetricTensor ()
+SymmetricTensor<rank,dim>::SymmetricTensor ()
 {}
 
 
@@ -526,10 +529,10 @@ SymmetricTensor<2,3>::SymmetricTensor (const Tensor<2,3> &t)
 }
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator = (const SymmetricTensor<rank,dim> &t)
 {
   data = t.data;
   return *this;
@@ -537,30 +540,30 @@ SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t)
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
 bool
-SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const
+SymmetricTensor<rank,dim>::operator == (const SymmetricTensor<rank,dim> &t) const
 {
   return data == t.data;
 }
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
 bool
-SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const
+SymmetricTensor<rank,dim>::operator != (const SymmetricTensor<rank,dim> &t) const
 {
   return data != t.data;
 }
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator += (const SymmetricTensor<rank,dim> &t)
 {
   data += t.data;
   return *this;
@@ -568,10 +571,10 @@ SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t)
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator -= (const SymmetricTensor<rank,dim> &t)
 {
   data -= t.data;
   return *this;
@@ -579,10 +582,10 @@ SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t)
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator *= (const double d)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator *= (const double d)
 {
   data *= d;
   return *this;
@@ -590,10 +593,10 @@ SymmetricTensor<2,dim>::operator *= (const double d)
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim> &
-SymmetricTensor<2,dim>::operator /= (const double d)
+SymmetricTensor<rank,dim> &
+SymmetricTensor<rank,dim>::operator /= (const double d)
 {
   data /= d;
   return *this;
@@ -601,10 +604,10 @@ SymmetricTensor<2,dim>::operator /= (const double d)
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator + (const SymmetricTensor &t) const
 {
   SymmetricTensor tmp = *this;
   tmp.data += t.data;
@@ -613,10 +616,10 @@ SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator - (const SymmetricTensor &t) const
 {
   SymmetricTensor tmp = *this;
   tmp.data -= t.data;
@@ -625,10 +628,10 @@ SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim>
-SymmetricTensor<2,dim>::operator - () const
+SymmetricTensor<rank,dim>
+SymmetricTensor<rank,dim>::operator - () const
 {
   SymmetricTensor tmp = *this;
   tmp.data = -tmp.data;
@@ -637,59 +640,229 @@ SymmetricTensor<2,dim>::operator - () const
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
 void
-SymmetricTensor<2,dim>::clear ()
+SymmetricTensor<rank,dim>::clear ()
 {
   data.clear ();
 }
 
 
 
-template <int dim>
+template <int rank, int dim>
 inline
 unsigned int
-SymmetricTensor<2,dim>::memory_consumption ()
+SymmetricTensor<rank,dim>::memory_consumption ()
 {
-  return internal::SymmetricTensor::StorageType<2,dim>::memory_consumption ();
+  return
+    internal::SymmetricTensor::StorageType<rank,dim>::memory_consumption ();
 }
 
 
 
-template <int dim>
+template <>
 double
-SymmetricTensor<2,dim>::operator * (const SymmetricTensor &s) const
+SymmetricTensor<2,1>::operator * (const SymmetricTensor<2,1> &s) const
 {
-  double t = 0;
-  unsigned int i=0;
-  for (; i<dim; ++i)
-    t += data[i] * s.data[i];
+  return data[0] * s.data[0];
+}
 
-  for (; i<internal::SymmetricTensor::StorageType<2,dim>::n_independent_components; ++i)
-    t += 2 * data[i] * s.data[i];
 
-  return t;
+
+template <>
+double
+SymmetricTensor<2,2>::operator * (const SymmetricTensor<2,2> &s) const
+{
+  return (data[0] * s.data[0] +
+          data[1] * s.data[1] +
+          2*data[2] * s.data[2]);
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,3>::operator * (const SymmetricTensor<2,3> &s) const
+{
+  return (data[0] * s.data[0] +
+          data[1] * s.data[1] +
+          data[2] * s.data[2] +
+          2*data[3] * s.data[3] +
+          2*data[4] * s.data[4] +
+          2*data[5] * s.data[5]);
+}
+
+
+
+template <>
+double &
+SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices)
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 1;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+  return data[0];
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices) const
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 1;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+  return data[0];
+}
+
+
+
+template <>
+double &
+SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices)
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 2;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+                                   // first treat the main diagonal
+                                   // elements, which are stored
+                                   // consecutively at the beginning
+  if (indices[0] == indices[1])
+    return data[indices[0]];
+
+                                   // the rest is messier and requires a few
+                                   // switches. at least for the 2x2 case it
+                                   // is reasonably simple
+  Assert (((indices[0]==1) && (indices[1]==0)) ||
+          ((indices[0]==0) && (indices[1]==1)),
+          ExcInternalError());  
+  return data[2];
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices) const
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 2;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+                                   // first treat the main diagonal
+                                   // elements, which are stored
+                                   // consecutively at the beginning
+  if (indices[0] == indices[1])
+    return data[indices[0]];
+
+                                   // the rest is messier and requires a few
+                                   // switches. at least for the 2x2 case it
+                                   // is reasonably simple
+  Assert (((indices[0]==1) && (indices[1]==0)) ||
+          ((indices[0]==0) && (indices[1]==1)),
+          ExcInternalError());  
+  return data[2];
 }
+
+
+
+template <>
+double &
+SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices)
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 3;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+                                   // first treat the main diagonal
+                                   // elements, which are stored
+                                   // consecutively at the beginning
+  if (indices[0] == indices[1])
+    return data[indices[0]];
+
+                                   // the rest is messier and requires a few
+                                   // switches, but simpler if we just sort
+                                   // our indices
+  TableIndices<2> sorted_indices (indices);
+  sorted_indices.sort ();
   
+  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
+    return data[3];
+  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
+    return data[4];
+  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
+    return data[5];
+  else
+    Assert (false, ExcInternalError());
+
+  static double dummy_but_referenceable = 0;
+  return dummy_but_referenceable;
+}
+
 
 
-template <int dim>
-internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true>
-SymmetricTensor<2,dim>::operator [] (const unsigned int row) const
+template <>
+double
+SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices) const
+{
+  const unsigned int rank = 2;
+  const unsigned int dim  = 3;
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dim, ExcIndexRange (indices[r], 0, dim));
+
+                                   // first treat the main diagonal
+                                   // elements, which are stored
+                                   // consecutively at the beginning
+  if (indices[0] == indices[1])
+    return data[indices[0]];
+
+                                   // the rest is messier and requires a few
+                                   // switches, but simpler if we just sort
+                                   // our indices
+  TableIndices<2> sorted_indices (indices);
+  sorted_indices.sort ();
+  
+  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
+    return data[3];
+  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
+    return data[4];
+  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
+    return data[5];
+  else
+    Assert (false, ExcInternalError());
+
+  static double dummy_but_referenceable = 0;
+  return dummy_but_referenceable;
+}
+
+
+
+template <int rank, int dim>
+internal::SymmetricTensor::Accessor<rank,dim,true>
+SymmetricTensor<rank,dim>::operator [] (const unsigned int row) const
 {
   return
-    internal::SymmetricTensor::Rank2Accessors::Accessor<dim,true> (data, row);
+    internal::SymmetricTensor::Accessor<rank,dim,true> (*this, row);
 }
 
 
 
-template <int dim>
-internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false>
-SymmetricTensor<2,dim>::operator [] (const unsigned int row)
+template <int rank, int dim>
+internal::SymmetricTensor::Accessor<rank,dim,false>
+SymmetricTensor<rank,dim>::operator [] (const unsigned int row)
 {
   return
-    internal::SymmetricTensor::Rank2Accessors::Accessor<dim,false> (data, row);
+    internal::SymmetricTensor::Accessor<rank,dim,false> (*this, row);
 }
 
 
@@ -766,8 +939,8 @@ double determinant (const SymmetricTensor<2,3> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
-double trace (const SymmetricTensor<2,dim> &d)
+template <int rank, int dim>
+double trace (const SymmetricTensor<rank,dim> &d)
 {
   double t=0;
   for (unsigned int i=0; i<dim; ++i)
@@ -786,10 +959,10 @@ double trace (const SymmetricTensor<2,dim> &d)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int rank, int dim>
 inline
-SymmetricTensor<2,dim>
-transpose (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<rank,dim>
+transpose (const SymmetricTensor<rank,dim> &t)
 {
   return t;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.