From 75bc688322886f105242bd71e0a68b071070ecd3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 16 Jul 2009 03:37:44 +0000 Subject: [PATCH] Add a Tensor<0,dim> template specialization. git-svn-id: https://svn.dealii.org/trunk@19090 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/tensor_base.h | 437 +++++++++++++++++++++++- deal.II/doc/news/changes.h | 10 + tests/base/Makefile | 3 +- tests/base/scalar_01.cc | 68 ++++ tests/base/scalar_01/cmp/generic | 2 + 5 files changed, 516 insertions(+), 4 deletions(-) create mode 100644 tests/base/scalar_01.cc create mode 100644 tests/base/scalar_01/cmp/generic diff --git a/deal.II/base/include/base/tensor_base.h b/deal.II/base/include/base/tensor_base.h index 614f0936eb..7ee19d626c 100644 --- a/deal.II/base/include/base/tensor_base.h +++ b/deal.II/base/include/base/tensor_base.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -41,9 +41,239 @@ template class Point; // general template; specialized for rank==1; the general template is in // tensor.h template class Tensor; +template class Tensor<0,dim>; template class Tensor<1,dim>; + +/** + * This class is a specialized version of the + * Tensor class. It handles tensors of rank zero, + * i.e. scalars. The second template argument 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 double) for all purposes but + * is part of the Tensor template family. + * + * @ingroup geomprimitives + * @author Wolfgang Bangerth, 2009 + */ +template +class Tensor<0,dim> +{ + public: + /** + * Provide a way to get the + * dimension of an object without + * explicit knowledge of it's + * data type. Implementation is + * this way instead of providing + * a function 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; + + /** + * Type of stored objects. This + * is a double for a rank 1 tensor. + */ + + typedef double value_type; + + /** + * Constructor. Set to zero. + */ + Tensor (); + + /** + * Copy constructor, where the + * data is copied from a C-style + * array. + */ + Tensor (const value_type &initializer); + + /** + * Copy constructor. + */ + Tensor (const Tensor<0,dim> &); + + /** + * Conversion to double. Since + * rank-0 tensors are scalars, + * this is a natural operation. + */ + operator double () const; + + /** + * Conversion to double. Since + * rank-0 tensors are scalars, + * this is a natural operation. + * + * This is the non-const + * conversion operator that + * returns a writable reference. + */ + operator double& (); + + /** + * Assignment operator. + */ + Tensor<0,dim> & operator = (const Tensor<0,dim> &); + + /** + * Assignment operator. + */ + Tensor<0,dim> & operator = (const double d); + + /** + * Test for equality of two + * tensors. + */ + bool operator == (const Tensor<0,dim> &) const; + + /** + * Test for inequality of two + * tensors. + */ + bool operator != (const Tensor<0,dim> &) const; + + /** + * Add another vector, i.e. move + * this point by the given + * offset. + */ + Tensor<0,dim> & operator += (const Tensor<0,dim> &); + + /** + * Subtract another vector. + */ + Tensor<0,dim> & operator -= (const Tensor<0,dim> &); + + /** + * Scale the vector by + * factor, i.e. multiply all + * coordinates by factor. + */ + Tensor<0,dim> & operator *= (const double factor); + + /** + * Scale the vector by 1/factor. + */ + Tensor<0,dim> & operator /= (const double factor); + + /** + * Returns the scalar product of + * two vectors. + */ + double operator * (const Tensor<0,dim> &) const; + + /** + * Add two tensors. If possible, + * use operator += instead + * since this does not need to + * copy a point at least once. + */ + Tensor<0,dim> operator + (const Tensor<0,dim> &) const; + + /** + * Subtract two tensors. If + * possible, use operator += + * instead since this does not + * need to copy a point at least + * once. + */ + Tensor<0,dim> operator - (const Tensor<0,dim> &) const; + + /** + * Tensor with inverted entries. + */ + Tensor<0,dim> 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. + */ + double norm () const; + + /** + * 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. + */ + double 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 STL and of + * several other classes within deal.II + * which not only reset the values of + * stored elements to zero, but release + * all memory and return the object into + * a virginial state. However, since the + * size of objects of the present type is + * determined by its template parameters, + * resizing is not an option, and indeed + * the state where all elements have a + * zero value is the state right after + * construction of such an object. + */ + void clear (); + + /** + * Only tensors with a positive + * dimension are implemented. This + * exception is thrown by the + * constructor if the template + * argument dim is zero or + * less. + * + * @ingroup Exceptions + */ + DeclException1 (ExcDimTooSmall, + int, + << "dim must be positive, but was " << arg1); + + private: + /** + * The value of this scalar object. + */ + double value; +}; + + /** * This class is a specialized version of the Tensor class. * It handles tensors with one index, i.e. vectors, of fixed dimension and @@ -373,15 +603,216 @@ class Tensor<1,dim> /** - * Prints the values of this point in the + * Prints the value of this scalar. + */ +template +std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &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> &p); + #ifndef DOXYGEN -/*------------------------------- Inline functions: Tensor ---------------------------*/ +/*---------------------------- Inline functions: Tensor<0,dim> ------------------------*/ + +template +inline +Tensor<0,dim>::Tensor () +{ + Assert (dim>0, ExcDimTooSmall(dim)); + + value = 0; +} + + + +template +inline +Tensor<0,dim>::Tensor (const value_type &initializer) +{ + Assert (dim>0, ExcDimTooSmall(dim)); + + value = initializer; +} + + + +template +inline +Tensor<0,dim>::Tensor (const Tensor<0,dim> &p) +{ + Assert (dim>0, ExcDimTooSmall(dim)); + + value = p.value; +} + + + + +template +inline +Tensor<0,dim>::operator double () const +{ + return value; +} + + + +template +inline +Tensor<0,dim>::operator double & () +{ + return value; +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p) +{ + value = p.value; + return *this; +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator = (const double d) +{ + value = d; + return *this; +} + + + +template +inline +bool Tensor<0,dim>::operator == (const Tensor<0,dim> &p) const +{ + return (value == p.value); +} + + + +template +inline +bool Tensor<0,dim>::operator != (const Tensor<0,dim> &p) const +{ + return !((*this) == p); +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p) +{ + value += p.value; + return *this; +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p) +{ + value -= p.value; + return *this; +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s) +{ + value *= s; + return *this; +} + + + +template +inline +Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s) +{ + value /= s; + return *this; +} + + + +template +inline +double Tensor<0,dim>::operator * (const Tensor<0,dim> &p) const +{ + return value*p.value; +} + + + +template +inline +Tensor<0,dim> Tensor<0,dim>::operator + (const Tensor<0,dim> &p) const +{ + return value+p.value; +} + + + +template +inline +Tensor<0,dim> Tensor<0,dim>::operator - (const Tensor<0,dim> &p) const +{ + return value-p.value; +} + + + +template +inline +Tensor<0,dim> Tensor<0,dim>::operator - () const +{ + return -value; +} + + + +template +inline +double Tensor<0,dim>::norm () const +{ + return std::abs (value); +} + + + +template +inline +double Tensor<0,dim>::norm_square () const +{ + return value*value; +} + + + +template +inline +void Tensor<0,dim>::clear () +{ + value = 0; +} + + +/*---------------------------- Inline functions: Tensor<1,dim> ------------------------*/ template diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index b0fb3efc9d..cd6585b572 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,16 @@ inconvenience this causes.

