]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a Tensor<0,dim> template specialization.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 2009 03:37:44 +0000 (03:37 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 16 Jul 2009 03:37:44 +0000 (03:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@19090 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/tensor_base.h
deal.II/doc/news/changes.h
tests/base/Makefile
tests/base/scalar_01.cc [new file with mode: 0644]
tests/base/scalar_01/cmp/generic [new file with mode: 0644]

index 614f0936eb3ec7ffd9c6dfa3a7600e5d42df6cc5..7ee19d626ca02cc72b1489427a6afb9e636e91f9 100644 (file)
@@ -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 <int dim> class Point;
 // general template; specialized for rank==1; the general template is in
 // tensor.h
 template <int rank, int dim> class Tensor;
+template <int dim> class Tensor<0,dim>;
 template <int dim> class Tensor<1,dim>;
 
 
+
+/**
+ * This class is a specialized version of the
+ * <tt>Tensor<rank,dim></tt> 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@<spacedim-dim,dim@>, 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 <int dim>
+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 <tt>dimension()</tt>
+                                     * because now it is possible to
+                                     * get the dimension at compile
+                                     * time without the expansion and
+                                     * preevaluation of an inlined
+                                     * function; the compiler may
+                                     * therefore produce more
+                                     * efficient code and you may use
+                                     * this value to declare other
+                                     * data types.
+                                     */
+    static const unsigned int dimension = dim;
+
+                                    /**
+                                     * Publish the rank of this tensor to
+                                     * the outside world.
+                                     */
+    static const unsigned int rank      = 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
+                                     * <tt>factor</tt>, i.e. multiply all
+                                     * coordinates by <tt>factor</tt>.
+                                     */
+    Tensor<0,dim> & operator *= (const double factor);
+
+                                    /**
+                                     * Scale the vector by <tt>1/factor</tt>.
+                                     */
+    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 <tt>operator +=</tt> 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 <tt>operator +=</tt>
+                                     * 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
+                                      * <tt>l<sub>2</sub></tt> 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 <tt>dim</tt> 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 <tt>Tensor<rank,dim></tt> 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 <int dim>
+std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p);
+
+                                /**
+                                 *  Prints the values of this tensor in the
                                  *  form <tt>x1 x2 x3 etc</tt>.
                                  */
 template <int dim>
 std::ostream & operator << (std::ostream &out, const Tensor<1,dim> &p);
 
+
 #ifndef DOXYGEN
 
-/*------------------------------- Inline functions: Tensor ---------------------------*/
+/*---------------------------- Inline functions: Tensor<0,dim> ------------------------*/
+
+template <int dim>
+inline
+Tensor<0,dim>::Tensor ()
+{
+  Assert (dim>0, ExcDimTooSmall(dim));
+
+  value = 0;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim>::Tensor (const value_type &initializer)
+{
+  Assert (dim>0, ExcDimTooSmall(dim));
+
+  value = initializer;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim>::Tensor (const Tensor<0,dim> &p)
+{
+  Assert (dim>0, ExcDimTooSmall(dim));
+  
+  value = p.value;
+}
+
+
+
+
+template <int dim>
+inline
+Tensor<0,dim>::operator double () const
+{
+  return value;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim>::operator double & ()
+{
+  return value;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p)
+{
+  value = p.value;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator = (const double d)
+{
+  value = d;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+bool Tensor<0,dim>::operator == (const Tensor<0,dim> &p) const
+{
+  return (value == p.value);
+}
+
+
+
+template <int dim>
+inline
+bool Tensor<0,dim>::operator != (const Tensor<0,dim> &p) const
+{
+  return !((*this) == p);
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p)
+{
+  value += p.value;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p)
+{
+  value -= p.value;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s)
+{
+  value *= s;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s)
+{
+  value /= s;
+  return *this;
+}
+
+
+
+template <int dim>
+inline
+double Tensor<0,dim>::operator * (const Tensor<0,dim> &p) const
+{
+  return value*p.value;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> Tensor<0,dim>::operator + (const Tensor<0,dim> &p) const
+{
+  return value+p.value;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> Tensor<0,dim>::operator - (const Tensor<0,dim> &p) const
+{
+  return value-p.value;
+}
+
+
+
+template <int dim>
+inline
+Tensor<0,dim> Tensor<0,dim>::operator - () const
+{
+  return -value;
+}
+
+
+
+template <int dim>
+inline
+double Tensor<0,dim>::norm () const
+{
+  return std::abs (value);
+}
+
+
+
+template <int dim>
+inline
+double Tensor<0,dim>::norm_square () const
+{
+  return value*value;
+}
+
+
+
+template <int dim>
+inline
+void Tensor<0,dim>::clear ()
+{
+  value = 0;
+}
+
+
+/*---------------------------- Inline functions: Tensor<1,dim> ------------------------*/
 
 
 template <int dim>
index b0fb3efc9d96d49a772297642231ccf90dcad775..cd6585b57232f5d95fb255d9a116e10d4ebe0b4b 100644 (file)
@@ -124,6 +124,16 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li>
+  <p>
+  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.
+  <br>
+  (WB 2009/09/15)
+  </p>
+  </li>
+
   <li>
   <p>
   New: The GeometryInfo::jacobian_determinants_at_vertices can be used
index 92fcd19f85d3b1ec2e7aefbb1ffc05071b45beb5..ce67d3a7fa6cf2e31e270eca9be42df3dee72dcb 100644 (file)
@@ -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 (file)
index 0000000..272fcb2
--- /dev/null
@@ -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 <base/tensor.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <fstream>
+#include <iomanip>
+
+template <typename U, typename V>
+void compare (const U &u, const V &v)
+{
+  Assert (static_cast<double>(u)==static_cast<double>(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<double>(t1), 13.);
+  compare (static_cast<double&>(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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.