]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Eliminate specialization for Tensor<1,dim>
authorMatthias Maier <tamiko@43-1.org>
Wed, 2 Sep 2015 03:51:05 +0000 (22:51 -0500)
committerMatthias Maier <tamiko@43-1.org>
Mon, 7 Sep 2015 18:27:20 +0000 (13:27 -0500)
This is more or less an exact code copy of Tensor<rank_,dim>.

tensor_base.h now contains the full definition of Tensor<rank_,dim>. All
external functions operating on tensors that are not strictly needed reside
in tensor.h.

include/deal.II/base/tensor.h
include/deal.II/base/tensor_base.h

index 32eaae6b8b70f7c4335e659767d634c051688283..224e3fc5fdd139369b027f3ad473a527a7a37c98 100644 (file)
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/tensor_base.h>
-#include <deal.II/base/utilities.h>
 
 DEAL_II_NAMESPACE_OPEN
-template <int rank_, int dim, typename Number> class Tensor;
-template <int dim, typename Number> class Tensor<1,dim,Number>;
-
-/**
- * A general tensor class with an arbitrary rank, i.e. with an arbitrary
- * number of indices. The Tensor class provides an indexing operator and a bit
- * of infrastructure, but most functionality is recursively handed down to
- * tensors of rank 1 or put into external templated functions, e.g. the
- * <tt>contract</tt> family.
- *
- * 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. It also
- * makes the code easier to read because of the semantic difference between a
- * tensor (an object that relates to a coordinate system and has
- * transformation properties with regard to coordinate rotations and
- * transforms) and matrices (which we consider as operators on arbitrary
- * vector spaces related to linear algebra things).
- *
- * @tparam rank_ An integer that denotes the rank of this tensor. A rank-0
- * tensor is a scalar, a rank-1 tensor is a vector with @p dim components, a
- * rank-2 tensor is a matrix with dim-by-dim components, etc. There are
- * specializations of this class for rank-0 and rank-1 tensors. There is also
- * a related class SymmetricTensor for tensors of even rank whose elements are
- * symmetric.
- * @tparam dim An integer that denotes the dimension of the space in which
- * this tensor operates. This of course equals the number of coordinates that
- * identify a point and rank-1 tensor.
- * @tparam Number The data type in which the tensor elements are to be stored.
- * This will, in almost all cases, simply be the default @p double, but there
- * are cases where one may want to store elements in a different (and always
- * scalar) type. It can be used to base tensors on @p float or @p complex
- * numbers or any other data type that implements basic arithmetic operations.
- * Another example would be a type that allows for Automatic Differentiation
- * (see, for example, the Sacado type used in step-33) and thereby can
- * generate analytic (spatial) derivatives of a function that takes a tensor
- * as argument.
- *
- * @ingroup geomprimitives
- * @author Wolfgang Bangerth, 1998-2005
- */
-template <int rank_, int dim, typename Number>
-class Tensor
-{
-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 = rank_;
-
-  /**
-   * Number of independent components of a tensor of current rank. This is dim
-   * times the number of independent components of each sub-tensor.
-   */
-  static const unsigned int
-  n_independent_components = Tensor<rank_-1,dim>::n_independent_components *dim;
-
-  /**
-   * Type of stored objects. This is a tensor of lower rank.
-   */
-  typedef Tensor<rank_-1,dim,Number> value_type;
-
-  /**
-   * Declare a type that has holds real-valued numbers with the same precision
-   * as the template argument to this class. For std::complex<number>, this
-   * corresponds to type number, and it is equal to Number for all other
-   * cases. See also the respective field in Vector<Number>.
-   *
-   * This typedef is used to represent the return type of norms.
-   */
-  typedef typename numbers::NumberTraits<Number>::real_type real_type;
-
-  /**
-   * Declare an array type which can be used to initialize an object of this
-   * type statically.
-   */
-  typedef typename Tensor<rank_-1,dim,Number>::array_type array_type[dim];
-
-  /**
-   * Constructor. Initialize all entries to zero if
-   * <tt>initialize==true</tt>; this is the default behaviour.
-   */
-  explicit
-  Tensor (const bool initialize = true);
-
-  /**
-   * Copy constructor.
-   */
-  Tensor (const Tensor<rank_,dim,Number> &initializer);
-
-  /**
-   * Constructor, where the data is copied from a C-style array.
-   */
-  Tensor (const array_type &initializer);
-
-  /**
-   * Constructor from tensors with different underlying scalar type. This
-   * obviously requires that the @p OtherNumber type is convertible to @p
-   * Number.
-   */
-  template <typename OtherNumber>
-  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
-
-  /**
-   * Constructor that converts from a "tensor of tensors".
-   */
-  template <typename OtherNumber>
-  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
-
-  /**
-   * Conversion operator to tensor of tensors.
-   */
-  template <typename OtherNumber>
-  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
-
-  /**
-   * Read-Write access operator.
-   */
-  Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i);
-
-  /**
-   * Read-only access operator.
-   */
-  const Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i) const;
-
-  /**
-   * Read access using TableIndices <tt>indices</tt>
-   */
-  Number operator [] (const TableIndices<rank_> &indices) const;
-
-  /**
-   * Read and write access using TableIndices <tt>indices</tt>
-   */
-  Number &operator [] (const TableIndices<rank_> &indices);
-
-  /**
-   * Copy assignment operator.
-   */
-  Tensor &operator = (const Tensor<rank_,dim,Number> &rhs);
-
-  /**
-   * Assignment operator from tensors with different underlying scalar type.
-   * This obviously requires that the @p OtherNumber type is convertible to @p
-   * Number.
-   */
-  template <typename OtherNumber>
-  Tensor &operator = (const Tensor<rank_,dim,OtherNumber> &rhs);
-
-  /**
-   * This operator assigns a scalar to a tensor. To avoid confusion with what
-   * exactly it means to assign a scalar value to a tensor, zero is the only
-   * value allowed for <tt>d</tt>, allowing the intuitive notation
-   * <tt>t=0</tt> to reset all elements of the tensor to zero.
-   *
-   * @relates EnableIfScalar
-   */
-  template <typename OtherNumber,
-            typename = typename EnableIfScalar<OtherNumber>::type>
-  Tensor<rank_,dim,Number> &operator = (const OtherNumber d);
-
-  /**
-   * Test for equality of two tensors.
-   */
-  template <typename OtherNumber>
-  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
-
-  /**
-   * Test for inequality of two tensors.
-   */
-  template <typename OtherNumber>
-  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
-
-  /**
-   * Add another tensor.
-   */
-  template <typename OtherNumber>
-  Tensor<rank_,dim,Number> &operator += (const Tensor<rank_,dim,OtherNumber> &);
-
-  /**
-   * Subtract another tensor.
-   */
-  template <typename OtherNumber>
-  Tensor<rank_,dim,Number> &operator -= (const Tensor<rank_,dim,OtherNumber> &);
-
-  /**
-   * Scale the tensor by <tt>factor</tt>, i.e. multiply all components by
-   * <tt>factor</tt>.
-   */
-  template <typename OtherNumber>
-  Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
-
-  /**
-   * Scale the vector by <tt>1/factor</tt>.
-   */
-  template <typename OtherNumber>
-  Tensor<rank_,dim,Number> &operator /= (const OtherNumber factor);
-
-  /**
-   * Unary minus operator. Negate all entries of a tensor.
-   */
-  Tensor<rank_,dim,Number>   operator - () const;
-
-  /**
-   * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
-   * squares of all entries.
-   */
-  real_type norm () const;
-
-  /**
-   * Return the square of the Frobenius-norm of a tensor, i.e. the sum of
-   * squares of all entries.
-   *
-   * This function mainly exists because it makes computing the norm simpler
-   * recursively, but may also be useful in other contexts.
-   */
-  real_type norm_square () const;
-
-  /**
-   * Fill a vector with all tensor elements.
-   *
-   * This function unrolls all tensor entries into a single, linearly numbered
-   * vector. As usual in C++, the rightmost index of the tensor marches
-   * fastest.
-   */
-  template <typename OtherNumber>
-  void unroll (Vector<OtherNumber> &result) const;
-
-  /**
-   * Returns an unrolled index in the range [0,dim^rank-1] for the element of
-   * the tensor indexed by the argument to the function.
-   */
-  static
-  unsigned int
-  component_to_unrolled_index(const TableIndices<rank_> &indices);
-
-  /**
-   * Opposite of  component_to_unrolled_index: For an index in the range
-   * [0,dim^rank-1], return which set of indices it would correspond to.
-   */
-  static
-  TableIndices<rank_> unrolled_to_component_indices(const unsigned int i);
-
-  /**
-   * Reset all values to zero.
-   *
-   * Note that this is partly inconsistent with the semantics of the @p
-   * clear() member functions of the standard library containers 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 std::size_t memory_consumption ();
-
-  /**
-   * Exception.
-   */
-  DeclException1 (ExcInvalidTensorContractionIndex,
-                  int,
-                  << "You have requested contraction of tensors over index "
-                  << arg1
-                  << ", but this is not possible for tensors of the current type.");
-
-  /**
-   * Read or write the data of this object to or from a stream for the purpose
-   * of serialization
-   */
-  template <class Archive>
-  void serialize(Archive &ar, const unsigned int version);
-
-private:
-  /**
-   * Array of tensors holding the subelements.
-   */
-  Tensor<rank_-1,dim,Number> subtensor[dim];
-
-  /**
-   * Help function for unroll.
-   */
-  template <typename OtherNumber>
-  void unroll_recursion(Vector<OtherNumber> &result,
-                        unsigned int    &start_index) const;
-
-  // make the following class a
-  // friend to this class. in principle,
-  // it would suffice if otherrank==rank+1,
-  // but then the compiler complains
-  // that this be an explicit specialization
-  // which is not what we want
-  //
-  // also, it would be sufficient to make
-  // the function unroll_loops a friend,
-  // but that seems to be impossible as well.
-  template <int, int, typename> friend class Tensor;
-};
-
-
-/*--------------------------- Inline functions -----------------------------*/
-
-#ifndef DOXYGEN
-
-template <int rank_, int dim, typename Number>
-inline
-Tensor<rank_,dim,Number>::Tensor (const bool initialize)
-{
-  if (initialize)
-    // need to create an object Number() to initialize to zero to avoid
-    // confusion with Tensor::operator=(scalar) when using something like
-    // Tensor<1,dim,Tensor<1,dim,Number> >.
-    for (unsigned int i=0; i!=dim; ++i)
-      subtensor[i] = Tensor<rank_-1,dim,Number>();
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Tensor<rank_,dim,Number>::Tensor (const Tensor<rank_,dim,Number> &initializer)
-{
-  for (unsigned int i=0; i!=dim; ++i)
-    subtensor[i] = initializer[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] = initializer[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number>::Tensor (const Tensor<rank_,dim,OtherNumber> &initializer)
-{
-  for (unsigned int i=0; i!=dim; ++i)
-    subtensor[i] = initializer[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number>::Tensor
-(const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] = initializer[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number>::operator
-Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
-{
-  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (subtensor);
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-typename Tensor<rank_,dim,Number>::value_type &
-Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
-{
-  Assert (i<dim, ExcIndexRange(i, 0, dim));
-  return subtensor[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-const typename Tensor<rank_,dim,Number>::value_type &
-Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
-{
-  Assert (i<dim, ExcIndexRange(i, 0, dim));
-
-  return subtensor[i];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Number
-Tensor<rank_,dim,Number>::operator[] (const TableIndices<rank_> &indices) const
-{
-  const unsigned int inner_ind = indices[0];
-  Assert (inner_ind<dim, ExcIndexRange(inner_ind, 0, dim));
-
-  TableIndices<rank_-1> indices1;
-  for (unsigned int i = 0; i < rank_-1; i++)
-    indices1[i] = indices[i+1];
-  return (subtensor[inner_ind])[indices1];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Number &
-Tensor<rank_,dim,Number>::operator[] (const TableIndices<rank_> &indices)
-{
-  const unsigned int inner_ind = indices[0];
-  Assert (inner_ind<dim, ExcIndexRange(inner_ind, 0, dim));
-
-  TableIndices<rank_-1> indices1;
-  for (unsigned int i = 0; i < rank_-1; i++)
-    indices1[i] = indices[i+1];
-  return (subtensor[inner_ind])[indices1];
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,Number> &t)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] = t.subtensor[i];
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,OtherNumber> &t)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] = t.subtensor[i];
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber, typename>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator = (const OtherNumber d)
-{
-  Assert (d == OtherNumber(), ExcMessage ("Only assignment with zero is allowed"));
-  (void) d;
-
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] = Number();
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-bool
-Tensor<rank_,dim,Number>::operator == (const Tensor<rank_,dim,OtherNumber> &p) const
-{
-  for (unsigned int i=0; i<dim; ++i)
-    if (subtensor[i] != p.subtensor[i])
-      return false;
-  return true;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-bool
-Tensor<rank_,dim,Number>::operator != (const Tensor<rank_,dim,OtherNumber> &p) const
-{
-  return !((*this) == p);
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator += (const Tensor<rank_,dim,OtherNumber> &p)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] += p.subtensor[i];
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator -= (const Tensor<rank_,dim,OtherNumber> &p)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] -= p.subtensor[i];
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator *= (const OtherNumber s)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] *= s;
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-Tensor<rank_,dim,Number> &
-Tensor<rank_,dim,Number>::operator /= (const OtherNumber s)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] /= s;
-  return *this;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-Tensor<rank_,dim,Number>
-Tensor<rank_,dim,Number>::operator - () const
-{
-  Tensor<rank_,dim,Number> tmp;
-
-  for (unsigned int i=0; i<dim; ++i)
-    tmp.subtensor[i] = -subtensor[i];
-
-  return tmp;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-typename Tensor<rank_,dim,Number>::real_type
-Tensor<rank_,dim,Number>::norm () const
-{
-  return std::sqrt (norm_square());
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-typename Tensor<rank_,dim,Number>::real_type
-Tensor<rank_,dim,Number>::norm_square () const
-{
-  real_type s = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    s += subtensor[i].norm_square();
-
-  return s;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-void
-Tensor<rank_, dim, Number>::unroll (Vector<OtherNumber> &result) const
-{
-  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
-
-  unsigned int index = 0;
-  unroll_recursion (result, index);
-}
 
 
 
-template <int rank_, int dim, typename Number>
-template <typename OtherNumber>
-inline
-void
-Tensor<rank_, dim, Number>::unroll_recursion (Vector<OtherNumber> &result,
-                                              unsigned int    &index) const
-{
-  for (unsigned int i=0; i<dim; ++i)
-    {
-      operator[](i).unroll_recursion(result, index);
-    }
-}
-
-template <int rank_, int dim, typename Number>
-inline
-unsigned int
-Tensor<rank_, dim, Number>::component_to_unrolled_index(const TableIndices<rank_> &indices)
-{
-  TableIndices<rank_-1> indices1;
-  for (unsigned int i = 0; i < rank_-1; i++)
-    indices1[i] = indices[i];
-
-  Assert (indices[rank_-1] < dim,
-          ExcIndexRange (indices[rank_-1], 0, dim));
-  return ( Tensor<rank_-1,dim,Number>::component_to_unrolled_index(indices1) * dim + indices[rank_-1]);
-}
-
-template <int rank_, int dim, typename Number>
-inline
-TableIndices<rank_>
-Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
-{
-  Assert (i < n_independent_components,
-          ExcIndexRange (i, 0, n_independent_components));
-
-  TableIndices<rank_>   indices;
-
-  unsigned int remainder = i;
-  for (int r=rank_-1; r>=0; --r)
-    {
-      indices[r] = (remainder % dim);
-      remainder /= dim;
-    }
-  Assert (remainder == 0, ExcInternalError());
-
-  return indices;
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-void Tensor<rank_,dim,Number>::clear ()
-{
-  for (unsigned int i=0; i<dim; ++i)
-    subtensor[i].clear();
-}
-
-
-
-template <int rank_, int dim, typename Number>
-inline
-std::size_t
-Tensor<rank_,dim,Number>::memory_consumption ()
-{
-  return sizeof(Tensor<rank_,dim,Number>);
-}
-
-
-
-template <int rank_, int dim, typename Number>
-template <class Archive>
-inline
-void
-Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
-{
-  ar &subtensor;
-}
-
-#endif // DOXYGEN
 /* ----------------- Non-member functions operating on tensors. ------------ */
 
 
index 75735b46f3e2ba1896a37c398f00b29c20aca8ae..ada55ccbd85e3351504debb8325b7b2ca6818c15 100644 (file)
@@ -27,6 +27,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/table_indices.h>
 #include <deal.II/base/template_constraints.h>
+#include <deal.II/base/utilities.h>
 #include <vector>
 
 #include <cmath>
@@ -43,11 +44,9 @@ template <typename number> class Vector;
 // specified)
 template <int dim, typename Number> class Point;
 
-// general template; specialized for rank==1; the general template is in
-// tensor.h
-template <int rank_, int dim, typename Number=double> class Tensor;
+// general template; specialized for rank == 0
+template <int rank_, int dim, typename Number = double> class Tensor;
 template <int dim, typename Number> class Tensor<0,dim,Number>;
-template <int dim, typename Number> class Tensor<1,dim,Number>;
 
 
 
@@ -98,13 +97,12 @@ public:
   /**
    * Publish the rank of this tensor to the outside world.
    */
-  static const unsigned int rank      = 0;
+  static const unsigned int rank = 0;
 
   /**
-   * Type of stored objects. This is a Number for a rank 1 tensor.
+   * Number of independent components of a tensor of rank 0.
    */
-
-  typedef Number value_type;
+  static const unsigned int n_independent_components = 1;
 
   /**
    * Declare a type that has holds real-valued numbers with the same precision
@@ -116,10 +114,22 @@ public:
    */
   typedef typename numbers::NumberTraits<Number>::real_type real_type;
 
+  /**
+   * The tensor type this object represents. In the special case of a
+   * tensor or rank 0 we strip the tensor class and just store the scalar
+   * object.
+   */
+  typedef Number tensor_type;
+
+  /**
+   * Type of stored objects. This is a Number for a rank 0 tensor.
+   */
+  typedef Number value_type;
+
   /**
    * Declare an array type which can be used to initialize an object of this
-   * type statically. In case of a a tensor of rank 0 this is just a scalar
-   * number type
+   * type statically. In case of a a tensor of rank 0 this is just the scalar
+   * number type Number.
    */
   typedef Number array_type;
 
@@ -222,6 +232,20 @@ public:
    */
   Tensor<0,dim,Number>   operator - () const;
 
+  /**
+   * Reset all values to zero.
+   *
+   * Note that this is partly inconsistent with the semantics of the @p
+   * clear() member functions of the standard library containers 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 ();
+
   /**
    * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
    * squares of all entries. For the present case of rank-1 tensors, this
@@ -239,18 +263,11 @@ public:
   real_type norm_square () const;
 
   /**
-   * Reset all values to zero.
-   *
-   * Note that this is partly inconsistent with the semantics of the @p
-   * clear() member functions of the standard library containers 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.
+   * Read or write the data of this object to or from a stream for the purpose
+   * of serialization
    */
-  void clear ();
+  template <class Archive>
+  void serialize(Archive &ar, const unsigned int version);
 
   /**
    * Only tensors with a positive dimension are implemented. This exception is
@@ -263,44 +280,55 @@ public:
                   int,
                   << "dim must be positive, but was " << arg1);
 
-  /**
-   * Read or write the data of this object to or from a stream for the purpose
-   * of serialization
-   */
-  template <class Archive>
-  void serialize(Archive &ar, const unsigned int version);
-
 private:
   /**
    * The value of this scalar object.
    */
   Number value;
+
+  /**
+   * Help function for unroll.
+   */
+  template <typename OtherNumber>
+  void unroll_recursion(Vector<OtherNumber> &result,
+                        unsigned int        &start_index) const;
+
+  /**
+   * Allow an arbitrary Tensor to access the underlying values.
+   */
+  template <int, int, typename> friend class Tensor;
 };
 
 
 
 /**
- * This class is a specialized version of the <tt>Tensor<rank,dim,Number></tt>
- * class.  It handles tensors with one index, i.e. vectors, of fixed dimension
- * and provides the basis for the functionality needed for tensors of higher
- * rank.
+ * A general tensor class with an arbitrary rank, i.e. with an arbitrary
+ * number of indices. The Tensor class provides an indexing operator and a bit
+ * of infrastructure, but most functionality is recursively handed down to
+ * tensors of rank 1 or put into external templated functions, e.g. the
+ * <tt>contract</tt> family.
  *
- * Within deal.II, the distinction between this class and its derived class
- * <tt>Point</tt> is that we use the <tt>Point</tt> class mainly to denote the
- * points that make up geometric objects. As such, they have a small number of
- * additional operations over general tensors of rank 1 for which we use the
- * <tt>Tensor<1,dim,Number></tt> class. In particular, there is a distance()
- * function to compute the Euclidean distance between two points in space.
+ * 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. It also
+ * makes the code easier to read because of the semantic difference between a
+ * tensor (an object that relates to a coordinate system and has
+ * transformation properties with regard to coordinate rotations and
+ * transforms) and matrices (which we consider as operators on arbitrary
+ * vector spaces related to linear algebra things).
  *
- * However, the <tt>Point</tt> class is really only used where the coordinates
- * of an object can be thought to possess the dimension of a length. For all
- * other uses, such as the gradient of a scalar function (which is a tensor of
- * rank 1, or vector, with as many elements as a point object, but with
- * different physical units), we use the <tt>Tensor<1,dim,Number></tt> class.
+ * @tparam rank_ An integer that denotes the rank of this tensor. A rank-0
+ * tensor is a scalar, a rank-1 tensor is a vector with @p dim components, a
+ * rank-2 tensor is a matrix with dim-by-dim components, etc. There are
+ * specializations of this class for rank-0 and rank-1 tensors. There is also
+ * a related class SymmetricTensor for tensors of even rank whose elements are
+ * symmetric.
  *
  * @tparam dim An integer that denotes the dimension of the space in which
  * this tensor operates. This of course equals the number of coordinates that
  * identify a point and rank-1 tensor.
+ *
  * @tparam Number The data type in which the tensor elements are to be stored.
  * This will, in almost all cases, simply be the default @p double, but there
  * are cases where one may want to store elements in a different (and always
@@ -314,8 +342,8 @@ private:
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 1998-2005, Matthias Maier, 2015
  */
-template <int dim,typename Number>
-class Tensor<1,dim,Number>
+template <int rank_, int dim, typename Number>
+class Tensor
 {
 public:
   /**
@@ -331,22 +359,17 @@ public:
   /**
    * Publish the rank of this tensor to the outside world.
    */
-  static const unsigned int rank = 1;
+  static const unsigned int rank = rank_;
 
   /**
-   * Number of independent components of a tensor of rank 1.
+   * Number of independent components of a tensor of current rank. This is dim
+   * times the number of independent components of each sub-tensor.
    */
   static const unsigned int
-  n_independent_components = dim;
-
-  /**
-   * Type of stored objects. This is a Number for a rank 1 tensor.
-   */
-
-  typedef Number value_type;
+  n_independent_components = Tensor<rank_-1,dim>::n_independent_components * dim;
 
   /**
-   * Declare a type that has holds real-valued numbers with the same precision
+   * Declare a type that holds real-valued numbers with the same precision
    * as the template argument to this class. For std::complex<number>, this
    * corresponds to type number, and it is equal to Number for all other
    * cases. See also the respective field in Vector<Number>.
@@ -356,17 +379,26 @@ public:
   typedef typename numbers::NumberTraits<Number>::real_type real_type;
 
   /**
-   * Declare an array type which can be used to initialize statically an
-   * object of this type.
+   * The tensor type this object represents.
+   */
+  typedef Tensor<rank_,dim,Number> tensor_type;
+
+  /**
+   * Type of stored objects (i.e., the object returned by operator[]()).
+   * This is a tensor of lower rank.
    */
-  // Avoid a bogus warning in case of dim==0, and always provide a type
-  // with positive array size. The constructor will take care that no
-  // Tensor with dim==0 will be constructed.
-  typedef Number array_type[(dim!=0) ? dim : 1];
+  typedef typename Tensor<rank_-1,dim,Number>::tensor_type value_type;
 
   /**
-   * Constructor. Initialize all entries to zero if <tt>initialize==true</tt>;
-   * this is the default behaviour.
+   * Declare an array type which can be used to initialize an object of this
+   * type statically.
+   */
+  typedef typename Tensor<rank_-1,dim,Number>::array_type
+    array_type[(dim != 0) ? dim : 1];
+
+  /**
+   * Constructor. Initialize all entries to zero if
+   * <tt>initialize==true</tt>; this is the default behaviour.
    */
   explicit
   Tensor (const bool initialize = true);
@@ -374,51 +406,57 @@ public:
   /**
    * Copy constructor.
    */
-  Tensor (const Tensor<1,dim,Number> &initializer);
+  Tensor (const Tensor<rank_,dim,Number> &initializer);
 
   /**
-   * Copy constructor, where the data is copied from a C-style array.
+   * Constructor, where the data is copied from a C-style array.
    */
   Tensor (const array_type &initializer);
 
   /**
-   * Copy constructor from tensors with different underlying scalar type. This
+   * Constructor from tensors with different underlying scalar type. This
    * obviously requires that the @p OtherNumber type is convertible to @p
    * Number.
    */
   template <typename OtherNumber>
-  Tensor (const Tensor<1,dim,OtherNumber> &initializer);
+  Tensor (const Tensor<rank_,dim,OtherNumber> &initializer);
 
   /**
-   * Read access to the <tt>index</tt>th coordinate.
-   *
-   * Note that the derived <tt>Point</tt> class also provides access through
-   * the <tt>()</tt> operator for backcompatibility.
+   * Constructor that converts from a "tensor of tensors".
+   */
+  template <typename OtherNumber>
+  Tensor (const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer);
+
+  /**
+   * Conversion operator to tensor of tensors.
+   */
+  template <typename OtherNumber>
+  operator Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const;
+
+  /**
+   * Read-Write access operator.
    */
-  Number operator [] (const unsigned int index) const;
+  value_type & operator [] (const unsigned int i);
 
   /**
-   * Read and write access to the <tt>index</tt>th coordinate.
-   *
-   * Note that the derived <tt>Point</tt> class also provides access through
-   * the <tt>()</tt> operator for backcompatibility.
+   * Read-only access operator.
    */
-  Number &operator [] (const unsigned int index);
+  const value_type & operator[](const unsigned int i) const;
 
   /**
    * Read access using TableIndices <tt>indices</tt>
    */
-  Number operator [] (const TableIndices<1> &indices) const;
+  Number operator [] (const TableIndices<rank_> &indices) const;
 
   /**
    * Read and write access using TableIndices <tt>indices</tt>
    */
-  Number &operator [] (const TableIndices<1> &indices);
+  Number &operator [] (const TableIndices<rank_> &indices);
 
   /**
    * Copy assignment operator.
    */
-  Tensor<1,dim,Number> &operator = (const Tensor<1,dim,Number> &rhs);
+  Tensor &operator = (const Tensor<rank_,dim,Number> &rhs);
 
   /**
    * Assignment operator from tensors with different underlying scalar type.
@@ -426,74 +464,61 @@ public:
    * Number.
    */
   template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator = (const Tensor<1,dim,OtherNumber> &rhs);
+  Tensor &operator = (const Tensor<rank_,dim,OtherNumber> &rhs);
 
   /**
    * This operator assigns a scalar to a tensor. To avoid confusion with what
    * exactly it means to assign a scalar value to a tensor, zero is the only
    * value allowed for <tt>d</tt>, allowing the intuitive notation
    * <tt>t=0</tt> to reset all elements of the tensor to zero.
+   *
+   * @relates EnableIfScalar
    */
-  template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator = (const OtherNumber d);
+  template <typename OtherNumber,
+            typename = typename EnableIfScalar<OtherNumber>::type>
+  Tensor<rank_,dim,Number> &operator = (const OtherNumber d);
 
   /**
    * Test for equality of two tensors.
    */
   template <typename OtherNumber>
-  bool operator == (const Tensor<1,dim,OtherNumber> &rhs) const;
+  bool operator == (const Tensor<rank_,dim,OtherNumber> &) const;
 
   /**
    * Test for inequality of two tensors.
    */
   template <typename OtherNumber>
-  bool operator != (const Tensor<1,dim,OtherNumber> &rhs) const;
+  bool operator != (const Tensor<rank_,dim,OtherNumber> &) const;
 
   /**
-   * Add another vector to this vector.
+   * Add another tensor.
    */
   template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator += (const Tensor<1,dim,OtherNumber> &rhs);
+  Tensor<rank_,dim,Number> &operator += (const Tensor<rank_,dim,OtherNumber> &);
 
   /**
-   * Subtract another vector.
+   * Subtract another tensor.
    */
   template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator -= (const Tensor<1,dim,OtherNumber> &rhs);
+  Tensor<rank_,dim,Number> &operator -= (const Tensor<rank_,dim,OtherNumber> &);
 
   /**
-   * Scale the vector by <tt>factor</tt>, i.e., multiply all coordinates by
+   * Scale the tensor by <tt>factor</tt>, i.e. multiply all components by
    * <tt>factor</tt>.
    */
   template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator *= (const OtherNumber factor);
+  Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
 
   /**
    * Scale the vector by <tt>1/factor</tt>.
    */
   template <typename OtherNumber>
-  Tensor<1,dim,Number> &operator /= (const OtherNumber factor);
-
-  /**
-   * Tensor with inverted entries.
-   */
-  Tensor<1,dim,Number> operator - () const;
-
-  /**
-   * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
-   * squares of all entries. For the present case of rank-1 tensors, this
-   * equals the usual <tt>l<sub>2</sub></tt> norm of the vector.
-   */
-  real_type norm () const;
+  Tensor<rank_,dim,Number> &operator /= (const OtherNumber factor);
 
   /**
-   * Return the square of the Frobenius-norm of a tensor, i.e. the square root
-   * of the sum of squares of all entries.
-   *
-   * This function mainly exists because it makes computing the norm simpler
-   * recursively, but may also be useful in other contexts.
+   * Unary minus operator. Negate all entries of a tensor.
    */
-  real_type norm_square () const;
+  Tensor<rank_,dim,Number>   operator - () const;
 
   /**
    * Reset all values to zero.
@@ -509,37 +534,45 @@ public:
    */
   void clear ();
 
+  /**
+   * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
+   * squares of all entries.
+   */
+  real_type norm () const;
+
+  /**
+   * Return the square of the Frobenius-norm of a tensor, i.e. the sum of
+   * squares of all entries.
+   *
+   * This function mainly exists because it makes computing the norm simpler
+   * recursively, but may also be useful in other contexts.
+   */
+  real_type norm_square () const;
+
   /**
    * Fill a vector with all tensor elements.
    *
    * This function unrolls all tensor entries into a single, linearly numbered
-   * vector. As usual in C++, the rightmost index marches fastest.
+   * vector. As usual in C++, the rightmost index of the tensor marches
+   * fastest.
    */
-  template <typename Number2>
-  void unroll (Vector<Number2> &result) const;
+  template <typename OtherNumber>
+  void unroll (Vector<OtherNumber> &result) const;
 
   /**
-   * Returns an unrolled index in the range [0,dim-1] for the element of the
-   * tensor indexed by the argument to the function.
-   *
-   * Given that this is a rank-1 object, the returned value is simply the
-   * value of the only index stored by the argument.
+   * Returns an unrolled index in the range [0,dim^rank-1] for the element of
+   * the tensor indexed by the argument to the function.
    */
   static
   unsigned int
-  component_to_unrolled_index(const TableIndices<1> &indices);
+  component_to_unrolled_index(const TableIndices<rank_> &indices);
 
   /**
    * Opposite of  component_to_unrolled_index: For an index in the range
-   * [0,dim-1], return which set of indices it would correspond to.
-   *
-   * Given that this is a rank-1 object, the returned set of indices consists
-   * of only a single element with value equal to the argument to this
-   * function.
+   * [0,dim^rank-1], return which set of indices it would correspond to.
    */
   static
-  TableIndices<1> unrolled_to_component_indices(const unsigned int i);
-
+  TableIndices<rank_> unrolled_to_component_indices(const unsigned int i);
 
   /**
    * Determine an estimate for the memory consumption (in bytes) of this
@@ -547,6 +580,22 @@ public:
    */
   static std::size_t memory_consumption ();
 
+  /**
+   * Read or write the data of this object to or from a stream for the purpose
+   * of serialization
+   */
+  template <class Archive>
+  void serialize(Archive &ar, const unsigned int version);
+
+  /**
+   * Exception.
+   */
+  DeclException1 (ExcInvalidTensorContractionIndex,
+                  int,
+                  << "You have requested contraction of tensors over index "
+                  << arg1
+                  << ", but this is not possible for tensors of the current type.");
+
   /**
    * Only tensors with a positive dimension are implemented. This exception is
    * thrown by the constructor if the template argument <tt>dim</tt> is zero
@@ -558,39 +607,23 @@ public:
                   int,
                   << "dim must be positive, but was " << arg1);
 
-  /**
-   * Read or write the data of this object to or from a stream for the purpose
-   * of serialization
-   */
-  template <class Archive>
-  void serialize(Archive &ar, const unsigned int version);
-
 private:
   /**
-   * Store the values in a simple array.  For <tt>dim==0</tt> store one
-   * element, because otherwise the compiler would choke.  We catch this case
-   * in the constructor to disallow the creation of such an object.
+   * Array of tensors holding the subelements.
    */
-  array_type values;
+  value_type values[(dim != 0) ? dim : 1];
 
   /**
-   * Help function for unroll. If we have detected an access control bug in
-   * the compiler, this function is declared public, otherwise private. Do not
-   * attempt to use this function from outside in any case, even if it should
-   * be public for your compiler.
+   * Help function for unroll.
    */
-  template <typename Number2>
-  void unroll_recursion (Vector<Number2> &result,
-                         unsigned int    &start_index) const;
+  template <typename OtherNumber>
+  void unroll_recursion(Vector<OtherNumber> &result,
+                        unsigned int        &start_index) const;
 
   /**
-   * Make the following classes friends to this class. In principle, it would
-   * suffice if otherrank==2, but that is not possible in C++ at present.
-   *
-   * Also, it would be sufficient to make the function unroll_loops a friend,
-   * but that seems to be impossible as well.
+   * Allow an arbitrary Tensor to access the underlying values.
    */
-  template <int otherrank, int otherdim, typename OtherNumber>  friend class dealii::Tensor;
+  template <int, int, typename> friend class Tensor;
 
   /**
    * Point is allowed access to the coordinates. This is supposed to improve
@@ -600,23 +633,8 @@ private:
 };
 
 
-/**
- * Prints the value of this scalar.
- */
-template <int dim,typename Number>
-std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p);
-
-/**
- * Prints the values of this tensor in the form <tt>x1 x2 x3 etc</tt>.
- */
-template <int dim,typename Number>
-std::ostream &operator << (std::ostream &out, const Tensor<1,dim,Number> &p);
-
 
 #ifndef DOXYGEN
-
-
-
 /*---------------------- Inline functions: Tensor<0,dim> ---------------------*/
 
 
@@ -631,7 +649,6 @@ Tensor<0,dim,Number>::Tensor ()
 }
 
 
-
 template <int dim, typename Number>
 inline
 Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p)
@@ -642,7 +659,6 @@ Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p)
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -654,7 +670,6 @@ Tensor<0,dim,Number>::Tensor (const OtherNumber initializer)
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -666,7 +681,6 @@ Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,OtherNumber> &p)
 }
 
 
-
 template <int dim, typename Number>
 inline
 Tensor<0,dim,Number>::operator Number () const
@@ -675,7 +689,6 @@ Tensor<0,dim,Number>::operator Number () const
 }
 
 
-
 template <int dim, typename Number>
 inline
 Tensor<0,dim,Number>::operator Number &()
@@ -684,7 +697,6 @@ Tensor<0,dim,Number>::operator Number &()
 }
 
 
-
 template <int dim, typename Number>
 inline
 Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p)
@@ -694,7 +706,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Numbe
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -705,7 +716,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Other
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -716,7 +726,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const OtherNumber d)
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -726,7 +735,6 @@ bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,OtherNumber> &p) cons
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -736,7 +744,6 @@ bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,OtherNumber> &p) cons
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -747,7 +754,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator += (const Tensor<0,dim,Othe
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -758,7 +764,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Othe
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -769,7 +774,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s)
 }
 
 
