]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a class implementing symmetric tensors of rank 2.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 28 Mar 2005 17:31:03 +0000 (17:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 28 Mar 2005 17:31:03 +0000 (17:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@10249 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/symmetric_tensor.h [new file with mode: 0644]

diff --git a/deal.II/base/include/base/symmetric_tensor.h b/deal.II/base/include/base/symmetric_tensor.h
new file mode 100644 (file)
index 0000000..3af4dcd
--- /dev/null
@@ -0,0 +1,625 @@
+//----------------------------  symmetric_tensor.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 by the deal authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  symmetric_tensor.h  ---------------------------
+#ifndef __deal2__symmetric_tensor_h
+#define __deal2__symmetric_tensor_h
+
+
+template <int rank, int dim> class SymmetricTensor;
+
+
+
+/**
+ * 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.
+ *
+ * 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.
+ *
+ * @author Wolfgang Bangerth, 2005
+ */
+template <int dim>
+class SymmetricTensor<2,dim>
+{
+  public:
+                                    /**
+                                     * Provide a way to get the
+                                     * dimension of an object without
+                                     * explicit knowledge of it's
+                                     * data type. Implementation is
+                                     * this way instead of providing
+                                     * a function <tt>dimension()</tt>
+                                     * because now it is possible to
+                                     * get the dimension at compile
+                                     * time without the expansion and
+                                     * preevaluation of an inlined
+                                     * function; the compiler may
+                                     * therefore produce more
+                                     * efficient code and you may use
+                                     * this value to declare other
+                                     * 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
+                                      * tensor.
+                                      */
+    SymmetricTensor ();
+
+                                     /**
+                                      * Constructor. Generate a symmetric
+                                      * tensor from a general one. Assumes
+                                      * that @p t is already symmetric, but
+                                      * this is not checked: we simply copy
+                                      * only a subset of elements.
+                                      */
+    SymmetricTensor (const Tensor<2,dim> &t);
+
+                                    /**
+                                     *  Assignment operator.
+                                     */
+    SymmetricTensor & operator = (const SymmetricTensor &);
+
+                                    /**
+                                     *  Test for equality of two tensors.
+                                     */
+    bool operator == (const SymmetricTensor &) const;
+
+                                    /**
+                                     *  Test for inequality of two tensors.
+                                     */
+    bool operator != (const SymmetricTensor &) const;
+
+                                    /**
+                                     *  Add another tensor.
+                                     */
+    SymmetricTensor & operator += (const SymmetricTensor &);
+    
+                                    /**
+                                     *  Subtract another tensor.
+                                     */
+    SymmetricTensor & operator -= (const SymmetricTensor &);
+
+                                    /**
+                                     *  Scale the tensor by <tt>factor</tt>,
+                                     *  i.e. multiply all components by
+                                     *  <tt>factor</tt>.
+                                     */
+    SymmetricTensor & operator *= (const double factor);
+
+                                    /**
+                                     *  Scale the vector by
+                                     *  <tt>1/factor</tt>.
+                                     */
+    SymmetricTensor & operator /= (const double factor);
+
+                                    /**
+                                     *  Add two tensors. If possible, you
+                                     *  should use <tt>operator +=</tt>
+                                     *  instead since this does not need the
+                                     *  creation of a temporary.
+                                     */
+    SymmetricTensor   operator + (const SymmetricTensor &) const;
+
+                                    /**
+                                     *  Subtract two tensors. If possible,
+                                     *  you should use <tt>operator -=</tt>
+                                     *  instead since this does not need the
+                                     *  creation of a temporary.
+                                     */
+    SymmetricTensor   operator - (const SymmetricTensor &) const;
+
+                                    /**
+                                     * Unary minus operator. Negate all
+                                     * entries of a tensor.
+                                     */
+    SymmetricTensor   operator - () const;
+
+                                     /**
+                                      * Return the Frobenius-norm of a tensor,
+                                      * i.e. the square root of the sum of
+                                      * squares of all entries.
+                                      */
+    double norm () const;
+    
+                                    /**
+                                     * Reset all values to zero.
+                                     *
+                                     * Note that this is partly inconsistent
+                                     * with the semantics of the @p clear()
+                                     * member functions of the STL and of
+                                     * several other classes within deal.II
+                                     * which not only reset the values of
+                                     * stored elements to zero, but release
+                                     * all memory and return the object into
+                                     * a virginial state. However, since the
+                                     * size of objects of the present type is
+                                     * determined by its template parameters,
+                                     * resizing is not an option, and indeed
+                                     * the state where all elements have a
+                                     * zero value is the state right after
+                                     * construction of such an object.
+                                     */
+    void clear ();
+
+                                    /**
+                                     * Determine an estimate for
+                                     * the memory consumption (in
+                                     * bytes) of this
+                                     * object.
+                                     */
+    static unsigned int memory_consumption ();
+
+    
+  private:
+                                     /**
+                                      * Number of independent components of a
+                                      * symmetric tensor of rank 2. We store
+                                      * only the upper right half of it.
+                                      */
+    static const unsigned int
+    n_tensor_components = (dim*dim + dim)/2;
+
+                                     /**
+                                      * Declare the type in which we actually
+                                      * store the data.
+                                      */
+    typedef Tensor<1,n_tensor_components> StorageType;
+    
+                                     /**
+                                      * Data storage for a symmetric tensor.
+                                      */
+    StorageType data;
+};
+
+
+
+// ------------------------- inline functions ------------------------
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>::SymmetricTensor ()
+{}
+
+
+
+template <>
+inline
+SymmetricTensor<2,2>::SymmetricTensor (const Tensor<2,2> &t)
+{
+  Assert (t[0][1] == t[1][0], ExcInternalError());
+
+  data[0] = t[0][0];
+  data[1] = t[1][1];
+  data[2] = t[0][1];
+}
+
+
+
+template <>
+inline
+SymmetricTensor<2,3>::SymmetricTensor (const Tensor<2,3> &t)
+{
+  Assert (t[0][1] == t[1][0], ExcInternalError());
+  Assert (t[0][2] == t[2][0], ExcInternalError());
+  Assert (t[1][2] == t[2][1], ExcInternalError());
+  
+  data[0] = t[0][0];
+  data[1] = t[1][1];
+  data[2] = t[2][2];
+  data[3] = t[0][1];
+  data[4] = t[0][2];
+  data[5] = t[1][2];
+}
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator = (const SymmetricTensor<2,dim> &t)
+{
+  data = t.data;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+bool
+SymmetricTensor<2,dim>::operator == (const SymmetricTensor<2,dim> &t) const
+{
+  return data == t.data;
+}
+
+
+
+template <int dim>
+inline
+bool
+SymmetricTensor<2,dim>::operator != (const SymmetricTensor<2,dim> &t) const
+{
+  return data != t.data;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator += (const SymmetricTensor<2,dim> &t)
+{
+  data += t.data;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator -= (const SymmetricTensor<2,dim> &t)
+{
+  data -= t.data;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator *= (const double d)
+{
+  data *= d;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim> &
+SymmetricTensor<2,dim>::operator /= (const double d)
+{
+  data /= d;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator + (const SymmetricTensor &t) const
+{
+  SymmetricTensor tmp = *this;
+  tmp.data += t.data;
+  return tmp;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator - (const SymmetricTensor &t) const
+{
+  SymmetricTensor tmp = *this;
+  tmp.data -= t.data;
+  return tmp;
+}
+
+
+
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+SymmetricTensor<2,dim>::operator - () const
+{
+  SymmetricTensor tmp = *this;
+  tmp.data = -tmp.data;
+  return tmp;
+}
+
+
+
+template <int dim>
+inline
+void
+SymmetricTensor<2,dim>::clear ()
+{
+  data.clear ();
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+SymmetricTensor<2,dim>::memory_consumption ()
+{
+  return StorageType::memory_consumption ();
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,1>::norm () const
+{
+  return std::sqrt(data[0]*data[0]);
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,2>::norm () const
+{
+  return std::sqrt(data[0]*data[0] + data[1]*data[1] + 2*data[2]*data[2]);
+}
+
+
+
+template <>
+double
+SymmetricTensor<2,3>::norm () const
+{
+  return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2] +
+                   2*data[3]*data[3] + 2*data[4]*data[4] + 2*data[5]*data[5]);
+}
+
+
+/* ----------------- Non-member functions operating on tensors. ------------ */
+
+/**
+ * Compute the determinant of a tensor of rank one and dimension
+ * one. Since this is a number, the return value is, of course, the
+ * number itself.
+ *
+ * @relates Tensor
+ * @author Wolfgang Bangerth, 2005
+ */
+inline
+double determinant (const SymmetricTensor<1,1> &t)
+{
+  Assert (false, ExcNotImplemented());
+  return 0;
+  
+//  return t[0];
+}
+
+
+
+/**
+ * Compute the determinant of a tensor or rank 2, here for <tt>dim==2</tt>.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+inline
+double determinant (const SymmetricTensor<2,2> &t)
+{
+  Assert (false, ExcNotImplemented());
+  return 0;
+  
+//   return ((t[0][0] * t[1][1]) -
+//       (t[1][0] * t[0][1]));
+}
+
+
+
+
+/**
+ * Compute the determinant of a tensor or rank 2, here for <tt>dim==3</tt>.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+inline
+double determinant (const SymmetricTensor<2,3> &t)
+{
+  Assert (false, ExcNotImplemented());
+  return 0;
+  
+//                                // get this using Maple:
+//                                // with(linalg);
+//                                // a := matrix(3,3);
+//                                // x := det(a);
+//                                // readlib(C);
+//                                // C(x, optimized);
+//   return ( t[0][0]*t[1][1]*t[2][2]
+//        -t[0][0]*t[1][2]*t[2][1]
+//        -t[1][0]*t[0][1]*t[2][2]
+//        +t[1][0]*t[0][2]*t[2][1]
+//        +t[2][0]*t[0][1]*t[1][2]
+//        -t[2][0]*t[0][2]*t[1][1] );
+}
+
+
+
+/**
+ * Compute and return the trace of a tensor of rank 2, i.e. the sum of
+ * its diagonal entries.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+template <int dim>
+double trace (const SymmetricTensor<2,dim> &d)
+{
+  Assert (false, ExcNotImplemented());
+  return 0;
+  
+//   double t=0;
+//   for (unsigned int i=0; i<dim; ++i)
+//     t += d[i][i];
+//   return t;
+}
+
+
+
+/**
+ * Compute and return the inverse of the given tensor. Since the
+ * compiler can perform the return value optimization, and since the
+ * size of the return object is known, it is acceptable to return the
+ * result by value, rather than by reference as a parameter.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+invert (const SymmetricTensor<2,dim> &t)
+{
+  Assert (false, ExcNotImplemented());
+  return SymmetricTensor<2,dim>();
+  
+//   SymmetricTensor<2,dim> return_tensor;
+//   switch (dim) 
+//     {
+//       case 1:
+//         return_tensor[0][0] = 1.0/t[0][0];
+//         return return_tensor;
+//       case 2:
+//                                          // this is Maple output,
+//                                          // thus a bit unstructured
+//       {
+//     const double t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]);
+//     return_tensor[0][0] = t[1][1]*t4;
+//     return_tensor[0][1] = -t[0][1]*t4;
+//     return_tensor[1][0] = -t[1][0]*t4;
+//     return_tensor[1][1] = t[0][0]*t4;
+//     return return_tensor;
+//       };
+      
+//       case 3:
+//       {
+//     const double t4 = t[0][0]*t[1][1],
+//                  t6 = t[0][0]*t[1][2],
+//                  t8 = t[0][1]*t[1][0],
+//                 t00 = t[0][2]*t[1][0],
+//                 t01 = t[0][1]*t[2][0],
+//                 t04 = t[0][2]*t[2][0],
+//                 t07 = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
+//                            t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
+//     return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07;
+//     return_tensor[0][1] = -(t[0][1]*t[2][2]-t[0][2]*t[2][1])*t07;
+//     return_tensor[0][2] = -(-t[0][1]*t[1][2]+t[0][2]*t[1][1])*t07;
+//     return_tensor[1][0] = -(t[1][0]*t[2][2]-t[1][2]*t[2][0])*t07;
+//     return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07;
+//     return_tensor[1][2] = -(t6-t00)*t07;
+//     return_tensor[2][0] = -(-t[1][0]*t[2][1]+t[1][1]*t[2][0])*t07;
+//     return_tensor[2][1] = -(t[0][0]*t[2][1]-t01)*t07;
+//     return_tensor[2][2] = (t4-t8)*t07;
+//     return return_tensor;
+//       };
+
+//                                     // if desired, take over the
+//                                     // inversion of a 4x4 tensor
+//                                     // from the FullMatrix
+       
+//       default:
+//         AssertThrow (false, ExcNotImplemented());
+//     };    
+//   return return_tensor;
+}
+
+
+
+/**
+ * Return the transpose of the given symmetric tensor. Since we are working
+ * with symmetric objects, the transpose is of course the same as the original
+ * tensor. This function mainly exists for compatibility with the Tensor
+ * class.
+ *
+ * @relates SymmetricTensor
+ * @author Wolfgang Bangerth, 2005
+ */
+template <int dim>
+inline
+SymmetricTensor<2,dim>
+transpose (const SymmetricTensor<2,dim> &t)
+{
+  return t;
+}
+
+
+/**
+ * Multiplication of a symmetric tensor of general rank with a scalar double
+ * from the right.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator * (const SymmetricTensor<rank,dim> &t,
+           const double            factor)
+{
+  SymmetricTensor<rank,dim> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Multiplication of a symmetric tensor of general rank with a scalar double
+ * from the left.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator * (const double            factor,
+           const SymmetricTensor<rank,dim> &t)
+{
+  SymmetricTensor<rank,dim> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Division of a symmetric tensor of general rank by a scalar double.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim>
+inline
+SymmetricTensor<rank,dim>
+operator / (const SymmetricTensor<rank,dim> &t,
+           const double            factor)
+{
+  SymmetricTensor<rank,dim> tt = t;
+  tt /= factor;
+  return tt;
+}
+
+
+
+#endif

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.