#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. ------------ */
#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>
// 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>;
/**
* 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
*/
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;
*/
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
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
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
* @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:
/**
/**
* 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>.
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);
/**
* 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.
* 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.
*/
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
*/
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
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
};
-/**
- * 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> ---------------------*/
}
-
template <int dim, typename Number>
inline
Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p)
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
inline
Tensor<0,dim,Number>::operator Number () const
}
-
template <int dim, typename Number>
inline
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)
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
template <typename OtherNumber>
inline
}
-
template <int dim, typename Number>
inline
Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const
}
-
template <int dim, typename Number>
inline
typename Tensor<0,dim,Number>::real_type
}
-
template <int dim, typename Number>
inline
typename Tensor<0,dim,Number>::real_type
}
+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
}
-
template <int dim, typename Number>
template <class Archive>
inline
}
-
-#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));
}
-
-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
}
-
-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])
}
-
+// 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
}
-
-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];
}
-
-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];
}
-
-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;
}
-
-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;
}
-
-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
projects / dealii.git / commitdiff