-
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline
@@ -780,7 +784,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator /= (const OtherNumber s)
 }
 
 
-
 template <int dim, typename Number>
 inline
 Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const
@@ -789,7 +792,6 @@ Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const
 }
 
 
-
 template <int dim, typename Number>
 inline
 typename Tensor<0,dim,Number>::real_type
@@ -799,7 +801,6 @@ Tensor<0,dim,Number>::norm () const
 }
 
 
-
 template <int dim, typename Number>
 inline
 typename Tensor<0,dim,Number>::real_type
@@ -809,6 +810,17 @@ Tensor<0,dim,Number>::norm_square () const
 }
 
 
+template <int dim, typename Number>
+template <typename OtherNumber>
+inline
+void
+Tensor<0, dim, Number>::unroll_recursion (Vector<OtherNumber> &result,
+                                          unsigned int        &index) const
+{
+  result[index] = value;
+  ++index;
+}
+
 
 template <int dim, typename Number>
 inline
@@ -818,7 +830,6 @@ void Tensor<0,dim,Number>::clear ()
 }
 
 
-
 template <int dim, typename Number>
 template <class Archive>
 inline
@@ -828,151 +839,58 @@ void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
 }
 
 
-
-#ifndef DEAL_II_WITH_CXX11
-
-template <typename T, typename U, int rank, int dim>
-struct ProductType<T,Tensor<rank,dim,U> >
-{
-  typedef Tensor<rank,dim,typename ProductType<T,U>::type> type;
-};
-
-template <typename T, typename U, int rank, int dim>
-struct ProductType<Tensor<rank,dim,T>,U>
-{
-  typedef Tensor<rank,dim,typename ProductType<T,U>::type> type;
-};
-
-#endif
-
-
-
-/**
- * TODO
- *
- * @relates Tensor
- * @relates EnableIfScalar
- */
-template <int dim,
-         typename Number,
-         typename OtherNumber,
-         typename = typename EnableIfScalar<OtherNumber>::type>
-inline
-Tensor<0,dim,typename ProductType<OtherNumber, Number>::type>
-operator * (const OtherNumber           factor,
-            const Tensor<0,dim,Number> &t)
-{
-  return factor * static_cast<Number>(t);
-}
-
-
-
-/**
- * TODO
- *
- * @relates Tensor
- * @relates EnableIfScalar
- */
-template <int dim,
-         typename Number,
-         typename OtherNumber,
-         typename = typename EnableIfScalar<OtherNumber>::type>
-inline
-Tensor<0,dim,typename ProductType<Number, OtherNumber>::type>
-operator * (const Tensor<0,dim,Number> &t,
-            const OtherNumber           factor)
-{
-  return static_cast<Number>(t) * factor;
-}
-
-
-
-/**
- * TODO
- *
- * @relates Tensor
- * @relates EnableIfScalar
- */
-template <int dim,
-         typename Number,
-         typename OtherNumber,
-         typename = typename EnableIfScalar<OtherNumber>::type>
-inline
-Tensor<0,dim,typename ProductType<Number, OtherNumber>::type>
-operator / (const Tensor<0,dim,Number> &t,
-            const OtherNumber           factor)
-{
-  return static_cast<Number>(t) / factor;
-}
-
-
-
-/**
- * Add two tensors of rank 0.
- *
- * @relates Tensor
- */
-template <int dim, typename Number, typename OtherNumber>
-inline
-Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
-operator+ (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
-{
-  return static_cast<Number>(p) + static_cast<OtherNumber>(q);
-}
-
-
-
-/**
- * Subtract two tensors of rank 0.
- *
- * @relates Tensor
- */
-template <int dim, typename Number, typename OtherNumber>
-inline
-Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
-operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
-{
-  return static_cast<Number>(p) - static_cast<OtherNumber>(q);
-}
-
+/*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
 
 
-/**
- * Returns the contraction of two Tensors of rank 0.
- *
- * @relates Tensor
- */
-template <int dim, typename Number, typename OtherNumber>
-inline
-typename ProductType<Number, OtherNumber>::type
-operator* (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
+namespace internal
 {
-  return static_cast<Number>(p) * static_cast<OtherNumber>(q);
-}
-
-
-
-/*---------------------- Inline functions: Tensor<1,dim> ---------------------*/
+  // TODO: Think about refactoring this into the TableIndices class as a
+  // general, polymorphic for extracting an item out of an object with
+  // nested identifiers.
+  template<int rank_> struct TensorIndicesHelper
+  {
+    // used for implementing Tensor<rank,dim>::operator[] with TableIndices
+    // tail recursive call to form up access to
+    //   tensor[indices[0]][indices[1]]...[indices[rank_]]
+    template<int rank, int dim, typename Number>
+    static inline
+    Number & extract(Tensor<rank_,dim,Number> &t, const TableIndices<rank> &indices)
+    {
+      Assert (indices[rank - rank_]<dim, ExcIndexRange (indices[rank - rank_], 0, dim));
+      return TensorIndicesHelper<rank_ - 1>::template extract<rank, dim, Number>(
+          t[indices[rank - rank_]], indices);
+    }
+  };
 
+  template<> struct TensorIndicesHelper<1>
+  {
+    template<int rank, int dim, typename Number>
+    static inline
+    Number & extract(Tensor<1,dim,Number> &t, const TableIndices<rank> &indices)
+    {
+      Assert (indices[rank - 1]<dim, ExcIndexRange (indices[rank - 1], 0, dim));
+      return t[indices[rank-1]];
+    }
+  };
+} /* internal */
 
 
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<1,dim,Number>::Tensor (const bool initialize)
+Tensor<rank_,dim,Number>::Tensor (const bool initialize)
 {
   if (initialize)
     // need to create an object Number() to initialize to zero to avoid
     // confusion with Tensor::operator=(scalar) when using something like
     // Tensor<1,dim,Tensor<1,dim,Number> >.
     for (unsigned int i=0; i!=dim; ++i)
-      values[i] = Number();
+      values[i] = value_type();
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<1,dim,Number>::Tensor (const array_type &initializer)
+Tensor<rank_,dim,Number>::Tensor (const Tensor<rank_,dim,Number> &initializer)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
@@ -981,32 +899,27 @@ Tensor<1,dim,Number>::Tensor (const array_type &initializer)
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,Number> &p)
+Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
   for (unsigned int i=0; i<dim; ++i)
-    values[i] = p.values[i];
+    values[i] = initializer[i];
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,OtherNumber> &p)
+Tensor<rank_,dim,Number>::Tensor (const Tensor<rank_,dim,OtherNumber> &initializer)
 {
-  Assert (dim>0, ExcDimTooSmall(dim));
-
-  for (unsigned int i=0; i<dim; ++i)
-    values[i] = Number(p.values[i]);
+  for (unsigned int i=0; i!=dim; ++i)
+    values[i] = initializer[i];
 }
 
 