base

    +
  1. +

    + New: There is now a specialization Tensor<0,dim> of tensors of rank 0. Since rank-0 + tensors are scalars, this class essentially acts like a scalar, but it allows for + some neat template tricks that involve tensors of arbitrary rank. +
    + (WB 2009/09/15) +

    +
    2. +

  2. New: The GeometryInfo::jacobian_determinants_at_vertices can be used diff --git a/tests/base/Makefile b/tests/base/Makefile index 92fcd19f85..ce67d3a7fa 100644 --- a/tests/base/Makefile +++ b/tests/base/Makefile @@ -94,7 +94,8 @@ tests_x = geometry_info_* \ thread_validity_* \ task_* \ parallel_* \ - mutex_* + mutex_* \ + scalar_* # add threading tests. note that we have # to list them individually, without wildcards, because the .cc files are diff --git a/tests/base/scalar_01.cc b/tests/base/scalar_01.cc new file mode 100644 index 0000000000..272fcb2213 --- /dev/null +++ b/tests/base/scalar_01.cc @@ -0,0 +1,68 @@ +//---------------------------- scalar_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- scalar_01.cc --------------------------- + +// check Tensor<0,dim> + +#include "../tests.h" +#include +#include +#include +#include +#include + +template +void compare (const U &u, const V &v) +{ + Assert (static_cast(u)==static_cast(v), ExcInternalError()); +} + + + +int main () +{ + std::ofstream logfile("scalar_01/output"); + deallog << std::setprecision(3); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + typedef Tensor<0,1> T; + T t1(13.), t2(42); + + compare (T(), 0.); + compare (T(13.), 13.); + compare (T(t1), 13.); + compare (static_cast(t1), 13.); + compare (static_cast(t1), 13.); + compare ((T() = t1), 13.); + compare ((T() = 13.), 13.); + compare ((t1==t1), true); + compare ((t1==t2), false); + compare ((t1!=t2), true); + compare ((t1!=t1), false); + compare ((T() += t1), t1); + compare ((T() -= t1), -t1); + compare (T(13.)*=3., 39.); + compare (T(39)/=3., 13.); + compare ((t1*t2), 13*42); + compare ((t1+t2), 13+42); + compare ((t1-t2), 13-42); + compare (-t1, -13.); + compare (T(-12).norm(), 12.); + compare (T(-12).norm_square(), 12*12.); + + t1.clear(); + compare (t1, 0.); + + deallog << "OK" << std::endl; +} diff --git a/tests/base/scalar_01/cmp/generic b/tests/base/scalar_01/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/scalar_01/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5