From 0accc9110366b159e64c42428ab5969e493ee3dd Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sun, 6 Sep 2015 21:15:41 -0500 Subject: [PATCH] Obsolete tensor_base.h, refactor everything to tensor.h --- include/deal.II/base/tensor.h | 1307 +++++++++++++++++++++++++- include/deal.II/base/tensor_base.h | 1361 ---------------------------- 2 files changed, 1305 insertions(+), 1363 deletions(-) delete mode 100644 include/deal.II/base/tensor_base.h diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index d386528603..4e0cf9fdf2 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -16,15 +16,1193 @@ #ifndef dealii__tensor_h #define dealii__tensor_h - #include -#include +#include +#include +#include +#include + +#include +#include +#include DEAL_II_NAMESPACE_OPEN +// Forward declarations: + +template class Point; +template class Tensor; + + + +/** + * This class is a specialized version of the + * Tensor class. It handles tensors of rank zero, + * i.e. scalars. The second template argument @param dim is ignored. + * + * This class exists because in some cases we want to construct objects of + * type Tensor@, which should expand to scalars, + * vectors, matrices, etc, depending on the values of the template arguments + * @p dim and @p spacedim. We therefore need a class that acts as a scalar + * (i.e. @p Number) for all purposes but is part of the Tensor template + * family. + * + * @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. Since the current object is a rank-0 + * tensor (a scalar), this template argument has no meaning for this class. + * + * @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, 2009, Matthias Maier, 2015 + */ +template +class Tensor<0,dim,Number> +{ +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 = 0; + + /** + * Number of independent components of a tensor of rank 0. + */ + static const unsigned int n_independent_components = 1; + + /** + * 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; + + /** + * Type of objects encapsulated by this container and returned by + * operator[](). This is a scalar number type 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 the scalar + * number type Number. + */ + typedef Number array_type; + + /** + * Constructor. Set to zero. + */ + Tensor (); + + /** + * Copy constructor. + */ + Tensor (const Tensor<0,dim,Number> &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<0,dim,OtherNumber> &initializer); + + /** + * Constructor, where the data is copied from a C-style array. + */ + template + Tensor (const OtherNumber initializer); + + /** + * Return a reference to the encapsulated Number object. Since rank-0 + * tensors are scalars, this is a natural operation. + * + * This is the non-const conversion operator that returns a writable + * reference. + */ + operator Number &(); + + /** + * Return a reference to the encapsulated Number object. Since rank-0 + * tensors are scalars, this is a natural operation. + * + * This is the const conversion operator that returns a read-only + * reference. + */ + operator const Number &() const; + + /** + * Copy assignment operator. + */ + Tensor<0,dim,Number> &operator = (const Tensor<0,dim,Number> &rhs); + + /** + * Assignment from tensors with different underlying scalar type. + * This obviously requires that the @p OtherNumber type is convertible to @p + * Number. + */ + template + Tensor<0,dim,Number> &operator = (const Tensor<0,dim,OtherNumber> &rhs); + + /** + * Test for equality of two tensors. + */ + template + bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const; + + /** + * Test for inequality of two tensors. + */ + template + bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const; + + /** + * Add another scalar + */ + template + Tensor<0,dim,Number> &operator += (const Tensor<0,dim,OtherNumber> &rhs); + + /** + * Subtract another scalar. + */ + template + Tensor<0,dim,Number> &operator -= (const Tensor<0,dim,OtherNumber> &rhs); + + /** + * Multiply the scalar with a factor. + */ + template + Tensor<0,dim,Number> &operator *= (const OtherNumber factor); + + /** + * Divide the scalar by factor. + */ + template + Tensor<0,dim,Number> &operator /= (const OtherNumber factor); + + /** + * Tensor with inverted entries. + */ + 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 the absolute 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; + + /** + * Return the square of the Frobenius-norm of a tensor, i.e. the sum of + * the absolute squares of all entries. + */ + real_type norm_square () const; + + /** + * 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: + /** + * Internal type declaration that is used to specialize the return type + * of operator[]() for Tensor<1,dim,Number> + */ + typedef Number tensor_type; + + /** + * 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; +}; + + + +/** + * 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, Matthias Maier, 2015 + */ +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; + + /** + * 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. + * + * This typedef is used to represent the return type of norms. + */ + typedef typename numbers::NumberTraits::real_type real_type; + + /** + * Type of objects encapsulated by this container and returned by + * operator[](). This is a tensor of lower rank for a general tensor, and + * a scalar number type for Tensor<1,dim,Number>. + */ + typedef typename Tensor::tensor_type value_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 != 0) ? dim : 1]; + + /** + * 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. + */ + value_type &operator [] (const unsigned int i); + + /** + * Read-only access operator. + */ + const value_type &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. + */ + Tensor &operator = (const Number 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; + + /** + * 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 the absolute 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; + + /** + * Return the square of the Frobenius-norm of a tensor, i.e. the sum of + * the absolute squares of all entries. + */ + 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); + + /** + * Determine an estimate for the memory consumption (in bytes) of this + * object. + */ + 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."); + +private: + /** + * Internal type declaration that is used to specialize the return type + * of operator[]() for Tensor<1,dim,Number> + */ + typedef Tensor tensor_type; + + /** + * Array of tensors holding the subelements. + */ + Tensor values[(dim != 0) ? dim : 1]; + + /** + * 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; + + /** + * Point is allowed access to the coordinates. This is supposed to improve + * speed. + */ + friend class Point; +}; + + + +#ifndef DOXYGEN +/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/ + + +template +inline +Tensor<0,dim,Number>::Tensor () +{ + value = value_type(); +} + + +template +inline +Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p) +{ + value = p.value; +} + + +template +template +inline +Tensor<0,dim,Number>::Tensor (const OtherNumber initializer) +{ + value = initializer; +} + + +template +template +inline +Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,OtherNumber> &p) +{ + value = p.value; +} + + +template +inline +Tensor<0,dim,Number>::operator Number &() +{ + return value; +} + + +template +inline +Tensor<0,dim,Number>::operator const Number &() const +{ + return value; +} + + +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p) +{ + value = p.value; + return *this; +} + + +template +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,OtherNumber> &p) +{ + value = p.value; + return *this; +} + + +template +template +inline +bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,OtherNumber> &p) const +{ + return (value == p.value); +} + + +template +template +inline +bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,OtherNumber> &p) const +{ + return !((*this) == p); +} + + +template +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator += (const Tensor<0,dim,OtherNumber> &p) +{ + value += p.value; + return *this; +} + + +template +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,OtherNumber> &p) +{ + value -= p.value; + return *this; +} + + +template +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s) +{ + value *= s; + return *this; +} + + +template +template +inline +Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator /= (const OtherNumber s) +{ + value /= s; + return *this; +} + + +template +inline +Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const +{ + return -value; +} + + +template +inline +typename Tensor<0,dim,Number>::real_type +Tensor<0,dim,Number>::norm () const +{ + return numbers::NumberTraits::abs (value); +} + + +template +inline +typename Tensor<0,dim,Number>::real_type +Tensor<0,dim,Number>::norm_square () const +{ + return numbers::NumberTraits::abs_square (value); +} + + +template +template +inline +void +Tensor<0, dim, Number>::unroll_recursion (Vector &result, + unsigned int &index) const +{ + result[index] = value; + ++index; +} + + +template +inline +void Tensor<0,dim,Number>::clear () +{ + value = value_type(); +} + + +template +template +inline +void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int) +{ + ar &value; +} + + +/*-------------------- Inline functions: Tensor --------------------*/ + + +namespace internal +{ + // 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] +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) + values[i] = value_type(); +} + + +template +inline +Tensor::Tensor (const Tensor &initializer) +{ + for (unsigned int i=0; i +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) + values[i] = initializer[i]; +} + + +// 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 <> +inline +Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &) +{ +} + + +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 > (values); +} + + +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 +{ + Assert (indices[0]::extract(*this, indices); +} + + +template +inline +Number & +Tensor::operator[] (const TableIndices &indices) +{ + Assert (indices[0]::extract(*this, indices); +} + + +template +inline +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> &) +{ + return *this; +} + + +template +template +inline +Tensor & +Tensor::operator = (const Tensor &t) +{ + for (unsigned int i=0; i +inline +Tensor & +Tensor::operator = (const Number d) +{ + Assert (d == Number(), 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 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 +bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const +{ + return true; +} + + +template +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) +{ + unsigned int index = 0; + for (int r = 0; r < rank_; ++r) + index = index * dim + indices[r]; + + return index; +} + + +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 &values; +} + +#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 + /** * @name Output functions for Tensor objects */ @@ -158,6 +1336,131 @@ operator- (const Tensor<0,dim,Number> &p, const Tensor<0,dim,OtherNumber> &q) } +/** + * Multiplication of a tensor of general rank with a scalar number from the + * right. + * + * The purpose of this operator is to enable only multiplication of a tensor + * by a scalar number (i.e., a floating point number, a complex floating point + * number, etc.). The function is written in a way that only allows the + * compiler to consider the function if the second argument is indeed a scalar + * number -- in other words, @p OtherNumber will not match, for example + * std::vector@ as the product of a tensor and a vector + * clearly would make no sense. The mechanism by which the compiler is + * prohibited of considering this operator for multiplication with non-scalar + * types are explained in the documentation of the EnableIfScalar class. + * + * The return type of the function is chosen so that it matches the types of + * both the tensor and the scalar argument. For example, if you multiply a + * Tensor@<1,dim,double@> by std::complex@, + * then the result will be a + * Tensor@<1,dim,std::complex@@>. In other words, the + * type with which the returned tensor stores its components equals the type + * you would get if you multiplied an individual component of the input tensor + * by the scalar factor. + * + * @relates Tensor + * @relates EnableIfScalar + */ +template +inline +Tensor::type>::type> +operator * (const Tensor &t, + const OtherNumber factor) +{ + // recurse over the base objects + Tensor::type> tt; + for (unsigned int d=0; d +inline +Tensor::type, OtherNumber>::type> +operator * (const Number factor, + const Tensor &t) +{ + // simply forward to the operator above + return t * factor; +} + + +/** + * Division of a tensor of general rank with a scalar number. See the + * discussion on operator*() above for more information about template + * arguments and the return type. + * + * @relates Tensor + * @relates EnableIfScalar + */ +template +inline +Tensor::type>::type> +operator / (const Tensor &t, + const OtherNumber factor) +{ + // recurse over the base objects + Tensor::type> tt; + for (unsigned int d=0; d +inline +Tensor::type> +operator+ (const Tensor &p, const Tensor &q) +{ + Tensor::type> tmp (p); + + for (unsigned int i=0; i +inline +Tensor::type> +operator- (const Tensor &p, const Tensor &q) +{ + Tensor::type> tmp (p); + + for (unsigned int i=0; i -#include -#include -#include -#include -#include - -#include -#include - -DEAL_II_NAMESPACE_OPEN - -template class Vector; - -// forward declare Point and Tensor. This is the first definition of these -// classes and here we set the default Number type to double (this means that -// this file must be included when using something like Tensor<1,dim>, and -// Point and Tensor must not be forward declared without the number type -// specified) -template class Point; - -// general template; specialized for rank == 0 -template class Tensor; -template class Tensor<0,dim,Number>; - - - -/** - * This class is a specialized version of the - * Tensor class. It handles tensors of rank zero, - * i.e. scalars. The second template argument @param dim is ignored. - * - * This class exists because in some cases we want to construct objects of - * type Tensor@, which should expand to scalars, - * vectors, matrices, etc, depending on the values of the template arguments - * @p dim and @p spacedim. We therefore need a class that acts as a scalar - * (i.e. @p Number) for all purposes but is part of the Tensor template - * family. - * - * @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. Since the current object is a rank-0 - * tensor (a scalar), this template argument has no meaning for this class. - * - * @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, 2009, Matthias Maier, 2015 - */ -template -class Tensor<0,dim,Number> -{ -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 = 0; - - /** - * Number of independent components of a tensor of rank 0. - */ - static const unsigned int n_independent_components = 1; - - /** - * 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; - - /** - * Type of objects encapsulated by this container and returned by - * operator[](). This is a scalar number type 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 the scalar - * number type Number. - */ - typedef Number array_type; - - /** - * Constructor. Set to zero. - */ - Tensor (); - - /** - * Copy constructor. - */ - Tensor (const Tensor<0,dim,Number> &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<0,dim,OtherNumber> &initializer); - - /** - * Constructor, where the data is copied from a C-style array. - */ - template - Tensor (const OtherNumber initializer); - - /** - * Return a reference to the encapsulated Number object. Since rank-0 - * tensors are scalars, this is a natural operation. - * - * This is the non-const conversion operator that returns a writable - * reference. - */ - operator Number &(); - - /** - * Return a reference to the encapsulated Number object. Since rank-0 - * tensors are scalars, this is a natural operation. - * - * This is the const conversion operator that returns a read-only - * reference. - */ - operator const Number &() const; - - /** - * Copy assignment operator. - */ - Tensor<0,dim,Number> &operator = (const Tensor<0,dim,Number> &rhs); - - /** - * Assignment from tensors with different underlying scalar type. - * This obviously requires that the @p OtherNumber type is convertible to @p - * Number. - */ - template - Tensor<0,dim,Number> &operator = (const Tensor<0,dim,OtherNumber> &rhs); - - /** - * Test for equality of two tensors. - */ - template - bool operator == (const Tensor<0,dim,OtherNumber> &rhs) const; - - /** - * Test for inequality of two tensors. - */ - template - bool operator != (const Tensor<0,dim,OtherNumber> &rhs) const; - - /** - * Add another scalar - */ - template - Tensor<0,dim,Number> &operator += (const Tensor<0,dim,OtherNumber> &rhs); - - /** - * Subtract another scalar. - */ - template - Tensor<0,dim,Number> &operator -= (const Tensor<0,dim,OtherNumber> &rhs); - - /** - * Multiply the scalar with a factor. - */ - template - Tensor<0,dim,Number> &operator *= (const OtherNumber factor); - - /** - * Divide the scalar by factor. - */ - template - Tensor<0,dim,Number> &operator /= (const OtherNumber factor); - - /** - * Tensor with inverted entries. - */ - 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 the absolute 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; - - /** - * Return the square of the Frobenius-norm of a tensor, i.e. the sum of - * the absolute squares of all entries. - */ - real_type norm_square () const; - - /** - * 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: - /** - * Internal type declaration that is used to specialize the return type - * of operator[]() for Tensor<1,dim,Number> - */ - typedef Number tensor_type; - - /** - * 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; -}; - - - -/** - * 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, Matthias Maier, 2015 - */ -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; - - /** - * 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. - * - * This typedef is used to represent the return type of norms. - */ - typedef typename numbers::NumberTraits::real_type real_type; - - /** - * Type of objects encapsulated by this container and returned by - * operator[](). This is a tensor of lower rank for a general tensor, and - * a scalar number type for Tensor<1,dim,Number>. - */ - typedef typename Tensor::tensor_type value_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 != 0) ? dim : 1]; - - /** - * 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. - */ - value_type &operator [] (const unsigned int i); - - /** - * Read-only access operator. - */ - const value_type &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. - */ - Tensor &operator = (const Number 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; - - /** - * 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 the absolute 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; - - /** - * Return the square of the Frobenius-norm of a tensor, i.e. the sum of - * the absolute squares of all entries. - */ - 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); - - /** - * Determine an estimate for the memory consumption (in bytes) of this - * object. - */ - 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."); - -private: - /** - * Internal type declaration that is used to specialize the return type - * of operator[]() for Tensor<1,dim,Number> - */ - typedef Tensor tensor_type; - - /** - * Array of tensors holding the subelements. - */ - Tensor values[(dim != 0) ? dim : 1]; - - /** - * 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; - - /** - * Point is allowed access to the coordinates. This is supposed to improve - * speed. - */ - friend class Point; -}; - - -#ifndef DOXYGEN -/*---------------------- Inline functions: Tensor<0,dim> ---------------------*/ - - -template -inline -Tensor<0,dim,Number>::Tensor () -{ - value = value_type(); -} - - -template -inline -Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p) -{ - Assert(dim != 0 || p.value == Number(), - ExcMessage("Creation of a Tensor<0,0,Number> object with a non-zero scalar requested.")); - value = p.value; -} - - -template -template -inline -Tensor<0,dim,Number>::Tensor (const OtherNumber initializer) -{ - Assert(dim != 0 || initializer == OtherNumber(), - ExcMessage("Creation of a Tensor<0,0,Number> object with a non-zero scalar requested.")); - value = initializer; -} - - -template -template -inline -Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,OtherNumber> &p) -{ - Assert(dim != 0 || p.value == OtherNumber(), - ExcMessage("Cannot return a non-zero scalar from a Tensor<0,0,Number> object.")); - value = p.value; -} - - -template -inline -Tensor<0,dim,Number>::operator Number &() -{ - Assert(dim != 0 || value == Number(), - ExcMessage("Cannot return a non-zero scalar from a Tensor<0,0,Number> object.")); - return value; -} - - -template -inline -Tensor<0,dim,Number>::operator const Number &() const -{ - Assert(dim != 0 || value == Number(), - ExcMessage("Cannot assign a non-zero scalar to a Tensor<0,0,Number> object.")); - return value; -} - - -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p) -{ - Assert(dim != 0 || p.value == Number(), - ExcMessage("Cannot assign a non-zero scalar to a Tensor<0,0,Number> object.")); - value = p.value; - return *this; -} - - -template -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator = (const Tensor<0,dim,OtherNumber> &p) -{ - Assert(dim != 0 || p.value == OtherNumber(), - ExcMessage("Cannot assign a non-zero scalar to a Tensor<0,0,Number> object.")); - value = p.value; - return *this; -} - - -template -template -inline -bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,OtherNumber> &p) const -{ - return (value == p.value); -} - - -template -template -inline -bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,OtherNumber> &p) const -{ - return !((*this) == p); -} - - -template -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator += (const Tensor<0,dim,OtherNumber> &p) -{ - value += p.value; - return *this; -} - - -template -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,OtherNumber> &p) -{ - value -= p.value; - return *this; -} - - -template -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s) -{ - value *= s; - return *this; -} - - -template -template -inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator /= (const OtherNumber s) -{ - value /= s; - return *this; -} - - -template -inline -Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const -{ - return -value; -} - - -template -inline -typename Tensor<0,dim,Number>::real_type -Tensor<0,dim,Number>::norm () const -{ - return numbers::NumberTraits::abs (value); -} - - -template -inline -typename Tensor<0,dim,Number>::real_type -Tensor<0,dim,Number>::norm_square () const -{ - return numbers::NumberTraits::abs_square (value); -} - - -template -template -inline -void -Tensor<0, dim, Number>::unroll_recursion (Vector &result, - unsigned int &index) const -{ - result[index] = value; - ++index; -} - - -template -inline -void Tensor<0,dim,Number>::clear () -{ - value = value_type(); -} - - -template -template -inline -void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int) -{ - ar &value; -} - - -/*-------------------- Inline functions: Tensor --------------------*/ - - -namespace internal -{ - // 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] -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) - values[i] = value_type(); -} - - -template -inline -Tensor::Tensor (const Tensor &initializer) -{ - for (unsigned int i=0; i -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) - values[i] = initializer[i]; -} - - -// 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 <> -inline -Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &) -{ -} - - -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 > (values); -} - - -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 -{ - Assert (indices[0]::extract(*this, indices); -} - - -template -inline -Number & -Tensor::operator[] (const TableIndices &indices) -{ - Assert (indices[0]::extract(*this, indices); -} - - -template -inline -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> &) -{ - return *this; -} - - -template -template -inline -Tensor & -Tensor::operator = (const Tensor &t) -{ - for (unsigned int i=0; i -inline -Tensor & -Tensor::operator = (const Number d) -{ - Assert (d == Number(), 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 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 -bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const -{ - return true; -} - - -template -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) -{ - unsigned int index = 0; - for (int r = 0; r < rank_; ++r) - index = index * dim + indices[r]; - - return index; -} - - -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 &values; -} - - -#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 - -/** - * @name Vector space operations on Tensor objects: - */ -//@{ - -/** - * Multiplication of a tensor of general rank with a scalar number from the - * right. - * - * The purpose of this operator is to enable only multiplication of a tensor - * by a scalar number (i.e., a floating point number, a complex floating point - * number, etc.). The function is written in a way that only allows the - * compiler to consider the function if the second argument is indeed a scalar - * number -- in other words, @p OtherNumber will not match, for example - * std::vector@ as the product of a tensor and a vector - * clearly would make no sense. The mechanism by which the compiler is - * prohibited of considering this operator for multiplication with non-scalar - * types are explained in the documentation of the EnableIfScalar class. - * - * The return type of the function is chosen so that it matches the types of - * both the tensor and the scalar argument. For example, if you multiply a - * Tensor@<1,dim,double@> by std::complex@, - * then the result will be a - * Tensor@<1,dim,std::complex@@>. In other words, the - * type with which the returned tensor stores its components equals the type - * you would get if you multiplied an individual component of the input tensor - * by the scalar factor. - * - * @relates Tensor - * @relates EnableIfScalar - */ -template -inline -Tensor::type>::type> -operator * (const Tensor &t, - const OtherNumber factor) -{ - // recurse over the base objects - Tensor::type> tt; - for (unsigned int d=0; d -inline -Tensor::type, OtherNumber>::type> -operator * (const Number factor, - const Tensor &t) -{ - // simply forward to the operator above - return t * factor; -} - - -/** - * Division of a tensor of general rank with a scalar number. See the - * discussion on operator*() above for more information about template - * arguments and the return type. - * - * @relates Tensor - * @relates EnableIfScalar - */ -template -inline -Tensor::type>::type> -operator / (const Tensor &t, - const OtherNumber factor) -{ - // recurse over the base objects - Tensor::type> tt; - for (unsigned int d=0; d -inline -Tensor::type> -operator+ (const Tensor &p, const Tensor &q) -{ - Tensor::type> tmp (p); - - for (unsigned int i=0; i -inline -Tensor::type> -operator- (const Tensor &p, const Tensor &q) -{ - Tensor::type> tmp (p); - - for (unsigned int i=0; i