-
 // At some places in the library, we have Point<0> for formal reasons
 // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
 // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
@@ -1019,108 +932,123 @@ Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &)
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber>
 inline
-Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const
+Tensor<rank_,dim,Number>::Tensor
+(const Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > &initializer)
 {
-  Assert (index<dim, ExcIndexRange (index, 0, dim));
-  return values[index];
+  for (unsigned int i=0; i<dim; ++i)
+    values[i] = initializer[i];
 }
 
 
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber>
+inline
+Tensor<rank_,dim,Number>::operator
+Tensor<1,dim,Tensor<rank_-1,dim,OtherNumber> > () const
+{
+  return Tensor<1,dim,Tensor<rank_-1,dim,Number> > (values);
+}
+
 
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Number &Tensor<1,dim,Number>::operator [] (const unsigned int index)
+typename Tensor<rank_,dim,Number>::value_type &
+Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
 {
-  Assert (index<dim, ExcIndexRange (index, 0, dim));
-  return values[index];
+  Assert (i<dim, ExcIndexRange(i, 0, dim));
+  return values[i];
 }
 
 
+template <int rank_, int dim, typename Number>
+inline
+const typename Tensor<rank_,dim,Number>::value_type &
+Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
+{
+  Assert (i<dim, ExcIndexRange(i, 0, dim));
+  return values[i];
+}
 
