From 123952b525a6f25c81f03519a4276465713b7bde Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Tue, 1 Sep 2015 22:51:05 -0500 Subject: [PATCH] Eliminate specialization for Tensor<1,dim> This is more or less an exact code copy of Tensor. tensor_base.h now contains the full definition of Tensor. All external functions operating on tensors that are not strictly needed reside in tensor.h. --- include/deal.II/base/tensor.h | 692 --------------------- include/deal.II/base/tensor_base.h | 933 ++++++++++++++++------------- 2 files changed, 503 insertions(+), 1122 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 32eaae6b8b..224e3fc5fd 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -19,703 +19,11 @@ #include #include -#include DEAL_II_NAMESPACE_OPEN -template class Tensor; -template 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 - * contract 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 -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 dimension() 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::n_independent_components *dim; - - /** - * Type of stored objects. This is a tensor of lower rank. - */ - typedef Tensor 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, this - * corresponds to type number, and it is equal to Number for all other - * cases. See also the respective field in Vector. - * - * This typedef is used to represent the return type of norms. - */ - typedef typename numbers::NumberTraits::real_type real_type; - - /** - * Declare an array type which can be used to initialize an object of this - * type statically. - */ - typedef typename Tensor::array_type array_type[dim]; - - /** - * Constructor. Initialize all entries to zero if - * initialize==true; this is the default behaviour. - */ - explicit - Tensor (const bool initialize = true); - - /** - * Copy constructor. - */ - Tensor (const Tensor &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 - Tensor (const Tensor &initializer); - - /** - * Constructor that converts from a "tensor of tensors". - */ - template - Tensor (const Tensor<1,dim,Tensor > &initializer); - - /** - * Conversion operator to tensor of tensors. - */ - template - operator Tensor<1,dim,Tensor > () const; - - /** - * Read-Write access operator. - */ - Tensor &operator [] (const unsigned int i); - - /** - * Read-only access operator. - */ - const Tensor &operator [] (const unsigned int i) const; - - /** - * Read access using TableIndices indices - */ - Number operator [] (const TableIndices &indices) const; - - /** - * Read and write access using TableIndices indices - */ - Number &operator [] (const TableIndices &indices); - - /** - * Copy assignment operator. - */ - Tensor &operator = (const Tensor &rhs); - - /** - * Assignment operator from tensors with different underlying scalar type. - * This obviously requires that the @p OtherNumber type is convertible to @p - * Number. - */ - template - Tensor &operator = (const Tensor &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 d, allowing the intuitive notation - * t=0 to reset all elements of the tensor to zero. - * - * @relates EnableIfScalar - */ - template ::type> - Tensor &operator = (const OtherNumber d); - - /** - * Test for equality of two tensors. - */ - template - bool operator == (const Tensor &) const; - - /** - * Test for inequality of two tensors. - */ - template - bool operator != (const Tensor &) const; - - /** - * Add another tensor. - */ - template - Tensor &operator += (const Tensor &); - - /** - * Subtract another tensor. - */ - template - Tensor &operator -= (const Tensor &); - - /** - * Scale the tensor by factor, i.e. multiply all components by - * factor. - */ - template - Tensor &operator *= (const OtherNumber factor); - - /** - * Scale the vector by 1/factor. - */ - template - Tensor &operator /= (const OtherNumber factor); - - /** - * Unary minus operator. Negate all entries of a tensor. - */ - Tensor 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 - void unroll (Vector &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 &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 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 - void serialize(Archive &ar, const unsigned int version); - -private: - /** - * Array of tensors holding the subelements. - */ - Tensor subtensor[dim]; - - /** - * Help function for unroll. - */ - template - void unroll_recursion(Vector &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 friend class Tensor; -}; - - -/*--------------------------- Inline functions -----------------------------*/ - -#ifndef DOXYGEN - -template -inline -Tensor::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(); -} - - - -template -inline -Tensor::Tensor (const Tensor &initializer) -{ - for (unsigned int i=0; i!=dim; ++i) - subtensor[i] = initializer[i]; -} - - - -template -inline -Tensor::Tensor (const array_type &initializer) -{ - for (unsigned int i=0; i -template -inline -Tensor::Tensor (const Tensor &initializer) -{ - for (unsigned int i=0; i!=dim; ++i) - subtensor[i] = initializer[i]; -} - - - -template -template -inline -Tensor::Tensor -(const Tensor<1,dim,Tensor > &initializer) -{ - for (unsigned int i=0; i -template -inline -Tensor::operator -Tensor<1,dim,Tensor > () const -{ - return Tensor<1,dim,Tensor > (subtensor); -} - - - -template -inline -typename Tensor::value_type & -Tensor::operator[] (const unsigned int i) -{ - Assert (i -inline -const typename Tensor::value_type & -Tensor::operator[] (const unsigned int i) const -{ - Assert (i -inline -Number -Tensor::operator[] (const TableIndices &indices) const -{ - const unsigned int inner_ind = indices[0]; - Assert (inner_ind indices1; - for (unsigned int i = 0; i < rank_-1; i++) - indices1[i] = indices[i+1]; - return (subtensor[inner_ind])[indices1]; -} - - - -template -inline -Number & -Tensor::operator[] (const TableIndices &indices) -{ - const unsigned int inner_ind = indices[0]; - Assert (inner_ind indices1; - for (unsigned int i = 0; i < rank_-1; i++) - indices1[i] = indices[i+1]; - return (subtensor[inner_ind])[indices1]; -} - - - -template -inline -Tensor & -Tensor::operator = (const Tensor &t) -{ - for (unsigned int i=0; i -template -inline -Tensor & -Tensor::operator = (const Tensor &t) -{ - for (unsigned int i=0; i -template -inline -Tensor & -Tensor::operator = (const OtherNumber d) -{ - Assert (d == OtherNumber(), ExcMessage ("Only assignment with zero is allowed")); - (void) d; - - for (unsigned int i=0; i -template -inline -bool -Tensor::operator == (const Tensor &p) const -{ - for (unsigned int i=0; i -template -inline -bool -Tensor::operator != (const Tensor &p) const -{ - return !((*this) == p); -} - - - -template -template -inline -Tensor & -Tensor::operator += (const Tensor &p) -{ - for (unsigned int i=0; i -template -inline -Tensor & -Tensor::operator -= (const Tensor &p) -{ - for (unsigned int i=0; i -template -inline -Tensor & -Tensor::operator *= (const OtherNumber s) -{ - for (unsigned int i=0; i -template -inline -Tensor & -Tensor::operator /= (const OtherNumber s) -{ - for (unsigned int i=0; i -inline -Tensor -Tensor::operator - () const -{ - Tensor tmp; - - for (unsigned int i=0; i -inline -typename Tensor::real_type -Tensor::norm () const -{ - return std::sqrt (norm_square()); -} - - - -template -inline -typename Tensor::real_type -Tensor::norm_square () const -{ - real_type s = 0; - for (unsigned int i=0; i -template -inline -void -Tensor::unroll (Vector &result) const -{ - AssertDimension (result.size(),(Utilities::fixed_power(dim))); - - unsigned int index = 0; - unroll_recursion (result, index); -} -template -template -inline -void -Tensor::unroll_recursion (Vector &result, - unsigned int &index) const -{ - for (unsigned int i=0; i -inline -unsigned int -Tensor::component_to_unrolled_index(const TableIndices &indices) -{ - TableIndices 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::component_to_unrolled_index(indices1) * dim + indices[rank_-1]); -} - -template -inline -TableIndices -Tensor::unrolled_to_component_indices(const unsigned int i) -{ - Assert (i < n_independent_components, - ExcIndexRange (i, 0, n_independent_components)); - - TableIndices 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 -inline -void Tensor::clear () -{ - for (unsigned int i=0; i -inline -std::size_t -Tensor::memory_consumption () -{ - return sizeof(Tensor); -} - - - -template -template -inline -void -Tensor::serialize(Archive &ar, const unsigned int) -{ - ar &subtensor; -} - -#endif // DOXYGEN /* ----------------- Non-member functions operating on tensors. ------------ */ diff --git a/include/deal.II/base/tensor_base.h b/include/deal.II/base/tensor_base.h index 75735b46f3..ada55ccbd8 100644 --- a/include/deal.II/base/tensor_base.h +++ b/include/deal.II/base/tensor_base.h @@ -27,6 +27,7 @@ #include #include #include +#include #include #include @@ -43,11 +44,9 @@ template class Vector; // specified) template class Point; -// general template; specialized for rank==1; the general template is in -// tensor.h -template class Tensor; +// general template; specialized for rank == 0 +template class Tensor; template class Tensor<0,dim,Number>; -template 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::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 + 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 - void serialize(Archive &ar, const unsigned int version); - private: /** * The value of this scalar object. */ Number value; + + /** + * Help function for unroll. + */ + template + void unroll_recursion(Vector &result, + unsigned int &start_index) const; + + /** + * Allow an arbitrary Tensor to access the underlying values. + */ + template friend class Tensor; }; /** - * This class is a specialized version of the Tensor - * 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 + * contract family. * - * Within deal.II, the distinction between this class and its derived class - * Point is that we use the Point 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 - * Tensor<1,dim,Number> 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 Point 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 Tensor<1,dim,Number> 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 -class Tensor<1,dim,Number> +template +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::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, this * corresponds to type number, and it is equal to Number for all other * cases. See also the respective field in Vector. @@ -356,17 +379,26 @@ public: typedef typename numbers::NumberTraits::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 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::tensor_type value_type; /** - * Constructor. Initialize all entries to zero if initialize==true; - * this is the default behaviour. + * Declare an array type which can be used to initialize an object of this + * type statically. + */ + typedef typename Tensor::array_type + array_type[(dim != 0) ? dim : 1]; + + /** + * Constructor. Initialize all entries to zero if + * initialize==true; 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 &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 - Tensor (const Tensor<1,dim,OtherNumber> &initializer); + Tensor (const Tensor &initializer); /** - * Read access to the indexth coordinate. - * - * Note that the derived Point class also provides access through - * the () operator for backcompatibility. + * Constructor that converts from a "tensor of tensors". + */ + template + Tensor (const Tensor<1,dim,Tensor > &initializer); + + /** + * Conversion operator to tensor of tensors. + */ + template + operator Tensor<1,dim,Tensor > () 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 indexth coordinate. - * - * Note that the derived Point class also provides access through - * the () 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 indices */ - Number operator [] (const TableIndices<1> &indices) const; + Number operator [] (const TableIndices &indices) const; /** * Read and write access using TableIndices indices */ - Number &operator [] (const TableIndices<1> &indices); + Number &operator [] (const TableIndices &indices); /** * Copy assignment operator. */ - Tensor<1,dim,Number> &operator = (const Tensor<1,dim,Number> &rhs); + Tensor &operator = (const Tensor &rhs); /** * Assignment operator from tensors with different underlying scalar type. @@ -426,74 +464,61 @@ public: * Number. */ template - Tensor<1,dim,Number> &operator = (const Tensor<1,dim,OtherNumber> &rhs); + Tensor &operator = (const Tensor &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 d, allowing the intuitive notation * t=0 to reset all elements of the tensor to zero. + * + * @relates EnableIfScalar */ - template - Tensor<1,dim,Number> &operator = (const OtherNumber d); + template ::type> + Tensor &operator = (const OtherNumber d); /** * Test for equality of two tensors. */ template - bool operator == (const Tensor<1,dim,OtherNumber> &rhs) const; + bool operator == (const Tensor &) const; /** * Test for inequality of two tensors. */ template - bool operator != (const Tensor<1,dim,OtherNumber> &rhs) const; + bool operator != (const Tensor &) const; /** - * Add another vector to this vector. + * Add another tensor. */ template - Tensor<1,dim,Number> &operator += (const Tensor<1,dim,OtherNumber> &rhs); + Tensor &operator += (const Tensor &); /** - * Subtract another vector. + * Subtract another tensor. */ template - Tensor<1,dim,Number> &operator -= (const Tensor<1,dim,OtherNumber> &rhs); + Tensor &operator -= (const Tensor &); /** - * Scale the vector by factor, i.e., multiply all coordinates by + * Scale the tensor by factor, i.e. multiply all components by * factor. */ template - Tensor<1,dim,Number> &operator *= (const OtherNumber factor); + Tensor &operator *= (const OtherNumber factor); /** * Scale the vector by 1/factor. */ template - 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 l2 norm of the vector. - */ - real_type norm () const; + Tensor &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 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 - void unroll (Vector &result) const; + template + void unroll (Vector &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 &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 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 + 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 dim 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 - void serialize(Archive &ar, const unsigned int version); - private: /** - * Store the values in a simple array. For dim==0 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 - void unroll_recursion (Vector &result, - unsigned int &start_index) const; + template + void unroll_recursion(Vector &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 friend class dealii::Tensor; + template 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 -std::ostream &operator << (std::ostream &out, const Tensor<0,dim,Number> &p); - -/** - * Prints the values of this tensor in the form x1 x2 x3 etc. - */ -template -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 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 template inline @@ -654,7 +670,6 @@ Tensor<0,dim,Number>::Tensor (const OtherNumber initializer) } - template template inline @@ -666,7 +681,6 @@ Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,OtherNumber> &p) } - template inline Tensor<0,dim,Number>::operator Number () const @@ -675,7 +689,6 @@ Tensor<0,dim,Number>::operator Number () const } - template inline Tensor<0,dim,Number>::operator Number &() @@ -684,7 +697,6 @@ Tensor<0,dim,Number>::operator Number &() } - template 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 template inline @@ -705,7 +716,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Other } - template template inline @@ -716,7 +726,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const OtherNumber d) } - template template inline @@ -726,7 +735,6 @@ bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,OtherNumber> &p) cons } - template template inline @@ -736,7 +744,6 @@ bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,OtherNumber> &p) cons } - template template inline @@ -747,7 +754,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator += (const Tensor<0,dim,Othe } - template template inline @@ -758,7 +764,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Othe } - template template inline @@ -769,7 +774,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s) } - template template inline @@ -780,7 +784,6 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator /= (const OtherNumber s) } - template 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 inline typename Tensor<0,dim,Number>::real_type @@ -799,7 +801,6 @@ Tensor<0,dim,Number>::norm () const } - template inline typename Tensor<0,dim,Number>::real_type @@ -809,6 +810,17 @@ Tensor<0,dim,Number>::norm_square () const } +template +template +inline +void +Tensor<0, dim, Number>::unroll_recursion (Vector &result, + unsigned int &index) const +{ + result[index] = value; + ++index; +} + template inline @@ -818,7 +830,6 @@ void Tensor<0,dim,Number>::clear () } - template template inline @@ -828,151 +839,58 @@ void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int) } - -#ifndef DEAL_II_WITH_CXX11 - -template -struct ProductType > -{ - typedef Tensor::type> type; -}; - -template -struct ProductType,U> -{ - typedef Tensor::type> type; -}; - -#endif - - - -/** - * TODO - * - * @relates Tensor - * @relates EnableIfScalar - */ -template ::type> -inline -Tensor<0,dim,typename ProductType::type> -operator * (const OtherNumber factor, - const Tensor<0,dim,Number> &t) -{ - return factor * static_cast(t); -} - - - -/** - * TODO - * - * @relates Tensor - * @relates EnableIfScalar - */ -template ::type> -inline -Tensor<0,dim,typename ProductType::type> -operator * (const Tensor<0,dim,Number> &t, - const OtherNumber factor) -{ - return static_cast(t) * factor; -} - - - -/** - * TODO - * - * @relates Tensor - * @relates EnableIfScalar - */ -template ::type> -inline -Tensor<0,dim,typename ProductType::type> -operator / (const Tensor<0,dim,Number> &t, - const OtherNumber factor) -{ - return static_cast(t) / factor; -} - - - -/** - * Add two tensors of rank 0. - * - * @relates Tensor - */ -template -inline -Tensor<0, dim, typename ProductType::type> -operator+ (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) -{ - return static_cast(p) + static_cast(q); -} - - - -/** - * Subtract two tensors of rank 0. - * - * @relates Tensor - */ -template -inline -Tensor<0, dim, typename ProductType::type> -operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) -{ - return static_cast(p) - static_cast(q); -} - +/*-------------------- Inline functions: Tensor --------------------*/ -/** - * Returns the contraction of two Tensors of rank 0. - * - * @relates Tensor - */ -template -inline -typename ProductType::type -operator* (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) +namespace internal { - return static_cast(p) * static_cast(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 struct TensorIndicesHelper + { + // used for implementing Tensor::operator[] with TableIndices + // tail recursive call to form up access to + // tensor[indices[0]][indices[1]]...[indices[rank_]] + template + static inline + Number & extract(Tensor &t, const TableIndices &indices) + { + Assert (indices[rank - rank_]::template extract( + t[indices[rank - rank_]], indices); + } + }; + template<> struct TensorIndicesHelper<1> + { + template + static inline + Number & extract(Tensor<1,dim,Number> &t, const TableIndices &indices) + { + Assert (indices[rank - 1] +template inline -Tensor<1,dim,Number>::Tensor (const bool initialize) +Tensor::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 +template inline -Tensor<1,dim,Number>::Tensor (const array_type &initializer) +Tensor::Tensor (const Tensor &initializer) { Assert (dim>0, ExcDimTooSmall(dim)); @@ -981,32 +899,27 @@ Tensor<1,dim,Number>::Tensor (const array_type &initializer) } - -template +template inline -Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,Number> &p) +Tensor::Tensor (const array_type &initializer) { Assert (dim>0, ExcDimTooSmall(dim)); for (unsigned int i=0; i +template template inline -Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,OtherNumber> &p) +Tensor::Tensor (const Tensor &initializer) { - Assert (dim>0, ExcDimTooSmall(dim)); - - for (unsigned int i=0; i for formal reasons // (e.g., we sometimes have Quadrature 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 +template +template inline -Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const +Tensor::Tensor +(const Tensor<1,dim,Tensor > &initializer) { - Assert (index +template +inline +Tensor::operator +Tensor<1,dim,Tensor > () const +{ + return Tensor<1,dim,Tensor > (values); +} + -template +template inline -Number &Tensor<1,dim,Number>::operator [] (const unsigned int index) +typename Tensor::value_type & +Tensor::operator[] (const unsigned int i) { - Assert (index +inline +const typename Tensor::value_type & +Tensor::operator[] (const unsigned int i) const +{ + Assert (i + +template inline -Number Tensor<1,dim,Number>::operator [] (const TableIndices<1> &indices) const +Number +Tensor::operator[] (const TableIndices &indices) const { Assert (indices[0]::extract(*this, indices); } - -template +template inline -Number &Tensor<1,dim,Number>::operator [] (const TableIndices<1> &indices) +Number & +Tensor::operator[] (const TableIndices &indices) { Assert (indices[0]::extract(*this, indices); } - -template +template inline -Tensor<1,dim,Number> & -Tensor<1,dim,Number>::operator = (const Tensor<1,dim,Number> &p) +Tensor & +Tensor::operator = (const Tensor &t) { for (unsigned int i=0; i for formal reasons +// (e.g., we sometimes have Quadrature 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 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 +template template inline -Tensor<1,dim,Number> & -Tensor<1,dim,Number>::operator = (const Tensor<1,dim,OtherNumber> &p) +Tensor & +Tensor::operator = (const Tensor &t) { for (unsigned int i=0; i -template +template +template inline -Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator = (const OtherNumber d) +Tensor & +Tensor::operator = (const OtherNumber d) { Assert (d == OtherNumber(), ExcMessage ("Only assignment with zero is allowed")); (void) d; for (unsigned int i=0; i +template template inline -bool Tensor<1,dim,Number>::operator == (const Tensor<1,dim,OtherNumber> &p) const +bool +Tensor::operator == (const Tensor &p) const { for (unsigned int i=0; i::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 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 +template template inline -bool Tensor<1,dim,Number>::operator != (const Tensor<1,dim,OtherNumber> &p) const +bool +Tensor::operator != (const Tensor &p) const { return !((*this) == p); } - -template +template template inline -Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator += (const Tensor<1,dim,OtherNumber> &p) +Tensor & +Tensor::operator += (const Tensor &p) { for (unsigned int i=0; i &Tensor<1,dim,Number>::operator += (const Tensor<1,dim,Othe } - -template +template template inline -Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,OtherNumber> &p) +Tensor & +Tensor::operator -= (const Tensor &p) { for (unsigned int i=0; i &Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,Othe } - -template +template template inline -Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator *= (const OtherNumber s) +Tensor & +Tensor::operator *= (const OtherNumber s) { for (unsigned int i=0; i &Tensor<1,dim,Number>::operator *= (const OtherNumber s) } - -template +template template inline -Tensor<1,dim,Number> &Tensor<1,dim,Number>::operator /= (const OtherNumber s) +Tensor & +Tensor::operator /= (const OtherNumber s) { for (unsigned int i=0; i &Tensor<1,dim,Number>::operator /= (const OtherNumber s) } - -template +template inline -Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - () const +Tensor +Tensor::operator - () const { - Tensor<1,dim,Number> result (false); + Tensor tmp; + for (unsigned int i=0; i +template inline -typename Tensor<1,dim,Number>::real_type -Tensor<1,dim,Number>::norm () const +typename Tensor::real_type +Tensor::norm () const { return std::sqrt (norm_square()); } - -template +template inline -typename Tensor<1,dim,Number>::real_type -Tensor<1,dim,Number>::norm_square () const +typename Tensor::real_type +Tensor::norm_square () const { - real_type s = numbers::NumberTraits::abs_square(values[0]); - for (unsigned int i=1; i::abs_square(values[i]); + real_type s = 0; + for (unsigned int i=0; i >(values[i]).norm_square(); return s; } - -template -template +template +template inline void -Tensor<1,dim,Number>::unroll (Vector &result) const +Tensor::unroll (Vector &result) const { - Assert (result.size()==dim, - ExcDimensionMismatch(dim, result.size())); + AssertDimension (result.size(),(Utilities::fixed_power(dim))); unsigned int index = 0; - unroll_recursion (result,index); + unroll_recursion (result, index); } - -template -template +template +template inline void -Tensor<1,dim,Number>::unroll_recursion (Vector &result, - unsigned int &index) const +Tensor::unroll_recursion (Vector &result, + unsigned int &index) const { for (unsigned int i=0; i >(values[i]). + unroll_recursion(result, index); } -template +template inline unsigned int -Tensor<1, dim, Number>::component_to_unrolled_index (const TableIndices<1> &indices) +Tensor::component_to_unrolled_index(const TableIndices &indices) { - return indices[0]; + unsigned int index = 0; + for (int r = 0; r < rank_; ++r) + index = index * dim + indices[r]; + + return index; } -template + +template inline -TableIndices<1> -Tensor<1, dim, Number>::unrolled_to_component_indices (const unsigned int i) +TableIndices +Tensor::unrolled_to_component_indices(const unsigned int i) { - return TableIndices<1>(i); -} + Assert (i < n_independent_components, + ExcIndexRange (i, 0, n_independent_components)); + TableIndices indices; + unsigned int remainder = i; + for (int r=rank_-1; r>=0; --r) + { + indices[r] = (remainder % dim); + remainder /= dim; + } + Assert (remainder == 0, ExcInternalError()); -template + return indices; +} + + +template inline -void Tensor<1,dim,Number>::clear () +void Tensor::clear () { for (unsigned int i=0; i +template inline std::size_t -Tensor<1,dim,Number>::memory_consumption () +Tensor::memory_consumption () { - return sizeof(Tensor<1,dim,Number>); + return sizeof(Tensor); } - -template +template template inline -void Tensor<1,dim,Number>::serialize(Archive &ar, const unsigned int) +void +Tensor::serialize(Archive &ar, const unsigned int) { ar &values; } -#endif // DOXYGEN +#endif /* DOXYGEN */ +/* ----------------- Non-member functions operating on tensors. ------------- */ + + + +#ifndef DEAL_II_WITH_CXX11 + +template +struct ProductType > +{ + typedef Tensor::type> type; +}; + +template +struct ProductType,U> +{ + typedef Tensor::type> type; +}; + +#endif + + +/** + * TODO + * + * @relates Tensor + * @relates EnableIfScalar + */ +template ::type> +inline +Tensor<0,dim,typename ProductType::type> +operator * (const OtherNumber factor, + const Tensor<0,dim,Number> &t) +{ + return factor * static_cast(t); +} + + + +/** + * TODO + * + * @relates Tensor + * @relates EnableIfScalar + */ +template ::type> +inline +Tensor<0,dim,typename ProductType::type> +operator * (const Tensor<0,dim,Number> &t, + const OtherNumber factor) +{ + return static_cast(t) * factor; +} + + + +/** + * TODO + * + * @relates Tensor + * @relates EnableIfScalar + */ +template ::type> +inline +Tensor<0,dim,typename ProductType::type> +operator / (const Tensor<0,dim,Number> &t, + const OtherNumber factor) +{ + return static_cast(t) / factor; +} + + + +/** + * Add two tensors of rank 0. + * + * @relates Tensor + */ +template +inline +Tensor<0, dim, typename ProductType::type> +operator+ (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) +{ + return static_cast(p) + static_cast(q); +} + + + +/** + * Subtract two tensors of rank 0. + * + * @relates Tensor + */ +template +inline +Tensor<0, dim, typename ProductType::type> +operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) +{ + return static_cast(p) - static_cast(q); +} + + + +/** + * Returns the contraction of two Tensors of rank 0. + * + * @relates Tensor + */ +template +inline +typename ProductType::type +operator* (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) +{ + return static_cast(p) * static_cast(q); +} + + +// TODO: + /** * Output operator for tensors of rank 0. Since such tensors are scalars, we -- 2.39.5