-template <int dim, typename Number>
+
+template <int rank_, int dim, typename Number>
 inline
-Number Tensor<1,dim,Number>::operator [] (const TableIndices<1> &indices) const
+Number
+Tensor<rank_,dim,Number>::operator[] (const TableIndices<rank_> &indices) const
 {
   Assert (indices[0]<dim, ExcIndexRange (indices[0], 0, dim));
-  return values[indices[0]];
+  return internal::TensorIndicesHelper<rank_>::extract(*this, indices);
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Number &Tensor<1,dim,Number>::operator [] (const TableIndices<1> &indices)
+Number &
+Tensor<rank_,dim,Number>::operator[] (const TableIndices<rank_> &indices)
 {
   Assert (indices[0]<dim, ExcIndexRange (indices[0], 0, dim));
-  return values[indices[0]];
+  return internal::TensorIndicesHelper<rank_>::extract(*this, indices);
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<1,dim,Number> &
-Tensor<1,dim,Number>::operator = (const Tensor<1,dim,Number> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,Number> &t)
 {
   for (unsigned int i=0; i<dim; ++i)
-    values[i] = p.values[i];
-
+    values[i] = t.values[i];
   return *this;
 }
 
 
-
+// At some places in the library, we have Point<0> for formal reasons
+// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
+// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
+// in the above function that the loop end check always fails, we
+// implement this function here
 template <>
 inline
 Tensor<1,0,double> &Tensor<1,0,double>::operator = (const Tensor<1,0,double> &)
 {
-  // at some places in the library, we have Point<0> for formal reasons
-  // (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
-  // Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
-  // in the above function that the loop end check always fails, we
-  // implement this function here
   return *this;
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number> &
-Tensor<1,dim,Number>::operator = (const Tensor<1,dim,OtherNumber> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,OtherNumber> &t)
 {
   for (unsigned int i=0; i<dim; ++i)
-    values[i] = Number(p.values[i]);
-
+    values[i] = t.values[i];
   return *this;
 }
 
 
-
-template <int dim, typename Number>
-template <typename OtherNumber>
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber, typename>
 inline
-Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator = (const OtherNumber d)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator = (const OtherNumber d)
 {
   Assert (d == OtherNumber(), ExcMessage ("Only assignment with zero is allowed"));
   (void) d;
 
   for (unsigned int i=0; i<dim; ++i)
     values[i] = Number();
-
   return *this;
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-bool Tensor<1,dim,Number>::operator == (const Tensor<1,dim,OtherNumber> &p) const
+bool
+Tensor<rank_,dim,Number>::operator == (const Tensor<rank_,dim,OtherNumber> &p) const
 {
   for (unsigned int i=0; i<dim; ++i)
     if (values[i] != p.values[i])
@@ -1129,7 +1057,11 @@ bool Tensor<1,dim,Number>::operator == (const Tensor<1,dim,OtherNumber> &p) cons
 }
 
 
-
+// At some places in the library, we have Point<0> for formal reasons
+// (e.g., we sometimes have Quadrature<dim-1> for faces, so we have
+// Quadrature<0> for dim=1, and then we have Point<0>). To avoid warnings
+// in the above function that the loop end check always fails, we
+// implement this function here
 template <>
 template <>
 inline
@@ -1139,21 +1071,21 @@ bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-bool Tensor<1,dim,Number>::operator != (const Tensor<1,dim,OtherNumber> &p) const
+bool
+Tensor<rank_,dim,Number>::operator != (const Tensor<rank_,dim,OtherNumber> &p) const
 {
   return !((*this) == p);
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator += (const Tensor<1,dim,OtherNumber> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator += (const Tensor<rank_,dim,OtherNumber> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] += p.values[i];
@@ -1161,11 +1093,11 @@ Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator += (const Tensor<1,dim,Othe
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,OtherNumber> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator -= (const Tensor<rank_,dim,OtherNumber> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] -= p.values[i];
@@ -1173,11 +1105,11 @@ Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,Othe
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator *= (const OtherNumber s)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator *= (const OtherNumber s)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] *= s;
@@ -1185,11 +1117,11 @@ Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator *= (const OtherNumber s)
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
-Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator /= (const OtherNumber s)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator /= (const OtherNumber s)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] /= s;
@@ -1197,118 +1129,259 @@ Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator /= (const OtherNumber s)
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - () const
+Tensor<rank_,dim,Number>
+Tensor<rank_,dim,Number>::operator - () const
 {
-  Tensor<1,dim,Number> result (false);
+  Tensor<rank_,dim,Number> tmp;
+
   for (unsigned int i=0; i<dim; ++i)
-    result.values[i] = -values[i];
-  return result;
-}
+    tmp.values[i] = -values[i];
 
+  return tmp;
+}
 
 
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-typename Tensor<1,dim,Number>::real_type
-Tensor<1,dim,Number>::norm () const
+typename Tensor<rank_,dim,Number>::real_type
+Tensor<rank_,dim,Number>::norm () const
 {
   return std::sqrt (norm_square());
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
-typename Tensor<1,dim,Number>::real_type
-Tensor<1,dim,Number>::norm_square () const
+typename Tensor<rank_,dim,Number>::real_type
+Tensor<rank_,dim,Number>::norm_square () const
 {
-  real_type s = numbers::NumberTraits<Number>::abs_square(values[0]);
-  for (unsigned int i=1; i<dim; ++i)
-    s += numbers::NumberTraits<Number>::abs_square(values[i]);
+  real_type s = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    s += static_cast<Tensor<rank_-1,dim,Number> >(values[i]).norm_square();
 
   return s;
 }
 
 
-
-template <int dim, typename Number>
-template <typename Number2>
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber>
 inline
 void
-Tensor<1,dim,Number>::unroll (Vector<Number2> &result) const
+Tensor<rank_, dim, Number>::unroll (Vector<OtherNumber> &result) const
 {
-  Assert (result.size()==dim,
-          ExcDimensionMismatch(dim, result.size()));
+  AssertDimension (result.size(),(Utilities::fixed_power<rank_, unsigned int>(dim)));
 
   unsigned int index = 0;
-  unroll_recursion (result,index);
+  unroll_recursion (result, index);
 }
 
 
-
-template<int dim, typename Number>
-template <typename Number2>
+template <int rank_, int dim, typename Number>
+template <typename OtherNumber>
 inline
 void
-Tensor<1,dim,Number>::unroll_recursion (Vector<Number2> &result,
-                                        unsigned int    &index) const
+Tensor<rank_, dim, Number>::unroll_recursion (Vector<OtherNumber> &result,
+                                              unsigned int        &index) const
 {
   for (unsigned int i=0; i<dim; ++i)
-    result(index++) = operator[](i);
+    static_cast<Tensor<rank_ - 1, dim, Number> >(values[i]).
+      unroll_recursion(result, index);
 }
 
 
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
 unsigned int
-Tensor<1, dim, Number>::component_to_unrolled_index (const TableIndices<1> &indices)
+Tensor<rank_, dim, Number>::component_to_unrolled_index(const TableIndices<rank_> &indices)
 {
-  return indices[0];
+  unsigned int index = 0;
+  for (int r = 0; r < rank_; ++r)
+    index = index * dim + indices[r];
+
+  return index;
 }
 
-template <int dim, typename Number>
+
+template <int rank_, int dim, typename Number>
 inline
-TableIndices<1>
-Tensor<1, dim, Number>::unrolled_to_component_indices (const unsigned int i)
+TableIndices<rank_>
+Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
 {
-  return TableIndices<1>(i);
-}
+  Assert (i < n_independent_components,
+          ExcIndexRange (i, 0, n_independent_components));
 
+  TableIndices<rank_>   indices;
 
+  unsigned int remainder = i;
+  for (int r=rank_-1; r>=0; --r)
+    {
+      indices[r] = (remainder % dim);
+      remainder /= dim;
+    }
+  Assert (remainder == 0, ExcInternalError());
 
-template <int dim, typename Number>
+  return indices;
+}
+
+
+template <int rank_, int dim, typename Number>
 inline
-void Tensor<1,dim,Number>::clear ()
+void Tensor<rank_,dim,Number>::clear ()
 {
   for (unsigned int i=0; i<dim; ++i)
-    values[i] = 0;
+    values[i] = value_type();
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 inline
 std::size_t
-Tensor<1,dim,Number>::memory_consumption ()
+Tensor<rank_,dim,Number>::memory_consumption ()
 {
-  return sizeof(Tensor<1,dim,Number>);
+  return sizeof(Tensor<rank_,dim,Number>);
 }
 
 
-
-template <int dim, typename Number>
+template <int rank_, int dim, typename Number>
 template <class Archive>
 inline
-void Tensor<1,dim,Number>::serialize(Archive &ar, const unsigned int)
+void
+Tensor<rank_,dim,Number>::serialize(Archive &ar, const unsigned int)
 {
   ar &values;
 }
-#endif // DOXYGEN
 
 
+#endif /* DOXYGEN */
+/* ----------------- Non-member functions operating on tensors. ------------- */
+
+
+
+#ifndef DEAL_II_WITH_CXX11
+
+template <typename T, typename U, int rank, int dim>
+struct ProductType<T,Tensor<rank,dim,U> >
+{
+  typedef Tensor<rank,dim,typename ProductType<T,U>::type> type;
+};
+
+template <typename T, typename U, int rank, int dim>
+struct ProductType<Tensor<rank,dim,T>,U>
+{
+  typedef Tensor<rank,dim,typename ProductType<T,U>::type> type;
+};
+
+#endif
+
+
+/**
+ * TODO
+ *
+ * @relates Tensor
+ * @relates EnableIfScalar
+ */
+template <int dim,
+         typename Number,
+         typename OtherNumber,
+         typename = typename EnableIfScalar<OtherNumber>::type>
+inline
+Tensor<0,dim,typename ProductType<OtherNumber, Number>::type>
+operator * (const OtherNumber           factor,
+            const Tensor<0,dim,Number> &t)
+{
+  return factor * static_cast<Number>(t);
+}
+
+
+
+/**
+ * TODO
+ *
+ * @relates Tensor
+ * @relates EnableIfScalar
+ */
+template <int dim,
+         typename Number,
+         typename OtherNumber,
+         typename = typename EnableIfScalar<OtherNumber>::type>
+inline
+Tensor<0,dim,typename ProductType<Number, OtherNumber>::type>
+operator * (const Tensor<0,dim,Number> &t,
+            const OtherNumber           factor)
+{
+  return static_cast<Number>(t) * factor;
+}
+
+
+
+/**
+ * TODO
+ *
+ * @relates Tensor
+ * @relates EnableIfScalar
+ */
+template <int dim,
+         typename Number,
+         typename OtherNumber,
+         typename = typename EnableIfScalar<OtherNumber>::type>
+inline
+Tensor<0,dim,typename ProductType<Number, OtherNumber>::type>
+operator / (const Tensor<0,dim,Number> &t,
+            const OtherNumber           factor)
+{
+  return static_cast<Number>(t) / factor;
+}
+
+
+
+/**
+ * Add two tensors of rank 0.
+ *
+ * @relates Tensor
+ */
+template <int dim, typename Number, typename OtherNumber>
+inline
+Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
+operator+ (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
+{
+  return static_cast<Number>(p) + static_cast<OtherNumber>(q);
+}
+
+
+
+/**
+ * Subtract two tensors of rank 0.
+ *
+ * @relates Tensor
+ */
+template <int dim, typename Number, typename OtherNumber>
+inline
+Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
+operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
+{
+  return static_cast<Number>(p) - static_cast<OtherNumber>(q);
+}
+
+
+
+/**
+ * Returns the contraction of two Tensors of rank 0.
+ *
+ * @relates Tensor
+ */
+template <int dim, typename Number, typename OtherNumber>
+inline
+typename ProductType<Number, OtherNumber>::type
+operator* (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q)
+{
+  return static_cast<Number>(p) * static_cast<OtherNumber>(q);
+}
+
+
+// TODO:
+
 
 /**
  * Output operator for tensors of rank 0. Since such tensors are scalars, we

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.