]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Extend classes Tensor, SymmetricTensor and Point to work with templates for the numbe...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 2 Aug 2011 09:47:19 +0000 (09:47 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 2 Aug 2011 09:47:19 +0000 (09:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@23989 0785d39b-7218-0410-832d-ea1e28bc413d

23 files changed:
deal.II/doc/doxygen/headers/coding_conventions.h
deal.II/doc/news/changes.h
deal.II/include/deal.II/base/function.h
deal.II/include/deal.II/base/function_parser.h
deal.II/include/deal.II/base/point.h
deal.II/include/deal.II/base/polynomials_raviart_thomas.h
deal.II/include/deal.II/base/symmetric_tensor.h
deal.II/include/deal.II/base/tensor.h
deal.II/include/deal.II/base/tensor_base.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/include/deal.II/fe/fe_tools.h
deal.II/include/deal.II/grid/grid_in.h
deal.II/include/deal.II/grid/tria_accessor.h
deal.II/include/deal.II/lac/full_matrix.h
deal.II/include/deal.II/numerics/vectors.h
deal.II/source/base/config.cc
deal.II/source/base/function_parser.cc
deal.II/source/base/symmetric_tensor.cc
deal.II/source/base/tensor.cc
deal.II/source/fe/fe_tools.cc
deal.II/source/fe/mapping_q1.cc
deal.II/source/hp/dof_handler.cc
deal.II/source/lac/full_matrix.cc

index 7265d69103703c1cea2590bf3d546d2c0cfe03b8..9136b101b5a6018791a5cf6418615b1f9c51353b 100644 (file)
@@ -102,6 +102,10 @@ code is written, without having to look up the exact definition of something.
   and DoFHandler (named cell_iterator, active_line_iterator, etc)
   to make the connection to the STL classes clear.</li>
 
+<li> For classes with multiple template arguments, the dimension is usually
+  put before the data type specifier, i.e., we use Point<dim,number> and not
+  Point<number,dim>.
+
 <li> Each class has to have at least 200 pages of documentation ;-)</li>
 
 </ol>
index be4b64a31e7b76175856a8ad313e542f7197d782..08d6c50b3e0bf5d7369deecdda4abc29362cc18e 100644 (file)
@@ -21,6 +21,19 @@ inconvenience this causes.
 </p>
 
 <ol>
+<li> Changed: The classes Tensor, SymmetricTensor and Point now have an
+additional template argument for the number type. While a default template
+value of double ensures that all old code is still valid, this
+change invalidates forward declarations of the form <code>template
+<int dim@> class Point</code> that might be present in user-defined header
+files. Now forward declarations need to specify the type as well, i.e.,
+<code>template <int dim, typename Number> class Point</code>. However,
+nothing changes if the full declarations in <code>deal.II/base/tensor.h,
+deal.II/base/symmetric_tensor.h</code> and <code>deal.II/base/point.h</code>
+are included.
+<br>
+(Martin Kronbichler, 2011/08/02)
+
 <li> Removed: deal.II no longer supports Trilinos versions prior to 10.0.
 <br>
 (Wolfgang Bangerth, 2011/06/29)
@@ -50,7 +63,7 @@ corresponding elements of namespace std::placeholders into the global namespace
 if your compiler supported this part of the C++ 1x standard, or otherwise using
 the BOOST counterparts which are already in the global namespace. However,
 this leads to a conflict if one has a C++ 1x enabled compiler (e.g. GCC 4.6)
-<i>and</i> #includes certain BOOST headers, since the importation of symbols
+<i>and</i> includes certain BOOST headers, since the importation of symbols
 into the global namespace now leads to ambiguous names. The only solution to
 the problem is to not import names into the global namespace, but rather
 import the names from either BOOST or namespace std into the deal.II namespace
@@ -65,7 +78,18 @@ changed to use std_cxx1x::_1, std_cxx1x::_2, etc from now on.
 <h3>General</h3>
 
 <ol>
-<li> Fixed: The function VectorTools::create_right_hand_side now also works for objects of type hp::DoFHandler with different finite elements.
+<li> Extended: The classes Tensor, SymmetricTensor and Point now have an
+additional template argument for the number type. It is now possible to base
+these classes on any abstract data type that implements basic arithmetic
+operations, like <code>Tensor<1,dim,std::complex<double> ></code>. deal.II
+uses a default template argument <code>double</code> that ensures that all
+code using e.g. <code>Tensor<1,dim></code> remains valid.
+<br>
+(Martin Kronbichler, 2011/08/02)
+
+
+<li> Fixed: The function VectorTools::create_right_hand_side now also works
+for objects of type hp::DoFHandler with different finite elements.
 <br>
 (Daniel Gerecht, 2011/07/20)
 
@@ -383,12 +407,10 @@ more than <code>rhs.size()</code> CG iterations if the number of degrees of free
 is very small.
 <br>
 (Jichao Yin, WB, 2011/04/06)
-</li>
 
 <li> New: There is now a new function ConditionalOStream::get_stream().
 <br>
 (WB, 2011/03/09)
-</li>
 
 <li> Fixed: FESystem::get_unit_face_support_points would refuse to return
 anything if one of the base elements did not have support points. This
@@ -397,7 +419,6 @@ base element has no support points and also has degrees of freedom on
 the face.
 <br>
 (WB, 2011/03/07)
-</li>
 
 <li> Fixed: Objects of type FE_Nothing could be generated with multiple vector components
 by passing an argument to the constructor. Yet, the FE_Nothing::get_name() function
@@ -405,18 +426,15 @@ always just returned the string <code>FE_Nothing@<dim@>()</code> independently o
 number of components. This is now fixed.
 <br>
 (WB, 2011/03/07)
-</li>
 
 <li> Fixed: PETScWrappers:MPI:SparseMatrix and apply_boundary_values() produced an error in debug mode about non-existant SparsityPattern entries. Reason: clear_rows() also eliminated the whole row in the PETSc-internal SparsityPattern, which resulted in an error in the next assembly process.
 <br>
 (Timo Heister, 2011/02/23)
-</li>
 
 <li> Fixed: It wasn't possible to use the FE_Nothing element inside an FESystem
 object and hand the result over to an FEValues object. This is now fixed.
 <br>
 (Wolfgang Bangerth, 2011/02/18)
-</li>
 
 <li> New: There is now a function DataOutBase::write_visit_record that does
 the equivalent for VisIt that DataOutBase::write_pvtu_record does for ParaView:
@@ -424,25 +442,21 @@ generate a file that contains a list of all other VTK or VTU files of which the
 current parallel simulation consists.
 <br>
 (Wolfgang Bangerth, 2011/02/16)
-</li>
 
 <li> New: There is now a function TrilinosWrappers::VectorBase::minimal_value.
 <br>
 (Wolfgang Bangerth, 2011/02/16)
-</li>
 
 <li> Fixed: TableBase::operator= could not be compiled if the type of the
 elements of the table was <code>bool</code>. This is now fixed.
 <br>
 (Wolfgang Bangerth, 2011/02/16)
-</li>
 
 <li> Improved: The GridGenerator::hyper_shell function generated meshes in 3d
 that are valid but of poor quality upon refinement. There is now an additional
 option to generate a coarse mesh of 96 cells that has a much better quality.
 <br>
 (Wolfgang Bangerth, 2011/02/12)
-</li>
 
 <li> Fixed: There are systems where the <code>libz</code> library is installed
 but the <code>zlib.h</code> header file is not available. Since the latter
@@ -450,7 +464,6 @@ condition was not tested, this would result in compiler errors. This is now
 fixed.
 <br>
 (Wolfgang Bangerth, 2011/02/09)
-</li>
 
 <li> Fixed: Prolongation and restriction matrices were not computed at all
 for elements of type FE_DGQ if <code>dim@<spacedim</code>. Consequently,
@@ -459,7 +472,6 @@ the DoFCellAccess::set_dof_values_by_interpolation function did not
 work either and simply returned zero results. This is now fixed.
 <br>
 (M. Sebastian Pauletti, Wolfgang Bangerth, 2011/02/09)
-</li>
 
 <li> Fixed: When refining a mesh with codimension one, edges were refined using
 the same manifold description as adjacent cells, but this ignored that a
@@ -467,63 +479,58 @@ boundary indicator might have been purposefully set for edges that are truly at
 the boundary of the mesh. For such edges, the boundary indicator is now honored.
 <br>
 (M. Sebastian Pauletti, Wolfgang Bangerth, 2011/02/09)
-</li>
 
 <li> Fixed: The functions VectorTools::compute_mean_value and
 VectorTools::integrate_difference now also work for distributed
 triangulations of type parallel::distributed::Triangulation.
 <br>
 (Wolfgang Bangerth, 2011/02/07)
-</li>
 
 <li> Changed: If the <code>libz</code> library was detected during library
 configuration, the function DataOutBase::write_vtu now writes data in compressed
 format, saving a good fraction of disk space (80-90% for big output files).
 <br>
 (Wolfgang Bangerth, 2011/01/28)
-</li>
 
 <li> New: Trilinos and PETSc vectors now have a function has_ghost_elements().
 <br>
 (Timo Heister, 2011/01/26)
-</li>
 
 <li> Changed: The TrilinosWrappers::MPI::BlockVector::compress function now takes an
 argument (with a default value) in exactly the same way as the
 TrilinosWrappers::MPI::Vector::compress function already did.
 <br>
 (Wolfgang Bangerth, 2011/01/21)
-</li>
 
 <li> Fixed: When calling DataOut::build_patches with a mapping, requesting more
 than one subdivision, and when <code>dim@<spacedim</code>, then some cells
 were not properly mapped. This is now fixed.
 <br>
 (Wolfgang Bangerth, 2011/01/18)
-</li>
 
 <li> New: Restructured the internals of PETScWrappers::Precondition* to allow a
-PETSc PC object to exist without a solver. New: use Precondition*::vmult() to
+PETSc PC object to exist without a solver. New: use Precondition*::@vmult() to
 apply the preconditioner once. Preconditioners now have a default constructor
 and an initialize() function and are no longer initialized in the solver call,
 but in the constructor or initialize().
 <br>
 (Timo Heister, 2011/01/17)
-</li>
 
 <li> Fixed: Boundary conditions in the step-23 tutorial program are now
 applied correctly. Matrix columns get eliminated with the used method
 and introduce some contribution to the right hand side coming from
 inhomogeneous boundary values. The old implementation did not reset the
-matrix columns before applying new boundary values.<br>
+matrix columns before applying new boundary values.
+<br>
 (Martin Stoll, Martin Kronbichler, 2011/01/14)
-</li>
 
 <li> Extended: <code>base/tensor.h</code> has an additional collection of
 contractions between three tensors (<i>ie</i>. <code>contract3</code>).
 This can be useful for writing matrix/vector assembly in a more compact
-form than before.<br>
+form than before.
+<br>
 (Toby D. Young, 2011/01/12)
+
 </ol>
 
 
index bd2096a9c7f68f960883a12b223f205fc06cfa00..c9daa892be8e5c2c57854ef389f25c79366d9d82 100644 (file)
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/function_time.h>
 #include <deal.II/base/subscriptor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/point.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
-template <int dim> class Point;
-template <int rank_, int dim> class Tensor;
-template <int dim> class Tensor<1,dim>;
 template <typename number> class Vector;
 
 /**
@@ -109,7 +108,7 @@ class Function : public FunctionTime,
                                      * template programming.
                                      */
     static const unsigned int dimension = dim;
-    
+
                                     /**
                                      * Number of vector components.
                                      */
@@ -126,7 +125,7 @@ class Function : public FunctionTime,
                                      */
     Function (const unsigned int n_components = 1,
              const double       initial_time = 0.0);
-    
+
                                     /**
                                      * Virtual destructor; absolutely
                                      * necessary in this case.
@@ -182,7 +181,7 @@ class Function : public FunctionTime,
                                       * to.
                                       */
     Function & operator= (const Function &f);
-    
+
                                     /**
                                      * Return the value of the
                                      * function at the given
@@ -208,7 +207,7 @@ class Function : public FunctionTime,
                                      */
     virtual void vector_value (const Point<dim>   &p,
                                Vector<double>     &values) const;
-    
+
                                     /**
                                      * Set <tt>values</tt> to the point
                                      * values of the specified
@@ -245,7 +244,7 @@ class Function : public FunctionTime,
                                      */
     virtual void vector_value_list (const std::vector<Point<dim> > &points,
                                    std::vector<Vector<double> >   &values) const;
-    
+
                                     /**
                                      * For each component of the
                                      * function, fill a vector of
@@ -277,7 +276,7 @@ class Function : public FunctionTime,
                                      */
     virtual void vector_gradient (const Point<dim>            &p,
                                   std::vector<Tensor<1,dim> > &gradients) const;
-    
+
                                     /**
                                      * Set <tt>gradients</tt> to the
                                      * gradients of the specified
@@ -290,7 +289,7 @@ class Function : public FunctionTime,
     virtual void gradient_list (const std::vector<Point<dim> > &points,
                                std::vector<Tensor<1,dim> >    &gradients,
                                const unsigned int              component = 0) const;
-    
+
                                     /**
                                      * For each component of the
                                      * function, fill a vector of
@@ -306,12 +305,12 @@ class Function : public FunctionTime,
                                      */
     virtual void vector_gradients (const std::vector<Point<dim> >            &points,
                                   std::vector<std::vector<Tensor<1,dim> > > &gradients) const;
-      
+
                                     /**
                                      * Set <tt>gradients</tt> to the gradients of
                                      * the function at the <tt>points</tt>,
                                      * for all components.
-                                     * It is assumed that <tt>gradients</tt> 
+                                     * It is assumed that <tt>gradients</tt>
                                      * already has the right size, i.e.
                                      * the same size as the <tt>points</tt> array.
                                      *
@@ -338,7 +337,7 @@ class Function : public FunctionTime,
                                      */
     virtual void vector_laplacian (const Point<dim>   &p,
                                   Vector<double>     &values) const;
-    
+
                                     /**
                                      * Compute the Laplacian of one
                                      * component at a set of points.
@@ -396,13 +395,13 @@ class ZeroFunction : public Function<dim>
                                      * components is preset to one.
                                      */
     ZeroFunction (const unsigned int n_components = 1);
-    
+
                                     /**
                                      * Virtual destructor; absolutely
                                      * necessary in this case.
                                      */
     virtual ~ZeroFunction ();
-    
+
                                     /**
                                      * Return the value of the function
                                      * at the given point for one
@@ -443,7 +442,7 @@ class ZeroFunction : public Function<dim>
                                      */
     virtual void vector_value_list (const std::vector<Point<dim> > &points,
                                    std::vector<Vector<double> >   &values) const;
-    
+
                                     /**
                                      * Return the gradient of the function
                                      * at the given point, for the
@@ -460,7 +459,7 @@ class ZeroFunction : public Function<dim>
                                      */
     virtual void vector_gradient (const Point<dim>            &p,
                                   std::vector<Tensor<1,dim> > &gradients) const;
-    
+
                                     /**
                                      * Set <tt>gradients</tt> to the gradients of
                                      * the function at the <tt>points</tt>,
@@ -472,12 +471,12 @@ class ZeroFunction : public Function<dim>
     virtual void gradient_list (const std::vector<Point<dim> > &points,
                                std::vector<Tensor<1,dim> >    &gradients,
                                const unsigned int              component = 0) const;
-    
+
                                     /**
                                      * Set <tt>gradients</tt> to the gradients of
                                      * the function at the <tt>points</tt>,
                                      * for all components.
-                                     * It is assumed that <tt>gradients</tt> 
+                                     * It is assumed that <tt>gradients</tt>
                                      * already has the right size, i.e.
                                      * the same size as the <tt>points</tt> array.
                                      *
@@ -527,7 +526,7 @@ class ConstantFunction : public ZeroFunction<dim>
                                      */
     ConstantFunction (const double       value,
                      const unsigned int n_components = 1);
-    
+
                                     /**
                                      * Virtual destructor; absolutely
                                      * necessary in this case.
@@ -656,8 +655,8 @@ class ComponentSelectFunction : public ConstantFunction<dim>
                                      * number of vector components.
                                      */
     ComponentSelectFunction (const std::pair<unsigned int, unsigned int> &selected,
-                            const unsigned int n_components);    
-    
+                            const unsigned int n_components);
+
                                     /**
                                      * Return the value of the function
                                      * at the given point for all
index 745f9b347f74812240d7af7b03057fd1c1e81ce8..40317b10f17297e8122af49b651491f9ec40951e 100644 (file)
@@ -20,6 +20,8 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/function.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/point.h>
 #include <vector>
 #include <map>
 
@@ -31,9 +33,6 @@ namespace fparser
 DEAL_II_NAMESPACE_OPEN
 
 
-template <int> class Point;
-template <int, int> class Tensor;
-template <int dim> class Tensor<1,dim>;
 template <typename> class Vector;
 
 
@@ -59,23 +58,23 @@ template <typename> class Vector;
   // Define some constants that will be used by the function parser
   std::map<std::string> constants;
   constants["pi"] = numbers::PI;
-  
+
   // Define the variables that will be used inside the expressions
   std::string variables = "x,y,z";
+
   // Define the expressions of the individual components of a
   // vector valued function with two components:
   std::vector<std::string> expressions(2);
   expressions[0] = "sin(2*pi*x)+sinh(pi*z)";
   expressions[1] = "sin(2*pi*y)*exp(x^2)";
-  
+
   // Generate an empty function for these two components.
   ParsedFunction<3> vector_function(2);
+
   // And populate it with the newly created objects.
   vector_function.initialize(variables,
                            expressions,
-                           constants); 
+                           constants);
   @endverbatim
 
  * FunctionParser also provides an option to use <b>units</b> in expressions.
@@ -135,12 +134,12 @@ template <typename> class Vector;
  * For more information on this feature, please see
  * contrib/functionparser/fparser.txt
 
- * 
+
+ *
  * See http://warp.povusers.org/FunctionParser/ for an
- * explanation on how the underlying library works. 
+ * explanation on how the underlying library works.
  *
+
  From the fparser.txt file:
  @verbatim
 
@@ -182,9 +181,9 @@ template <typename> class Vector;
       example the following expression is valid: x*-y
       Note that the '=' comparison can be inaccurate due to floating point
       precision problems (eg. "sqrt(100)=10" probably returns 0, not 1).
-      
+
       The class supports these functions:
-      
+
       abs(A)    : Absolute value of A. If A is negative, returns -A otherwise
                  returns A.
       acos(A)   : Arc-cosine of A. Returns the angle, measured in radians,
@@ -241,9 +240,9 @@ template <typename> class Vector;
     "x-1"
     "-sin(sqrt(x^2+y^2))"
     "sqrt(XCoord*XCoord + YCoord*YCoord)"
-    
+
     An example of a recursive function is the factorial function:
-    
+
     "if(n>1, n*eval(n-1), 1)"
 
   Note that a recursive call has some overhead, which makes it a bit slower
@@ -263,53 +262,53 @@ template <typename> class Vector;
  *
  * An example of time dependent scalar function is the following:
       @verbatim
-     
+
       // Empty constants object
       std::map<std::string> constants;
-      
+
       // Variables that will be used inside the expressions
       std::string variables = "x,y,t";
-      
+
       // Define the expression of the scalar time dependent function.
       std::string expression = "exp(y*x)*exp(-t)";
-      
+
       // Generate an empty scalar function
       FunctionParser<2> function;
-      
+
       // And populate it with the newly created objects.
       function.initialize(variables,
                          expression,
                          constants,
-                         true);        // This tells the parser that 
+                         true);        // This tells the parser that
                                        // it is a time-dependent function
                                        // and there is another variable
                                        // to be taken into account (t).
+
      @endverbatim
-     
+
  * The following is another example of how to instantiate a
  * vector valued function by using a single string:
      @verbatim
-     
+
       // Empty constants object
       std::map<std::string> constants;
-      
+
       // Variables that will be used inside the expressions
       std::string variables = "x,y";
-      
+
       // Define the expression of the vector valued  function.
       std::string expression = "cos(2*pi*x)*y^2; sin(2*pi*x)*exp(y)";
-      
+
       // Generate an empty vector valued function
       FunctionParser<2> function(2);
-      
+
       // And populate it with the newly created objects.
       function.initialize(variables,
                          expression,
-                         constants);   
+                         constants);
+
      @endverbatim
- *    
+ *
  *
  * @ingroup functions
  * @author Luca Heltai, 2005
@@ -333,8 +332,8 @@ class FunctionParser : public Function<dim>
                                      * exception is thrown.
                                      */
     FunctionParser (const unsigned int n_components = 1,
-                   const double       initial_time = 0.0); 
-  
+                   const double       initial_time = 0.0);
+
                                     /**
                                      * Destructor. Explicitly delete
                                      * the FunctionParser objects
@@ -342,14 +341,14 @@ class FunctionParser : public Function<dim>
                                      * component of the function).
                                      */
     ~FunctionParser();
-    
+
                                     /**
                                      * Type for the constant
                                      * map. Used by the initialize()
                                      * method.
                                      */
     typedef std::map<std::string, double> ConstMap;
-    
+
                                     /**
                                      * Iterator for the constants
                                      * map. Used by the initialize()
@@ -380,7 +379,7 @@ class FunctionParser : public Function<dim>
                                       * is different from dim (if this
                                       * function is not time-dependent) or
                                       * from dim+1 (if it is time-dependent).
-                                      * 
+                                      *
                                       * <b>expressions</b>: a list of strings
                                       * containing the expressions that will
                                       * be byte compiled by the internal
@@ -404,7 +403,7 @@ class FunctionParser : public Function<dim>
                                       * name is not valid (eg:
                                       * <tt>constants["sin"] = 1.5;</tt>) an
                                       * exception is thrown.
-                                      * 
+                                      *
                                       * <b>time_dependent</b>. If this is a
                                       * time dependent function, then the last
                                       * variable declared in <b>vars</b> is
@@ -448,7 +447,7 @@ class FunctionParser : public Function<dim>
                       const ConstMap                 &units,
                       const bool time_dependent = false,
                       const bool use_degrees = false);
-  
+
                                      /**
                                       * Initialize the function. Same as
                                       * above, but accepts a string rather
@@ -521,24 +520,24 @@ class FunctionParser : public Function<dim>
                                      /** @addtogroup Exceptions
                                       * @{ */
     DeclException2 (ExcParseError,
-                    int, char*, 
+                    int, char*,
                     << "Parsing Error at Column " << arg1
                     << ". The parser said: " << arg2);
-    
+
     DeclException2 (ExcInvalidExpressionSize,
-                    int, int, 
+                    int, int,
                     << "The number of components (" << arg1
-                    << ") is not equal to the number of expressions (" 
+                    << ") is not equal to the number of expressions ("
                     << arg2 << ").");
-    
+
                                      //@}
-  private: 
+  private:
                                     /**
                                      * A pointer to the actual
                                      * function parsers.
                                      */
     fparser::FunctionParser * fp;
-    
+
                                     /**
                                      * State of usability. This
                                      * variable is checked every time
@@ -547,7 +546,7 @@ class FunctionParser : public Function<dim>
                                      * in the initialize() methods.
                                      */
     bool initialized;
-    
+
                                     /**
                                      * Number of variables. If this
                                      * is also a function of time,
@@ -584,10 +583,10 @@ FunctionParser<dim>::default_variable_names ()
   return "";
 }
 
-           
+
 
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
 
-  
+
index 600b18275cdcd36c7f57e5a751b1c9c50f93b85a..72cfc9e22eb5c4e09c94d760dd8f5dfddcb55466 100644 (file)
@@ -47,8 +47,8 @@ DEAL_II_NAMESPACE_OPEN
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 1997
  */
-template <int dim>
-class Point : public Tensor<1,dim>
+template <int dim, typename Number>
+class Point : public Tensor<1,dim,Number>
 {
   public:
                                     /**
@@ -56,7 +56,7 @@ class Point : public Tensor<1,dim>
                                      * an origin.
                                      */
     Point ();
-    
+
                                     /**
                                      * Constructor. Initialize all
                                      * entries to zero if
@@ -69,8 +69,8 @@ class Point : public Tensor<1,dim>
                                      * additional data is inside a point,
                                      * this is ok.
                                      */
-    Point (const Tensor<1,dim> &);
-    
+    Point (const Tensor<1,dim,Number> &);
+
                                     /**
                                      *  Constructor for one dimensional
                                      *  points. This function is only
@@ -78,7 +78,7 @@ class Point : public Tensor<1,dim>
                                      *  the usage is considered unsafe for
                                      *  points with <tt>dim!=1</tt>.
                                      */
-    explicit Point (const double x);
+    explicit Point (const Number x);
 
                                     /**
                                      *  Constructor for two dimensional
@@ -87,8 +87,8 @@ class Point : public Tensor<1,dim>
                                      *  the usage is considered unsafe for
                                      *  points with <tt>dim!=2</tt>.
                                      */
-    Point (const double x, const double y);
-    
+    Point (const Number x, const Number y);
+
                                     /**
                                      *  Constructor for three dimensional
                                      *  points. This function is only
@@ -96,38 +96,38 @@ class Point : public Tensor<1,dim>
                                      *  the usage is considered unsafe for
                                      *  points with <tt>dim!=3</tt>.
                                      */
-    Point (const double x, const double y, const double z);
+    Point (const Number x, const Number y, const Number z);
 
                                     /**
                                      * Return a unit vector in
                                      * coordinate direction <tt>i</tt>.
                                      */
-    static Point<dim> unit_vector(const unsigned int i);
-    
+    static Point<dim,Number> unit_vector(const unsigned int i);
+
                                     /**
                                      *  Read access to the <tt>index</tt>th
                                      *  coordinate.
                                      */
-    double   operator () (const unsigned int index) const;
+    Number   operator () (const unsigned int index) const;
 
                                     /**
                                      *  Read and write access to the
                                      *  <tt>index</tt>th coordinate.
                                      */
-    double & operator () (const unsigned int index);
+    Number & operator () (const unsigned int index);
 
 /*
  * Plus and minus operators are re-implemented from Tensor<1,dim>
  * to avoid additional casting.
  */
-                                        
+
                                     /**
                                      *  Add two point vectors. If possible,
                                      *  use <tt>operator +=</tt> instead
                                      *  since this does not need to copy a
                                      *  point at least once.
                                      */
-    Point<dim>   operator + (const Tensor<1,dim>&) const;
+    Point<dim,Number>   operator + (const Tensor<1,dim,Number>&) const;
 
                                     /**
                                      *  Subtract two point vectors. If
@@ -135,13 +135,13 @@ class Point : public Tensor<1,dim>
                                      *  instead since this does not need to
                                      *  copy a point at least once.
                                      */
-    Point<dim>   operator - (const Tensor<1,dim>&) const;
+    Point<dim,Number>   operator - (const Tensor<1,dim,Number>&) const;
 
                                     /**
                                      * The opposite vector.
                                      */
-    Point<dim>   operator - () const;
-    
+    Point<dim,Number>   operator - () const;
+
                                     /**
                                      *  Multiply by a factor. If possible,
                                      *  use <tt>operator *=</tt> instead
@@ -151,13 +151,13 @@ class Point : public Tensor<1,dim>
                                      * There is a commutative complement to this
                                      * function also
                                      */
-    Point<dim>   operator * (const double) const;
+    Point<dim,Number>   operator * (const Number) const;
 
                                     /**
                                      *  Returns the scalar product of two
                                      *  vectors.
                                      */
-    double       operator * (const Tensor<1,dim> &) const;
+    Number       operator * (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      *  Divide by a factor. If possible, use
@@ -165,15 +165,15 @@ class Point : public Tensor<1,dim>
                                      *  this does not need to copy a point at
                                      *  least once.
                                      */
-    Point<dim>   operator / (const double) const;
+    Point<dim,Number>   operator / (const Number) const;
 
                                     /**
                                      *  Returns the scalar product of this
                                      *  point vector with itself, i.e. the
                                      *  square, or the square of the norm.
                                      */
-    double              square () const;
-    
+    Number              square () const;
+
                                     /**
                                      * Returns the Euclidian distance of
                                      * <tt>this</tt> point to the point
@@ -181,12 +181,12 @@ class Point : public Tensor<1,dim>
                                      * of the difference between the vectors
                                      * representing the two points.
                                      */
-    double distance (const Point<dim> &p) const;
+    Number distance (const Point<dim,Number> &p) const;
 
                      /**
-                      * Read or write the data of this object to or 
+                      * Read or write the data of this object to or
                       * from a stream for the purpose of serialization
-                      */ 
+                      */
     template <class Archive>
     void serialize(Archive & ar, const unsigned int version);
 };
@@ -195,34 +195,34 @@ class Point : public Tensor<1,dim>
 
 #ifndef DOXYGEN
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point ()
+Point<dim,Number>::Point ()
 {}
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point (const bool initialize)
-                :
-               Tensor<1,dim>(initialize) 
+Point<dim,Number>::Point (const bool initialize)
+                          :
+                         Tensor<1,dim,Number>(initialize)
 {}
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point (const Tensor<1,dim> &t)
-                :
-               Tensor<1,dim>(t) 
+Point<dim,Number>::Point (const Tensor<1,dim,Number> &t)
+                          :
+                          Tensor<1,dim,Number>(t)
 {}
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point (const double x)
+Point<dim,Number>::Point (const Number x)
 {
   Assert (dim==1, StandardExceptions::ExcInvalidConstructorCall());
   this->values[0] = x;
@@ -230,9 +230,9 @@ Point<dim>::Point (const double x)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point (const double x, const double y)
+Point<dim,Number>::Point (const Number x, const Number y)
 {
   Assert (dim==2, StandardExceptions::ExcInvalidConstructorCall());
   this->values[0] = x;
@@ -241,9 +241,9 @@ Point<dim>::Point (const double x, const double y)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>::Point (const double x, const double y, const double z)
+Point<dim,Number>::Point (const Number x, const Number y, const Number z)
 {
   Assert (dim==3, StandardExceptions::ExcInvalidConstructorCall());
   this->values[0] = x;
@@ -252,20 +252,21 @@ Point<dim>::Point (const double x, const double y, const double z)
 }
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim>
-Point<dim>::unit_vector(unsigned int i)
+Point<dim,Number>
+Point<dim,Number>::unit_vector(unsigned int i)
 {
-  Point<dim> p;
+  Point<dim,Number> p;
   p[i] = 1.;
   return p;
 }
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Point<dim>::operator () (const unsigned int index) const
+Number
+Point<dim,Number>::operator () (const unsigned int index) const
 {
   AssertIndexRange((int) index, dim);
   return this->values[index];
@@ -273,9 +274,10 @@ double Point<dim>::operator () (const unsigned int index) const
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double & Point<dim>::operator () (const unsigned int index)
+Number &
+Point<dim,Number>::operator () (const unsigned int index)
 {
   AssertIndexRange((int) index, dim);
   return this->values[index];
@@ -283,29 +285,32 @@ double & Point<dim>::operator () (const unsigned int index)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim> Point<dim>::operator + (const Tensor<1,dim> &p) const
+Point<dim,Number>
+Point<dim,Number>::operator + (const Tensor<1,dim,Number> &p) const
 {
-  return (Point<dim>(*this) += p);
+  return (Point<dim,Number>(*this) += p);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim> Point<dim>::operator - (const Tensor<1,dim> &p) const
+Point<dim,Number>
+Point<dim,Number>::operator - (const Tensor<1,dim,Number> &p) const
 {
-  return (Point<dim>(*this) -= p);
+  return (Point<dim,Number>(*this) -= p);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim> Point<dim>::operator - () const
+Point<dim,Number>
+Point<dim,Number>::operator - () const
 {
-  Point<dim> result;
+  Point<dim,Number> result;
   for (unsigned int i=0; i<dim; ++i)
     result.values[i] = -this->values[i];
   return result;
@@ -313,67 +318,74 @@ Point<dim> Point<dim>::operator - () const
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Point<dim> Point<dim>::operator * (const double factor) const
+Point<dim,Number>
+Point<dim,Number>::operator * (const Number factor) const
 {
-  return (Point<dim>(*this) *= factor);
+  return (Point<dim,Number>(*this) *= factor);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Point<dim>::operator * (const Tensor<1,dim> &p) const
+Number
+Point<dim,Number>::operator * (const Tensor<1,dim,Number> &p) const
 {
                                   // simply pass down
-  return Tensor<1,dim>::operator * (p);
+  return Tensor<1,dim,Number>::operator * (p);
 }
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Point<dim>::square () const
+Number
+Point<dim,Number>::square () const
 {
-  double q=0;
+  Number q=0;
   for (unsigned int i=0; i<dim; ++i)
     q += this->values[i] * this->values[i];
   return q;
 }
 
 
-template <int dim>
+
+template <int dim, typename Number>
 inline
-double Point<dim>::distance (const Point<dim> &p) const
+Number
+Point<dim,Number>::distance (const Point<dim,Number> &p) const
 {
-  double sum=0;
+  Number sum=0;
   for (unsigned int i=0; i<dim; ++i)
     {
       const double diff=this->values[i]-p(i);
       sum += diff*diff;
     }
-  
+
   return std::sqrt(sum);
 }
 
 
-template <int dim>
+
+template <int dim, typename Number>
 inline
-Point<dim> Point<dim>::operator / (const double factor) const
+Point<dim,Number> Point<dim,Number>::operator / (const Number factor) const
 {
-  return (Point<dim>(*this) /= factor);
+  return (Point<dim,Number>(*this) /= factor);
 }
 
 
-template <int dim>
+
+template <int dim, typename Number>
 template <class Archive>
 inline
-void 
-Point<dim>::serialize(Archive & ar, const unsigned int)
+void
+Point<dim,Number>::serialize(Archive & ar, const unsigned int)
 {
                                      // forward to serialization
                                      // function in the base class
-  ar &  static_cast<Tensor<1,dim> &>(*this);
+  ar &  static_cast<Tensor<1,dim,Number> &>(*this);
 }
 
 #endif // DOXYGEN
@@ -382,27 +394,43 @@ Point<dim>::serialize(Archive & ar, const unsigned int)
 /*------------------------------- Global functions: Point ---------------------------*/
 
 
+/**
+ * Global operator scaling a point vector by a scalar.
+ * @relates Point
+ */
+template <int dim, typename Number>
+inline
+Point<dim,Number> operator * (const Number             factor,
+                             const Point<dim,Number> &p)
+{
+  return p*factor;
+}
+
+
+
 /**
  * Global operator scaling a point vector by a scalar.
  * @relates Point
  */
 template <int dim>
 inline
-Point<dim> operator * (const double factor, const Point<dim> &p)
+Point<dim,double> operator * (const double             factor,
+                             const Point<dim,double> &p)
 {
   return p*factor;
 }
 
 
+
 /**
  * Output operator for points. Print the elements consecutively,
  * with a space in between.
  * @relates Point
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-std::ostream & operator << (std::ostream     &out,
-                           const Point<dim> &p)
+std::ostream & operator << (std::ostream            &out,
+                           const Point<dim,Number> &p)
 {
   for (unsigned int i=0; i<dim-1; ++i)
     out << p[i] << ' ';
@@ -418,10 +446,10 @@ std::ostream & operator << (std::ostream     &out,
  * with a space in between.
  * @relates Point
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-std::istream & operator >> (std::istream &in,
-                           Point<dim>   &p)
+std::istream & operator >> (std::istream      &in,
+                           Point<dim,Number> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     in >> p[i];
@@ -432,13 +460,15 @@ std::istream & operator >> (std::istream &in,
 
 #ifndef DOXYGEN
 
-/** 
+/**
  * Output operator for points of dimension 1. This is implemented
  * specialized from the general template in order to avoid a compiler
- * warning that the loop is empty.  
+ * warning that the loop is empty.
  */
+template <typename Number>
 inline
-std::ostream & operator << (std::ostream &out, const Point<1> &p)
+std::ostream & operator << (std::ostream &out,
+                           const Point<1,Number> &p)
 {
   out << p[0];
 
index f6425982b1e62eddc45db2d5fb2f51e3a3465f1e..5e00ef9cb3936493b055b6e4292313a7671cd062 100644 (file)
@@ -59,7 +59,7 @@ class PolynomialsRaviartThomas
                                      * <i>Q<sub>k</sub></i> contained.
                                      */
     PolynomialsRaviartThomas (const unsigned int k);
-    
+
                                     /**
                                      * Computes the value and the
                                      * first and second derivatives
@@ -87,13 +87,13 @@ class PolynomialsRaviartThomas
     void compute (const Point<dim>            &unit_point,
                   std::vector<Tensor<1,dim> > &values,
                   std::vector<Tensor<2,dim> > &grads,
-                  std::vector<Tensor<3,dim> > &grad_grads) const;    
-    
+                  std::vector<Tensor<3,dim> > &grad_grads) const;
+
                                     /**
                                      * Returns the number of Raviart-Thomas polynomials.
                                      */
     unsigned int n () const;
-    
+
                                     /**
                                      * Returns the degree of the Raviart-Thomas
                                      * space, which is one less than
@@ -106,7 +106,7 @@ class PolynomialsRaviartThomas
                                      * which is <tt>RaviartThomas</tt>.
                                      */
     std::string name () const;
-    
+
                                     /**
                                      * Return the number of
                                      * polynomials in the space
@@ -117,14 +117,14 @@ class PolynomialsRaviartThomas
                                      * classes.
                                      */
     static unsigned int compute_n_pols(unsigned int degree);
-    
+
   private:
                                     /**
                                      * The degree of this object as
                                      * given to the constructor.
                                      */
     const unsigned int my_degree;
-    
+
                                     /**
                                      * An object representing the
                                      * polynomial space for a single
@@ -153,6 +153,7 @@ class PolynomialsRaviartThomas
 };
 
 
+
 template <int dim>
 inline unsigned int
 PolynomialsRaviartThomas<dim>::n() const
@@ -161,6 +162,7 @@ PolynomialsRaviartThomas<dim>::n() const
 }
 
 
+
 template <int dim>
 inline unsigned int
 PolynomialsRaviartThomas<dim>::degree() const
@@ -169,6 +171,7 @@ PolynomialsRaviartThomas<dim>::degree() const
 }
 
 
+
 template <int dim>
 inline std::string
 PolynomialsRaviartThomas<dim>::name() const
index 0ab2a688fc69d1f62e624bd84b15b6fa24c3ebb7..4029965e0b6678e023a4bb856538ad40edc57dd0 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-template <int rank, int dim> class SymmetricTensor;
+template <int rank, int dim, typename Number=double> class SymmetricTensor;
 
-template <int dim> SymmetricTensor<2,dim> unit_symmetric_tensor ();
-template <int dim> SymmetricTensor<4,dim> deviator_tensor ();
-template <int dim> SymmetricTensor<4,dim> identity_tensor ();
-template <int dim> SymmetricTensor<4,dim> invert (const SymmetricTensor<4,dim> &);
-template <int dim2> double trace (const SymmetricTensor<2,dim2> &);
+template <int dim, typename Number=double> SymmetricTensor<2,dim,Number>
+unit_symmetric_tensor ();
+template <int dim, typename Number=double> SymmetricTensor<4,dim,Number>
+deviator_tensor ();
+template <int dim, typename Number=double> SymmetricTensor<4,dim,Number>
+identity_tensor ();
+template <int dim, typename Number=double> SymmetricTensor<4,dim,Number>
+invert (const SymmetricTensor<4,dim,Number> &);
+template <int dim2, typename Number=double> Number
+trace (const SymmetricTensor<2,dim2,Number> &);
 
-template <int dim> SymmetricTensor<2,dim>
-deviator (const SymmetricTensor<2,dim> &);
-template <int dim> double
-determinant (const SymmetricTensor<2,dim> &);
+template <int dim, typename Number=double> SymmetricTensor<2,dim,Number>
+deviator (const SymmetricTensor<2,dim,Number> &);
+template <int dim, typename Number=double> Number
+determinant (const SymmetricTensor<2,dim,Number> &);
 
 
 namespace internal
@@ -116,15 +121,15 @@ namespace internal
                                      * tensor of rank
                                      * <tt>rank1+rank2-4</tt>, but if
                                      * this is zero it is a single
-                                     * scalar double. For this case,
+                                     * scalar Number. For this case,
                                      * we have a specialization.
                                      *
                                      * @author Wolfgang Bangerth, 2005
                                      */
-    template <int rank1, int rank2, int dim>
+    template <int rank1, int rank2, int dim, typename Number>
     struct double_contraction_result
     {
-       typedef ::dealii::SymmetricTensor<rank1+rank2-4,dim> type;
+      typedef ::dealii::SymmetricTensor<rank1+rank2-4,dim,Number> type;
     };
 
 
@@ -137,15 +142,15 @@ namespace internal
                                      * tensor of rank
                                      * <tt>rank1+rank2-4</tt>, but if
                                      * this is zero it is a single
-                                     * scalar double. For this case,
+                                     * scalar Number. For this case,
                                      * we have a specialization.
                                      *
                                      * @author Wolfgang Bangerth, 2005
                                      */
-    template <int dim>
-    struct double_contraction_result<2,2,dim>
+    template <int dim, typename Number>
+    struct double_contraction_result<2,2,dim,Number>
     {
-       typedef double type;
+       typedef Number type;
     };
 
 
@@ -169,15 +174,15 @@ namespace internal
                                       * about the storage format in your
                                       * application programs.
                                       */
-    template <int rank, int dim>
+    template <int rank, int dim, typename Number>
     struct StorageType;
 
                                      /**
                                       * Specialization of StorageType for
                                       * rank-2 tensors.
                                       */
-    template <int dim>
-    struct StorageType<2,dim>
+    template <int dim, typename Number>
+    struct StorageType<2,dim,Number>
     {
                                          /**
                                           * Number of independent components of a
@@ -191,7 +196,7 @@ namespace internal
                                           * Declare the type in which we actually
                                           * store the data.
                                           */
-        typedef Tensor<1,n_independent_components> base_tensor_type;
+        typedef Tensor<1,n_independent_components,Number> base_tensor_type;
     };
 
 
@@ -200,8 +205,8 @@ namespace internal
                                       * Specialization of StorageType for
                                       * rank-4 tensors.
                                       */
-    template <int dim>
-    struct StorageType<4,dim>
+    template <int dim, typename Number>
+    struct StorageType<4,dim,Number>
     {
                                          /**
                                           * Number of independent components
@@ -219,7 +224,7 @@ namespace internal
                                           */
        static const unsigned int
        n_independent_components = (n_rank2_components *
-                                   StorageType<2,dim>::n_independent_components);
+                                   StorageType<2,dim,Number>::n_independent_components);
 
                                          /**
                                           * Declare the type in which we
@@ -230,7 +235,7 @@ namespace internal
                                           * matrix if we represent the rank-2
                                           * tensors as vectors.
                                           */
-        typedef Tensor<2,n_rank2_components> base_tensor_type;
+        typedef Tensor<2,n_rank2_components,Number> base_tensor_type;
     };
 
 
@@ -241,7 +246,7 @@ namespace internal
                                       * switching on whether the tensor
                                       * should be constant or not.
                                       */
-    template <int rank, int dim, bool constness>
+    template <int rank, int dim, bool constness, typename Number>
     struct AccessorTypes;
 
                                      /**
@@ -252,12 +257,12 @@ namespace internal
                                       *
                                       * Specialization for constant tensors.
                                       */
-    template <int rank, int dim>
-    struct AccessorTypes<rank, dim,true>
+    template <int rank, int dim, typename Number>
+    struct AccessorTypes<rank,dim,true,Number>
     {
-        typedef const ::dealii::SymmetricTensor<rank,dim> tensor_type;
+        typedef const ::dealii::SymmetricTensor<rank,dim,Number> tensor_type;
 
-        typedef double reference;
+        typedef Number reference;
     };
 
                                      /**
@@ -269,12 +274,12 @@ namespace internal
                                       * Specialization for non-constant
                                       * tensors.
                                       */
-    template <int rank, int dim>
-    struct AccessorTypes<rank,dim,false>
+    template <int rank, int dim, typename Number>
+    struct AccessorTypes<rank,dim,false,Number>
     {
-        typedef ::dealii::SymmetricTensor<rank,dim> tensor_type;
+        typedef ::dealii::SymmetricTensor<rank,dim,Number> tensor_type;
 
-        typedef double &reference;
+        typedef Number &reference;
     };
 
 
@@ -313,7 +318,7 @@ namespace internal
  *
  * @author Wolfgang Bangerth, 2002, 2005
  */
-    template <int rank, int dim, bool constness, int P>
+    template <int rank, int dim, bool constness, int P, typename Number>
     class Accessor
     {
       public:
@@ -321,8 +326,8 @@ namespace internal
                                           * Import two typedefs from the
                                           * switch class above.
                                           */
-        typedef typename AccessorTypes<rank,dim,constness>::reference reference;
-        typedef typename AccessorTypes<rank,dim,constness>::tensor_type tensor_type;
+        typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
+        typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
 
       private:
                                          /**
@@ -387,7 +392,7 @@ namespace internal
                                          /**
                                           * Index operator.
                                           */
-        Accessor<rank,dim,constness,P-1> operator [] (const unsigned int i);
+        Accessor<rank,dim,constness,P-1,Number> operator [] (const unsigned int i);
 
       private:
                                          /**
@@ -402,16 +407,16 @@ namespace internal
                                          // work around bugs in some
                                          // compilers
 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
-        template <int,int> friend class SymmetricTensor;
-        template <int,int,bool,int>
+        template <int,int,typename> friend class SymmetricTensor;
+        template <int,int,bool,int,typename>
         friend class Accessor;
 #  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
-        friend class ::dealii::SymmetricTensor<rank,dim>;
-        friend class Accessor<rank,dim,constness,P+1>;
+        friend class ::dealii::SymmetricTensor<rank,dim,Number>;
+        friend class Accessor<rank,dim,constness,P+1,Number>;
 #  endif
 #else
-        friend class SymmetricTensor<rank,dim>;
-        friend class Accessor<rank,dim,constness,P+1>;
+        friend class SymmetricTensor<rank,dim,Number>;
+        friend class Accessor<rank,dim,constness,P+1,Number>;
 #endif
     };
 
@@ -427,16 +432,16 @@ namespace internal
  *
  * @author Wolfgang Bangerth, 2002, 2005
  */
-    template <int rank, int dim, bool constness>
-    class Accessor<rank,dim,constness,1>
+    template <int rank, int dim, bool constness, typename Number>
+    class Accessor<rank,dim,constness,1,Number>
     {
       public:
                                          /**
                                           * Import two typedefs from the
                                           * switch class above.
                                           */
-        typedef typename AccessorTypes<rank,dim,constness>::reference reference;
-        typedef typename AccessorTypes<rank,dim,constness>::tensor_type tensor_type;
+        typedef typename AccessorTypes<rank,dim,constness,Number>::reference reference;
+        typedef typename AccessorTypes<rank,dim,constness,Number>::tensor_type tensor_type;
 
       private:
                                          /**
@@ -521,16 +526,16 @@ namespace internal
                                          // work around bugs in some
                                          // compilers
 #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG
-        template <int,int> friend class SymmetricTensor;
-        template <int,int,bool,int>
+        template <int,int,typename> friend class SymmetricTensor;
+        template <int,int,bool,int,typename>
         friend class SymmetricTensorAccessors::Accessor;
 #  ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG
-        friend class ::dealii::SymmetricTensor<rank,dim>;
-        friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2>;
+        friend class ::dealii::SymmetricTensor<rank,dim,Number>;
+        friend class SymmetricTensorAccessors::Accessor<rank,dim,constness,2,Number>;
 #  endif
 #else
-        friend class SymmetricTensor<rank,dim>;
-        friend class Accessor<rank,dim,constness,2>;
+        friend class SymmetricTensor<rank,dim,Number>;
+        friend class Accessor<rank,dim,constness,2,Number>;
 #endif
     };
   }
@@ -604,7 +609,7 @@ namespace internal
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 2005
  */
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 class SymmetricTensor
 {
   public:
@@ -636,7 +641,7 @@ class SymmetricTensor
                                      * tensors of rank 2.
                                      */
     static const unsigned int n_independent_components
-    = internal::SymmetricTensorAccessors::StorageType<rank,dim>::
+    = internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::
       n_independent_components;
 
                                      /**
@@ -663,7 +668,7 @@ class SymmetricTensor
                                       * practice to check before calling
                                       * <tt>symmetrize</tt>.
                                       */
-    SymmetricTensor (const Tensor<2,dim> &t);
+    SymmetricTensor (const Tensor<2,dim,Number> &t);
 
                                     /**
                                      * A constructor that creates a
@@ -692,7 +697,7 @@ class SymmetricTensor
                                      * namespace is to work around
                                      * bugs in some older compilers.
                                      */
-    SymmetricTensor (const double (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim>::n_independent_components]);
+    SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::n_independent_components]);
 
                                     /**
                                      *  Assignment operator.
@@ -711,7 +716,7 @@ class SymmetricTensor
                                      * all elements of the tensor to
                                      * zero.
                                      */
-    SymmetricTensor & operator = (const double d);
+    SymmetricTensor & operator = (const Number d);
 
                                      /**
                                       * Convert the present symmetric tensor
@@ -719,7 +724,7 @@ class SymmetricTensor
                                       * elements, but using the different
                                       * storage scheme of full tensors.
                                       */
-    operator Tensor<rank,dim> () const;
+    operator Tensor<rank,dim,Number> () const;
 
                                     /**
                                      *  Test for equality of two tensors.
@@ -746,13 +751,13 @@ class SymmetricTensor
                                      *  i.e. multiply all components by
                                      *  <tt>factor</tt>.
                                      */
-    SymmetricTensor & operator *= (const double factor);
+    SymmetricTensor & operator *= (const Number factor);
 
                                     /**
                                      *  Scale the vector by
                                      *  <tt>1/factor</tt>.
                                      */
-    SymmetricTensor & operator /= (const double factor);
+    SymmetricTensor & operator /= (const Number factor);
 
                                     /**
                                      *  Add two tensors. If possible, you
@@ -827,8 +832,8 @@ class SymmetricTensor
                                      * into the first argument to the
                                      * function.
                                       */
-    typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim>::type
-    operator * (const SymmetricTensor<2,dim> &s) const;
+  typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
+  operator * (const SymmetricTensor<2,dim,Number> &s) const;
 
                                     /**
                                      * Contraction over two indices
@@ -836,14 +841,14 @@ class SymmetricTensor
                                      * rank-4 symmetric tensor given
                                      * as argument.
                                      */
-    typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim>::type
-    operator * (const SymmetricTensor<4,dim> &s) const;
+    typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
+    operator * (const SymmetricTensor<4,dim,Number> &s) const;
 
                                      /**
                                       * Return a read-write reference
                                       * to the indicated element.
                                       */
-    double & operator() (const TableIndices<rank> &indices);
+    Number & operator() (const TableIndices<rank> &indices);
 
                                      /**
                                       * Return the value of the
@@ -858,14 +863,14 @@ class SymmetricTensor
                                       * don't know here whether
                                       * copying is expensive or not.
                                       */
-    double operator() (const TableIndices<rank> &indices) const;
+    Number operator() (const TableIndices<rank> &indices) const;
 
                                      /**
                                       * Access the elements of a row of this
                                       * symmetric tensor. This function is
                                       * called for constant tensors.
                                       */
-    internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1>
+    internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
     operator [] (const unsigned int row) const;
 
                                      /**
@@ -873,7 +878,7 @@ class SymmetricTensor
                                       * symmetric tensor. This function is
                                       * called for non-constant tensors.
                                       */
-    internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1>
+    internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
     operator [] (const unsigned int row);
 
                                     /**
@@ -881,7 +886,7 @@ class SymmetricTensor
                                      * specify the entire set of
                                      * indices.
                                      */
-    double
+    Number
     operator [] (const TableIndices<rank> &indices) const;
 
                                     /**
@@ -889,7 +894,7 @@ class SymmetricTensor
                                      * specify the entire set of
                                      * indices.
                                      */
-    double &
+    Number &
     operator [] (const TableIndices<rank> &indices);
 
                                      /**
@@ -907,7 +912,7 @@ class SymmetricTensor
                                       * just one of them, although they are
                                       * equal for symmetric tensors).
                                       */
-    double norm () const;
+    Number norm () const;
 
                                     /**
                                      * Tensors can be unrolled by
@@ -980,7 +985,7 @@ class SymmetricTensor
                                      * properties of the base tensor.
                                      */
     typedef
-    internal::SymmetricTensorAccessors::StorageType<rank,dim>
+    internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>
     base_tensor_descriptor;
 
                                      /**
@@ -998,32 +1003,32 @@ class SymmetricTensor
                                     /**
                                      * Make all other symmetric tensors friends.
                                      */
-    template <int, int> friend class SymmetricTensor;
+    template <int, int, typename> friend class SymmetricTensor;
 
                                     /**
                                      * Make a few more functions friends.
                                      */
-    template <int dim2>
-    friend double trace (const SymmetricTensor<2,dim2> &d);
+    template <int dim2, typename Number2>
+    friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d);
 
-    template <int dim2>
-    friend double determinant (const SymmetricTensor<2,dim2> &t);
+    template <int dim2, typename Number2>
+    friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t);
 
-    template <int dim2>
-    friend SymmetricTensor<2,dim2>
-    deviator (const SymmetricTensor<2,dim2> &t);
+    template <int dim2, typename Number2>
+    friend SymmetricTensor<2,dim2,Number2>
+    deviator (const SymmetricTensor<2,dim2,Number2> &t);
 
-    template <int dim2>
-    friend SymmetricTensor<2,dim2> unit_symmetric_tensor ();
+    template <int dim2, typename Number2>
+    friend SymmetricTensor<2,dim2,Number2> unit_symmetric_tensor ();
 
-    template <int dim2>
-    friend SymmetricTensor<4,dim2> deviator_tensor ();
+    template <int dim2, typename Number2>
+    friend SymmetricTensor<4,dim2,Number2> deviator_tensor ();
 
-    template <int dim2>
-    friend SymmetricTensor<4,dim2> identity_tensor ();
+    template <int dim2, typename Number2>
+    friend SymmetricTensor<4,dim2,Number2> identity_tensor ();
 
-    template <int dim2>
-    friend SymmetricTensor<4,dim2> invert (const SymmetricTensor<4,dim2> &);
+    template <int dim2, typename Number2>
+    friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &);
 };
 
 
@@ -1036,8 +1041,8 @@ namespace internal
 {
   namespace SymmetricTensorAccessors
   {
-    template <int rank, int dim, bool constness, int P>
-    Accessor<rank,dim,constness,P>::
+    template <int rank, int dim, bool constness, int P, typename Number>
+    Accessor<rank,dim,constness,P,Number>::
     Accessor (tensor_type              &tensor,
              const TableIndices<rank> &previous_indices)
                    :
@@ -1047,18 +1052,18 @@ namespace internal
 
 
 
-    template <int rank, int dim, bool constness, int P>
-    Accessor<rank,dim,constness,P-1>
-    Accessor<rank,dim,constness,P>::operator[] (const unsigned int i)
+    template <int rank, int dim, bool constness, int P, typename Number>
+    Accessor<rank,dim,constness,P-1,Number>
+    Accessor<rank,dim,constness,P,Number>::operator[] (const unsigned int i)
     {
-      return Accessor<rank,dim,constness,P-1> (tensor,
-                                              merge (previous_indices, i, rank-P));
+      return Accessor<rank,dim,constness,P-1,Number> (tensor,
+                                                     merge (previous_indices, i, rank-P));
     }
 
 
 
-    template <int rank, int dim, bool constness>
-    Accessor<rank,dim,constness,1>::
+    template <int rank, int dim, bool constness, typename Number>
+    Accessor<rank,dim,constness,1,Number>::
     Accessor (tensor_type              &tensor,
              const TableIndices<rank> &previous_indices)
                    :
@@ -1068,9 +1073,9 @@ namespace internal
 
 
 
-    template <int rank, int dim, bool constness>
-    typename Accessor<rank,dim,constness,1>::reference
-    Accessor<rank,dim,constness,1>::operator[] (const unsigned int i)
+    template <int rank, int dim, bool constness, typename Number>
+    typename Accessor<rank,dim,constness,1,Number>::reference
+    Accessor<rank,dim,constness,1,Number>::operator[] (const unsigned int i)
     {
       return tensor(merge (previous_indices, i, rank-1));
     }
@@ -1081,55 +1086,61 @@ namespace internal
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>::SymmetricTensor ()
+SymmetricTensor<rank,dim,Number>::SymmetricTensor ()
 {}
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<2,2>::SymmetricTensor (const Tensor<2,2> &t)
+SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Tensor<2,dim,Number> &t)
 {
-  Assert (t[0][1] == t[1][0], ExcInternalError());
-
-  data[0] = t[0][0];
-  data[1] = t[1][1];
-  data[2] = t[0][1];
+  Assert (rank == 2, ExcNotImplemented());
+  switch (dim)
+    {
+    case 2:
+      Assert (t[0][1] == t[1][0], ExcInternalError());
+
+      data[0] = t[0][0];
+      data[1] = t[1][1];
+      data[2] = t[0][1];
+
+      break;
+    case 3:
+      Assert (t[0][1] == t[1][0], ExcInternalError());
+      Assert (t[0][2] == t[2][0], ExcInternalError());
+      Assert (t[1][2] == t[2][1], ExcInternalError());
+
+      data[0] = t[0][0];
+      data[1] = t[1][1];
+      data[2] = t[2][2];
+      data[3] = t[0][1];
+      data[4] = t[0][2];
+      data[5] = t[1][2];
+
+      break;
+    default:
+      Assert (false, ExcNotImplemented());
+    }
 }
 
 
 
-template <>
-inline
-SymmetricTensor<2,3>::SymmetricTensor (const Tensor<2,3> &t)
-{
-  Assert (t[0][1] == t[1][0], ExcInternalError());
-  Assert (t[0][2] == t[2][0], ExcInternalError());
-  Assert (t[1][2] == t[2][1], ExcInternalError());
-
-  data[0] = t[0][0];
-  data[1] = t[1][1];
-  data[2] = t[2][2];
-  data[3] = t[0][1];
-  data[4] = t[0][2];
-  data[5] = t[1][2];
-}
-
-
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>::SymmetricTensor (const double (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim>::n_independent_components])
+SymmetricTensor<rank,dim,Number>::SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::n_independent_components])
                :
                data (array)
 {}
 
 
-template <int rank, int dim>
+
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator = (const SymmetricTensor<rank,dim> &t)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator = (const SymmetricTensor<rank,dim,Number> &t)
 {
   data = t.data;
   return *this;
@@ -1137,10 +1148,10 @@ SymmetricTensor<rank,dim>::operator = (const SymmetricTensor<rank,dim> &t)
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator = (const double d)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator = (const Number d)
 {
   Assert (d==0, ExcMessage ("Only assignment with zero is allowed"));
 
@@ -1151,66 +1162,81 @@ SymmetricTensor<rank,dim>::operator = (const double d)
 
 
 
-template <>
-inline
-SymmetricTensor<2,1>::
-operator Tensor<2,1> () const
+                               // helper function to convert symmetric tensor
+                               // to generic tensor
+namespace internal
 {
-  const double t[1][1] = {{data[0]}};
-  return Tensor<2,1>(t);
-}
-
+  template <typename Number>
+  inline
+  Tensor<2,1,Number>
+  conversion (const Tensor<1,1,Number> &data)
+  {
+    const Number t[1][1] = {{data[0]}};
+    return Tensor<2,1,Number>(t);
+  }
 
+  template <typename Number>
+  inline
+  Tensor<2,2,Number>
+  conversion (const Tensor<1,3,Number> &data)
+  {
+    const Number t[2][2] = {{data[0], data[2]},
+                           {data[2], data[1]}};
+    return Tensor<2,2,Number>(t);
+  }
 
-template <>
-inline
-SymmetricTensor<2,2>::
-operator Tensor<2,2> () const
-{
-  const double t[2][2] = {{data[0], data[2]},
-                          {data[2], data[1]}};
-  return Tensor<2,2>(t);
+  template <typename Number>
+  inline
+  Tensor<2,3,Number>
+  conversion (const Tensor<1,6,Number> &data)
+  {
+    const Number t[3][3] = {{data[0], data[3], data[4]},
+                           {data[3], data[1], data[5]},
+                           {data[4], data[5], data[2]}};
+    return Tensor<2,3,Number>(t);
+  }
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<2,3>::
-operator Tensor<2,3> () const
+SymmetricTensor<rank,dim,Number>::
+operator Tensor<rank,dim,Number> () const
 {
-  const double t[3][3] = {{data[0], data[3], data[4]},
-                          {data[3], data[1], data[5]},
-                          {data[4], data[5], data[2]}};
-  return Tensor<2,3>(t);
+  Assert (rank == 2, ExcNotImplemented());
+  return internal::conversion(data);
 }
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
 bool
-SymmetricTensor<rank,dim>::operator == (const SymmetricTensor<rank,dim> &t) const
+SymmetricTensor<rank,dim,Number>::operator ==
+(const SymmetricTensor<rank,dim,Number> &t) const
 {
   return data == t.data;
 }
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
 bool
-SymmetricTensor<rank,dim>::operator != (const SymmetricTensor<rank,dim> &t) const
+SymmetricTensor<rank,dim,Number>::operator !=
+(const SymmetricTensor<rank,dim,Number> &t) const
 {
   return data != t.data;
 }
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator += (const SymmetricTensor<rank,dim> &t)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator +=
+(const SymmetricTensor<rank,dim,Number> &t)
 {
   data += t.data;
   return *this;
@@ -1218,10 +1244,11 @@ SymmetricTensor<rank,dim>::operator += (const SymmetricTensor<rank,dim> &t)
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator -= (const SymmetricTensor<rank,dim> &t)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator -=
+(const SymmetricTensor<rank,dim,Number> &t)
 {
   data -= t.data;
   return *this;
@@ -1229,10 +1256,10 @@ SymmetricTensor<rank,dim>::operator -= (const SymmetricTensor<rank,dim> &t)
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator *= (const double d)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator *= (const Number d)
 {
   data *= d;
   return *this;
@@ -1240,10 +1267,10 @@ SymmetricTensor<rank,dim>::operator *= (const double d)
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim> &
-SymmetricTensor<rank,dim>::operator /= (const double d)
+SymmetricTensor<rank,dim,Number> &
+SymmetricTensor<rank,dim,Number>::operator /= (const Number d)
 {
   data /= d;
   return *this;
@@ -1251,10 +1278,10 @@ SymmetricTensor<rank,dim>::operator /= (const double d)
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>
-SymmetricTensor<rank,dim>::operator + (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim,Number>
+SymmetricTensor<rank,dim,Number>::operator + (const SymmetricTensor &t) const
 {
   SymmetricTensor tmp = *this;
   tmp.data += t.data;
@@ -1263,10 +1290,10 @@ SymmetricTensor<rank,dim>::operator + (const SymmetricTensor &t) const
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>
-SymmetricTensor<rank,dim>::operator - (const SymmetricTensor &t) const
+SymmetricTensor<rank,dim,Number>
+SymmetricTensor<rank,dim,Number>::operator - (const SymmetricTensor &t) const
 {
   SymmetricTensor tmp = *this;
   tmp.data -= t.data;
@@ -1275,10 +1302,10 @@ SymmetricTensor<rank,dim>::operator - (const SymmetricTensor &t) const
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>
-SymmetricTensor<rank,dim>::operator - () const
+SymmetricTensor<rank,dim,Number>
+SymmetricTensor<rank,dim,Number>::operator - () const
 {
   SymmetricTensor tmp = *this;
   tmp.data = -tmp.data;
@@ -1287,417 +1314,335 @@ SymmetricTensor<rank,dim>::operator - () const
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
 void
-SymmetricTensor<rank,dim>::clear ()
+SymmetricTensor<rank,dim,Number>::clear ()
 {
   data.clear ();
 }
 
 
 
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
 std::size_t
-SymmetricTensor<rank,dim>::memory_consumption ()
+SymmetricTensor<rank,dim,Number>::memory_consumption ()
 {
   return
-    internal::SymmetricTensorAccessors::StorageType<rank,dim>::memory_consumption ();
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,2,1>::type
-SymmetricTensor<2,1>::operator * (const SymmetricTensor<2,1> &s) const
-{
-  return data[0] * s.data[0];
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,2,2>::type
-SymmetricTensor<2,2>::operator * (const SymmetricTensor<2,2> &s) const
-{
-  return (data[0] * s.data[0] +
-          data[1] * s.data[1] +
-          2*data[2] * s.data[2]);
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,2,3>::type
-SymmetricTensor<2,3>::operator * (const SymmetricTensor<2,3> &s) const
-{
-  return (data[0] * s.data[0] +
-          data[1] * s.data[1] +
-          data[2] * s.data[2] +
-          2*data[3] * s.data[3] +
-          2*data[4] * s.data[4] +
-          2*data[5] * s.data[5]);
+    internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>::memory_consumption ();
 }
 
 
 
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,2,1>::type
-SymmetricTensor<4,1>::
-operator * (const SymmetricTensor<2,1> &s) const
-{
-  const unsigned int dim = 1;
-  SymmetricTensor<2,dim> tmp;
-  tmp.data[0] = data[0][0] * s.data[0];
-  return tmp;
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,4,1>::type
-SymmetricTensor<2,1>::
-operator * (const SymmetricTensor<4,1> &s) const
-{
-  const unsigned int dim = 1;
-  SymmetricTensor<2,dim> tmp;
-  tmp[0][0] = data[0] * s.data[0][0];
-  return tmp;
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,2,2>::type
-SymmetricTensor<4,2>::
-operator * (const SymmetricTensor<2,2> &s) const
-{
-  const unsigned int dim = 2;
-  SymmetricTensor<2,dim> tmp;
-
-  for (unsigned int i=0; i<3; ++i)
-    tmp.data[i] = data[i][0] * s.data[0] +
-                 data[i][1] * s.data[1] +
-                 2 * data[i][2] * s.data[2];
-
-  return tmp;
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,4,2>::type
-SymmetricTensor<2,2>::
-operator * (const SymmetricTensor<4,2> &s) const
-{
-  const unsigned int dim = 2;
-  SymmetricTensor<2,dim> tmp;
-
-  for (unsigned int i=0; i<3; ++i)
-    tmp.data[i] = data[0] * s.data[0][i] +
-                 data[1] * s.data[1][i] +
-                 2 * data[2] * s.data[2][i];
-
-  return tmp;
-}
-
-
-
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,2,3>::type
-SymmetricTensor<4,3>::
-operator * (const SymmetricTensor<2,3> &s) const
+namespace internal
 {
-  const unsigned int dim = 3;
-  SymmetricTensor<2,dim> tmp;
 
-  for (unsigned int i=0; i<6; ++i)
-    tmp.data[i] = data[i][0] * s.data[0] +
-                 data[i][1] * s.data[1] +
-                 data[i][2] * s.data[2] +
-                 2 * data[i][3] * s.data[3] +
-                 2 * data[i][4] * s.data[4] +
-                 2 * data[i][5] * s.data[5];
+  template <int dim, typename Number>
+  inline
+  typename SymmetricTensorAccessors::double_contraction_result<2,2,dim,Number>::type
+  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
+                             const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
+  {
+    switch (dim)
+      {
+      case 1:
+       return data[0] * sdata[0];
+      case 2:
+       return (data[0] * sdata[0] +
+               data[1] * sdata[1] +
+               2*data[2] * sdata[2]);
+      case 3:
+       return (data[0] * sdata[0] +
+               data[1] * sdata[1] +
+               data[2] * sdata[2] +
+               2*data[3] * sdata[3] +
+               2*data[4] * sdata[4] +
+               2*data[5] * sdata[5]);
+      default:
+       Number sum = 0;
+       for (unsigned int d=0; d<dim; ++d)
+         sum += data[d] * sdata[d];
+       for (unsigned int d=dim; d<(dim*(dim+1)/2); ++d)
+         sum += Number(2.) * data[d] * sdata[d];
+       return sum;
+      }
+  }
 
-  return tmp;
-}
 
 
+  template <int dim, typename Number>
+  inline
+  typename SymmetricTensorAccessors::double_contraction_result<4,2,dim,Number>::type
+  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
+                             const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &sdata)
+  {
+    Number tmp [SymmetricTensorAccessors::StorageType<2,dim,Number>::n_independent_components];
+    switch (dim)
+      {
+      case 1:
+       tmp[0] = data[0][0] * sdata[0];
+       break;
+      case 2:
+       for (unsigned int i=0; i<3; ++i)
+         tmp[i] = (data[i][0] * sdata[0] +
+                   data[i][1] * sdata[1] +
+                   2 * data[i][2] * sdata[2]);
+       break;
+      case 3:
+       for (unsigned int i=0; i<6; ++i)
+         tmp[i] = (data[i][0] * sdata[0] +
+                   data[i][1] * sdata[1] +
+                   data[i][2] * sdata[2] +
+                   2 * data[i][3] * sdata[3] +
+                   2 * data[i][4] * sdata[4] +
+                   2 * data[i][5] * sdata[5]);
+       break;
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+    return SymmetricTensor<2,dim,Number>(tmp);
+  }
 
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<2,4,3>::type
-SymmetricTensor<2,3>::
-operator * (const SymmetricTensor<4,3> &s) const
-{
-  const unsigned int dim = 3;
-  SymmetricTensor<2,dim> tmp;
 
-  for (unsigned int i=0; i<6; ++i)
-    tmp.data[i] = data[0] * s.data[0][i] +
-                 data[1] * s.data[1][i] +
-                 data[2] * s.data[2][i] +
-                 2 * data[3] * s.data[3][i] +
-                 2 * data[4] * s.data[4][i] +
-                 2 * data[5] * s.data[5][i];
 
-  return tmp;
-}
+  template <int dim, typename Number>
+  inline
+  typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type
+  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data,
+                             const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
+  {
+    typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type tmp;
+    switch (dim)
+      {
+      case 1:
+       tmp[0] = data[0] * sdata[0][0];
+       break;
+      case 2:
+       for (unsigned int i=0; i<3; ++i)
+         tmp[i] = (data[0] * sdata[0][i] +
+                   data[1] * sdata[1][i] +
+                   2 * data[2] * sdata[2][i]);
+       break;
+      case 3:
+       for (unsigned int i=0; i<6; ++i)
+         tmp[i] = (data[0] * sdata[0][i] +
+                   data[1] * sdata[1][i] +
+                   data[2] * sdata[2][i] +
+                   2 * data[3] * sdata[3][i] +
+                   2 * data[4] * sdata[4][i] +
+                   2 * data[5] * sdata[5][i]);
+       break;
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+    return tmp;
+  }
 
 
 
+  template <int dim, typename Number>
+  inline
+  typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type
+  perform_double_contraction (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data,
+                             const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &sdata)
+  {
+    typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type tmp;
+    switch (dim)
+      {
+      case 1:
+       tmp[0][0] = data[0][0] * sdata[0][0];
+       break;
+      case 2:
+       for (unsigned int i=0; i<3; ++i)
+         for (unsigned int j=0; j<3; ++j)
+           tmp[i][j] = (data[i][0] * sdata[0][j] +
+                        data[i][1] * sdata[1][j] +
+                        2*data[i][2] * sdata[2][j]);
+       break;
+      case 3:
+       for (unsigned int i=0; i<6; ++i)
+         for (unsigned int j=0; j<6; ++j)
+           tmp[i][j] = (data[i][0] * sdata[0][j] +
+                        data[i][1] * sdata[1][j] +
+                        data[i][2] * sdata[2][j] +
+                        2*data[i][3] * sdata[3][j] +
+                        2*data[i][4] * sdata[4][j] +
+                        2*data[i][5] * sdata[5][j]);
+       break;
+      default:
+       Assert (false, ExcNotImplemented());
+      }
+    return tmp;
+  }
 
-template <>
-inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,4,1>::type
-SymmetricTensor<4,1>::
-operator * (const SymmetricTensor<4,1> &s) const
-{
-  const unsigned int dim = 1;
-  SymmetricTensor<4,dim> tmp;
-  tmp.data[0][0] = data[0][0] * s.data[0][0];
-  return tmp;
-}
+} // end of namespace internal
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,4,2>::type
-SymmetricTensor<4,2>::
-operator * (const SymmetricTensor<4,2> &s) const
+typename internal::SymmetricTensorAccessors::double_contraction_result<rank,2,dim,Number>::type
+SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<2,dim,Number> &s) const
 {
-  const unsigned int dim = 2;
-  SymmetricTensor<4,dim> tmp;
-  for (unsigned int i=0; i<base_tensor_descriptor::n_rank2_components; ++i)
-    for (unsigned int j=0; j<base_tensor_descriptor::n_rank2_components; ++j)
-      tmp.data[i][j] = data[i][0] * s.data[0][j] +
-                      data[i][1] * s.data[1][j] +
-                      2*data[i][2] * s.data[2][j];
-  return tmp;
+                               // need to have two different function calls
+                               // because a scalar and rank-2 tensor are not
+                               // the same data type (see internal function
+                               // above)
+  return internal::perform_double_contraction<dim,Number> (data, s.data);
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-internal::SymmetricTensorAccessors::double_contraction_result<4,4,3>::type
-SymmetricTensor<4,3>::
-operator * (const SymmetricTensor<4,3> &s) const
+typename internal::SymmetricTensorAccessors::double_contraction_result<rank,4,dim,Number>::type
+SymmetricTensor<rank,dim,Number>::operator * (const SymmetricTensor<4,dim,Number> &s) const
 {
-  const unsigned int dim = 3;
-  SymmetricTensor<4,dim> tmp;
-  for (unsigned int i=0; i<base_tensor_descriptor::n_rank2_components; ++i)
-    for (unsigned int j=0; j<base_tensor_descriptor::n_rank2_components; ++j)
-      tmp.data[i][j] = data[i][0] * s.data[0][j] +
-                      data[i][1] * s.data[1][j] +
-                      data[i][2] * s.data[2][j] +
-                      2*data[i][3] * s.data[3][j] +
-                      2*data[i][4] * s.data[4][j] +
-                      2*data[i][5] * s.data[5][j];
+  typename internal::SymmetricTensorAccessors::
+    double_contraction_result<rank,4,dim,Number>::type tmp;
+  tmp.data = internal::perform_double_contraction<dim,Number> (data,s.data);
   return tmp;
 }
 
 
 
-template <>
-inline
-double &
-SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices)
-{
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < 1, ExcIndexRange (indices[r], 0, 1));
-
-  return data[0];
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices) const
-{
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
-
-  return data[0];
-}
-
-
-
-template <>
-inline
-double &
-SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices)
+                               // internal namespace to switch between the
+                               // access of different tensors. There used to
+                               // be explicit instantiations before for
+                               // different ranks and dimensions, but since
+                               // we now allow for templates on the data
+                               // type, and since we cannot partially
+                               // specialize the implementation, this got
+                               // into a separate namespace
+namespace internal
 {
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+  template <int dim, typename Number>
+  inline
+  Number &
+  symmetric_tensor_access (const TableIndices<2> &indices,
+                          typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
+  {
+    switch (dim)
+      {
+      case 1:
+       return data[0];
 
+      case 2:
                                    // first treat the main diagonal
                                    // elements, which are stored
                                    // consecutively at the beginning
-  if (indices[0] == indices[1])
-    return data[indices[0]];
+       if (indices[0] == indices[1])
+         return data[indices[0]];
 
                                    // the rest is messier and requires a few
                                    // switches. at least for the 2x2 case it
                                    // is reasonably simple
-  Assert (((indices[0]==1) && (indices[1]==0)) ||
-          ((indices[0]==0) && (indices[1]==1)),
-          ExcInternalError());
-  return data[2];
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<2,2>::operator () (const TableIndices<2> &indices) const
-{
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+       Assert (((indices[0]==1) && (indices[1]==0)) ||
+               ((indices[0]==0) && (indices[1]==1)),
+               ExcInternalError());
+       return data[2];
 
+      case 3:
                                    // first treat the main diagonal
                                    // elements, which are stored
                                    // consecutively at the beginning
-  if (indices[0] == indices[1])
-    return data[indices[0]];
+       if (indices[0] == indices[1])
+         return data[indices[0]];
 
                                    // the rest is messier and requires a few
-                                   // switches. at least for the 2x2 case it
-                                   // is reasonably simple
-  Assert (((indices[0]==1) && (indices[1]==0)) ||
-          ((indices[0]==0) && (indices[1]==1)),
-          ExcInternalError());
-  return data[2];
-}
+                                   // switches, but simpler if we just sort
+                                   // our indices
+       {
+         TableIndices<2> sorted_indices (indices);
+         sorted_indices.sort ();
+
+         if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
+           return data[3];
+         else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
+           return data[4];
+         else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
+           return data[5];
+         else
+           Assert (false, ExcInternalError());
+       }
+      }
 
+    static Number dummy_but_referenceable = 0;
+    return dummy_but_referenceable;
+  }
 
 
-template <>
-inline
-double &
-SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices)
-{
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
 
+  template <int dim, typename Number>
+  inline
+  Number
+  symmetric_tensor_access (const TableIndices<2> &indices,
+                          const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
+  {
+    switch (dim)
+      {
+      case 1:
+       return data[0];
+
+      case 2:
                                    // first treat the main diagonal
                                    // elements, which are stored
                                    // consecutively at the beginning
-  if (indices[0] == indices[1])
-    return data[indices[0]];
+       if (indices[0] == indices[1])
+         return data[indices[0]];
 
                                    // the rest is messier and requires a few
-                                   // switches, but simpler if we just sort
-                                   // our indices
-  TableIndices<2> sorted_indices (indices);
-  sorted_indices.sort ();
-
-  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
-    return data[3];
-  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
-    return data[4];
-  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
-    return data[5];
-  else
-    Assert (false, ExcInternalError());
-
-  static double dummy_but_referenceable = 0;
-  return dummy_but_referenceable;
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<2,3>::operator () (const TableIndices<2> &indices) const
-{
-  const unsigned int rank = 2;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+                                   // switches. at least for the 2x2 case it
+                                   // is reasonably simple
+       Assert (((indices[0]==1) && (indices[1]==0)) ||
+               ((indices[0]==0) && (indices[1]==1)),
+               ExcInternalError());
+       return data[2];
 
+      case 3:
                                    // first treat the main diagonal
                                    // elements, which are stored
                                    // consecutively at the beginning
-  if (indices[0] == indices[1])
-    return data[indices[0]];
+       if (indices[0] == indices[1])
+         return data[indices[0]];
 
                                    // the rest is messier and requires a few
                                    // switches, but simpler if we just sort
                                    // our indices
-  TableIndices<2> sorted_indices (indices);
-  sorted_indices.sort ();
-
-  if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
-    return data[3];
-  else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
-    return data[4];
-  else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
-    return data[5];
-  else
-    Assert (false, ExcInternalError());
-
-  static double dummy_but_referenceable = 0;
-  return dummy_but_referenceable;
-}
-
-
-
-template <>
-inline
-double &
-SymmetricTensor<4,1>::operator () (const TableIndices<4> &indices)
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
-
-  return data[0][0];
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<4,1>::operator () (const TableIndices<4> &indices) const
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+       {
+         TableIndices<2> sorted_indices (indices);
+         sorted_indices.sort ();
+
+         if ((sorted_indices[0]==0) && (sorted_indices[1]==1))
+           return data[3];
+         else if ((sorted_indices[0]==0) && (sorted_indices[1]==2))
+           return data[4];
+         else if ((sorted_indices[0]==1) && (sorted_indices[1]==2))
+           return data[5];
+         else
+           Assert (false, ExcInternalError());
+       }
+      }
 
-  return data[0][0];
-}
+    static Number dummy_but_referenceable = 0;
+    return dummy_but_referenceable;
+  }
 
 
 
-template <>
-inline
-double &
-SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices)
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+  template <int dim, typename Number>
+  inline
+  Number &
+  symmetric_tensor_access (const TableIndices<4> &indices,
+                          typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
+  {
+    switch (dim)
+      {
+      case 1:
+       return data[0][0];
 
+      case 2:
                                   // each entry of the tensor can be
                                   // thought of as an entry in a
                                   // matrix that maps the rolled-out
@@ -1707,35 +1652,26 @@ SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices)
                                   // tensors. determine which
                                   // position the present entry is
                                   // stored in
-  unsigned int base_index[2] ;
-  if ((indices[0] == 0) && (indices[1] == 0))
-    base_index[0] = 0;
-  else if ((indices[0] == 1) && (indices[1] == 1))
-    base_index[0] = 1;
-  else
-    base_index[0] = 2;
-
-  if ((indices[2] == 0) && (indices[3] == 0))
-    base_index[1] = 0;
-  else if ((indices[2] == 1) && (indices[3] == 1))
-    base_index[1] = 1;
-  else
-    base_index[1] = 2;
-
-  return data[base_index[0]][base_index[1]];
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices) const
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+       {
+         unsigned int base_index[2] ;
+         if ((indices[0] == 0) && (indices[1] == 0))
+           base_index[0] = 0;
+         else if ((indices[0] == 1) && (indices[1] == 1))
+           base_index[0] = 1;
+         else
+           base_index[0] = 2;
+
+         if ((indices[2] == 0) && (indices[3] == 0))
+           base_index[1] = 0;
+         else if ((indices[2] == 1) && (indices[3] == 1))
+           base_index[1] = 1;
+         else
+           base_index[1] = 2;
+
+         return data[base_index[0]][base_index[1]];
+       }
 
+      case 3:
                                   // each entry of the tensor can be
                                   // thought of as an entry in a
                                   // matrix that maps the rolled-out
@@ -1745,35 +1681,66 @@ SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices) const
                                   // tensors. determine which
                                   // position the present entry is
                                   // stored in
-  unsigned int base_index[2] ;
-  if ((indices[0] == 0) && (indices[1] == 0))
-    base_index[0] = 0;
-  else if ((indices[0] == 1) && (indices[1] == 1))
-    base_index[0] = 1;
-  else
-    base_index[0] = 2;
-
-  if ((indices[2] == 0) && (indices[3] == 0))
-    base_index[1] = 0;
-  else if ((indices[2] == 1) && (indices[3] == 1))
-    base_index[1] = 1;
-  else
-    base_index[1] = 2;
-
-  return data[base_index[0]][base_index[1]];
-}
-
+       {
+         unsigned int base_index[2] ;
+         if ((indices[0] == 0) && (indices[1] == 0))
+           base_index[0] = 0;
+         else if ((indices[0] == 1) && (indices[1] == 1))
+           base_index[0] = 1;
+         else if ((indices[0] == 2) && (indices[1] == 2))
+           base_index[0] = 2;
+         else if (((indices[0] == 0) && (indices[1] == 1)) ||
+                  ((indices[0] == 1) && (indices[1] == 0)))
+           base_index[0] = 3;
+         else if (((indices[0] == 0) && (indices[1] == 2)) ||
+                  ((indices[0] == 2) && (indices[1] == 0)))
+           base_index[0] = 4;
+         else
+           {
+             Assert (((indices[0] == 1) && (indices[1] == 2)) ||
+                     ((indices[0] == 2) && (indices[1] == 1)),
+                     ExcInternalError());
+             base_index[0] = 5;
+           }
+
+         if ((indices[2] == 0) && (indices[3] == 0))
+           base_index[1] = 0;
+         else if ((indices[2] == 1) && (indices[3] == 1))
+           base_index[1] = 1;
+         else if ((indices[2] == 2) && (indices[3] == 2))
+           base_index[1] = 2;
+         else if (((indices[2] == 0) && (indices[3] == 1)) ||
+                  ((indices[2] == 1) && (indices[3] == 0)))
+           base_index[1] = 3;
+         else if (((indices[2] == 0) && (indices[3] == 2)) ||
+                  ((indices[2] == 2) && (indices[3] == 0)))
+           base_index[1] = 4;
+         else
+           {
+             Assert (((indices[2] == 1) && (indices[3] == 2)) ||
+                     ((indices[2] == 2) && (indices[3] == 1)),
+                     ExcInternalError());
+             base_index[1] = 5;
+           }
+
+         return data[base_index[0]][base_index[1]];
+       }
+      }
+  }
 
 
-template <>
-inline
-double &
-SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices)
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+  template <int dim, typename Number>
+  inline
+  Number
+  symmetric_tensor_access (const TableIndices<4> &indices,
+                          const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
+  {
+    switch (dim)
+      {
+      case 1:
+       return data[0][0];
 
+      case 2:
                                   // each entry of the tensor can be
                                   // thought of as an entry in a
                                   // matrix that maps the rolled-out
@@ -1783,61 +1750,26 @@ SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices)
                                   // tensors. determine which
                                   // position the present entry is
                                   // stored in
-  unsigned int base_index[2] ;
-  if ((indices[0] == 0) && (indices[1] == 0))
-    base_index[0] = 0;
-  else if ((indices[0] == 1) && (indices[1] == 1))
-    base_index[0] = 1;
-  else if ((indices[0] == 2) && (indices[1] == 2))
-    base_index[0] = 2;
-  else if (((indices[0] == 0) && (indices[1] == 1)) ||
-          ((indices[0] == 1) && (indices[1] == 0)))
-    base_index[0] = 3;
-  else if (((indices[0] == 0) && (indices[1] == 2)) ||
-          ((indices[0] == 2) && (indices[1] == 0)))
-    base_index[0] = 4;
-  else
-    {
-      Assert (((indices[0] == 1) && (indices[1] == 2)) ||
-             ((indices[0] == 2) && (indices[1] == 1)),
-             ExcInternalError());
-      base_index[0] = 5;
-    }
-
-  if ((indices[2] == 0) && (indices[3] == 0))
-    base_index[1] = 0;
-  else if ((indices[2] == 1) && (indices[3] == 1))
-    base_index[1] = 1;
-  else if ((indices[2] == 2) && (indices[3] == 2))
-    base_index[1] = 2;
-  else if (((indices[2] == 0) && (indices[3] == 1)) ||
-          ((indices[2] == 1) && (indices[3] == 0)))
-    base_index[1] = 3;
-  else if (((indices[2] == 0) && (indices[3] == 2)) ||
-          ((indices[2] == 2) && (indices[3] == 0)))
-    base_index[1] = 4;
-  else
-    {
-      Assert (((indices[2] == 1) && (indices[3] == 2)) ||
-             ((indices[2] == 2) && (indices[3] == 1)),
-             ExcInternalError());
-      base_index[1] = 5;
-    }
-
-  return data[base_index[0]][base_index[1]];
-}
-
-
-
-template <>
-inline
-double
-SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices) const
-{
-  const unsigned int rank = 4;
-  for (unsigned int r=0; r<rank; ++r)
-    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+       {
+         unsigned int base_index[2] ;
+         if ((indices[0] == 0) && (indices[1] == 0))
+           base_index[0] = 0;
+         else if ((indices[0] == 1) && (indices[1] == 1))
+           base_index[0] = 1;
+         else
+           base_index[0] = 2;
+
+         if ((indices[2] == 0) && (indices[3] == 0))
+           base_index[1] = 0;
+         else if ((indices[2] == 1) && (indices[3] == 1))
+           base_index[1] = 1;
+         else
+           base_index[1] = 2;
+
+         return data[base_index[0]][base_index[1]];
+       }
 
+      case 3:
                                   // each entry of the tensor can be
                                   // thought of as an entry in a
                                   // matrix that maps the rolled-out
@@ -1847,305 +1779,293 @@ SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices) const
                                   // tensors. determine which
                                   // position the present entry is
                                   // stored in
-  unsigned int base_index[2] ;
-  if ((indices[0] == 0) && (indices[1] == 0))
-    base_index[0] = 0;
-  else if ((indices[0] == 1) && (indices[1] == 1))
-    base_index[0] = 1;
-  else if ((indices[0] == 2) && (indices[1] == 2))
-    base_index[0] = 2;
-  else if (((indices[0] == 0) && (indices[1] == 1)) ||
-          ((indices[0] == 1) && (indices[1] == 0)))
-    base_index[0] = 3;
-  else if (((indices[0] == 0) && (indices[1] == 2)) ||
-          ((indices[0] == 2) && (indices[1] == 0)))
-    base_index[0] = 4;
-  else
-    {
-      Assert (((indices[0] == 1) && (indices[1] == 2)) ||
-             ((indices[0] == 2) && (indices[1] == 1)),
-             ExcInternalError());
-      base_index[0] = 5;
-    }
-
-  if ((indices[2] == 0) && (indices[3] == 0))
-    base_index[1] = 0;
-  else if ((indices[2] == 1) && (indices[3] == 1))
-    base_index[1] = 1;
-  else if ((indices[2] == 2) && (indices[3] == 2))
-    base_index[1] = 2;
-  else if (((indices[2] == 0) && (indices[3] == 1)) ||
-          ((indices[2] == 1) && (indices[3] == 0)))
-    base_index[1] = 3;
-  else if (((indices[2] == 0) && (indices[3] == 2)) ||
-          ((indices[2] == 2) && (indices[3] == 0)))
-    base_index[1] = 4;
-  else
-    {
-      Assert (((indices[2] == 1) && (indices[3] == 2)) ||
-             ((indices[2] == 2) && (indices[3] == 1)),
-             ExcInternalError());
-      base_index[1] = 5;
-    }
-
-  return data[base_index[0]][base_index[1]];
-}
-
-
-
-template <int rank, int dim>
-internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1>
-SymmetricTensor<rank,dim>::operator [] (const unsigned int row) const
-{
-  return
-    internal::SymmetricTensorAccessors::
-    Accessor<rank,dim,true,rank-1> (*this, TableIndices<rank> (row));
-}
-
-
-
-template <int rank, int dim>
-internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1>
-SymmetricTensor<rank,dim>::operator [] (const unsigned int row)
-{
-  return
-    internal::SymmetricTensorAccessors::
-    Accessor<rank,dim,false,rank-1> (*this, TableIndices<rank> (row));
-}
-
-
-
-template <int rank, int dim>
-inline
-double
-SymmetricTensor<rank,dim>::operator [] (const TableIndices<rank> &indices) const
-{
-  return data[component_to_unrolled_index(indices)];
-}
-
-
+       {
+         unsigned int base_index[2] ;
+         if ((indices[0] == 0) && (indices[1] == 0))
+           base_index[0] = 0;
+         else if ((indices[0] == 1) && (indices[1] == 1))
+           base_index[0] = 1;
+         else if ((indices[0] == 2) && (indices[1] == 2))
+           base_index[0] = 2;
+         else if (((indices[0] == 0) && (indices[1] == 1)) ||
+                  ((indices[0] == 1) && (indices[1] == 0)))
+           base_index[0] = 3;
+         else if (((indices[0] == 0) && (indices[1] == 2)) ||
+                  ((indices[0] == 2) && (indices[1] == 0)))
+           base_index[0] = 4;
+         else
+           {
+             Assert (((indices[0] == 1) && (indices[1] == 2)) ||
+                     ((indices[0] == 2) && (indices[1] == 1)),
+                     ExcInternalError());
+             base_index[0] = 5;
+           }
+
+         if ((indices[2] == 0) && (indices[3] == 0))
+           base_index[1] = 0;
+         else if ((indices[2] == 1) && (indices[3] == 1))
+           base_index[1] = 1;
+         else if ((indices[2] == 2) && (indices[3] == 2))
+           base_index[1] = 2;
+         else if (((indices[2] == 0) && (indices[3] == 1)) ||
+                  ((indices[2] == 1) && (indices[3] == 0)))
+           base_index[1] = 3;
+         else if (((indices[2] == 0) && (indices[3] == 2)) ||
+                  ((indices[2] == 2) && (indices[3] == 0)))
+           base_index[1] = 4;
+         else
+           {
+             Assert (((indices[2] == 1) && (indices[3] == 2)) ||
+                     ((indices[2] == 2) && (indices[3] == 1)),
+                     ExcInternalError());
+             base_index[1] = 5;
+           }
+
+         return data[base_index[0]][base_index[1]];
+       }
+      }
+  }
 
-template <int rank, int dim>
-inline
-double &
-SymmetricTensor<rank,dim>::operator [] (const TableIndices<rank> &indices)
-{
-  return data[component_to_unrolled_index(indices)];
-}
+} // end of namespace internal
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-double
-SymmetricTensor<2,1>::norm () const
+Number &
+SymmetricTensor<rank,dim,Number>::operator () (const TableIndices<rank> &indices)
 {
-  return std::fabs(data[0]);
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+  return internal::symmetric_tensor_access<dim,Number> (indices, data);
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-double
-SymmetricTensor<2,2>::norm () const
+Number
+SymmetricTensor<rank,dim,Number>::operator ()
+  (const TableIndices<rank> &indices) const
 {
-  return std::sqrt(data[0]*data[0] + data[1]*data[1] + 2*data[2]*data[2]);
+  for (unsigned int r=0; r<rank; ++r)
+    Assert (indices[r] < dimension, ExcIndexRange (indices[r], 0, dimension));
+  return internal::symmetric_tensor_access<dim,Number> (indices, data);
 }
 
 
 
-template <>
-inline
-double
-SymmetricTensor<2,3>::norm () const
+template <int rank, int dim, typename Number>
+internal::SymmetricTensorAccessors::Accessor<rank,dim,true,rank-1,Number>
+SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row) const
 {
-  return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2] +
-                   2*data[3]*data[3] + 2*data[4]*data[4] + 2*data[5]*data[5]);
+  return
+    internal::SymmetricTensorAccessors::
+    Accessor<rank,dim,true,rank-1,Number> (*this, TableIndices<rank> (row));
 }
 
 
 
-template <>
-inline
-double
-SymmetricTensor<4,1>::norm () const
+template <int rank, int dim, typename Number>
+internal::SymmetricTensorAccessors::Accessor<rank,dim,false,rank-1,Number>
+SymmetricTensor<rank,dim,Number>::operator [] (const unsigned int row)
 {
-  return std::fabs(data[0][0]);
+  return
+    internal::SymmetricTensorAccessors::
+    Accessor<rank,dim,false,rank-1,Number> (*this, TableIndices<rank> (row));
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-double
-SymmetricTensor<4,2>::norm () const
+Number
+SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices) const
 {
-  const unsigned int dim = 2;
-
-                                  // this is not really efficient and
-                                  // could be improved by counting
-                                  // how often each tensor entry is
-                                  // accessed, but this isn't a
-                                  // really frequent operation anyway
-  double t = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=0; j<dim; ++j)
-      for (unsigned int k=0; k<dim; ++k)
-       for (unsigned int l=0; l<dim; ++l)
-         {
-           const double a = (*this)[i][j][k][l];
-           t += a * a;
-         }
-  return std::sqrt(t);
+  return data[component_to_unrolled_index(indices)];
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-double
-SymmetricTensor<4,3>::norm () const
+Number &
+SymmetricTensor<rank,dim,Number>::operator [] (const TableIndices<rank> &indices)
 {
-  const unsigned int dim = 3;
-
-                                  // this is not really efficient and
-                                  // could be improved by counting
-                                  // how often each tensor entry is
-                                  // accessed, but this isn't a
-                                  // really frequent operation anyway
-  double t = 0;
-  for (unsigned int i=0; i<dim; ++i)
-    for (unsigned int j=0; j<dim; ++j)
-      for (unsigned int k=0; k<dim; ++k)
-       for (unsigned int l=0; l<dim; ++l)
-         {
-           const double a = (*this)[i][j][k][l];
-           t += a * a;
-         }
-  return std::sqrt(t);
+  return data[component_to_unrolled_index(indices)];
 }
 
 
 
-template <>
-inline
-unsigned int
-SymmetricTensor<2,1>::component_to_unrolled_index (const TableIndices<2> &indices)
+namespace internal
 {
-  Assert (indices[0] < dimension, ExcIndexRange(indices[0], 0, dimension));
-  Assert (indices[1] < dimension, ExcIndexRange(indices[1], 0, dimension));
-
-  return 0;
-}
-
+  template <int dim, typename Number>
+  inline
+  Number
+  compute_norm (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data)
+  {
+    Number return_value;
+    switch (dim)
+      {
+      case 1:
+       return_value = std::fabs(data[0]);
+       break;
+      case 2:
+       return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] +
+                                2*data[2]*data[2]);
+       break;
+      case 3:
+       return_value =  std::sqrt(data[0]*data[0] + data[1]*data[1] +
+                                 data[2]*data[2] + 2*data[3]*data[3] +
+                                 2*data[4]*data[4] + 2*data[5]*data[5]);
+       break;
+      default:
+       return_value = 0;
+       for (unsigned int d=0; d<dim; ++d)
+         return_value += data[d] * data[d];
+       for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d)
+         return_value += 2 * data[d] * data[d];
+       return_value = std::sqrt(return_value);
+      }
+    return return_value;
+  }
 
 
-template <>
-inline
-unsigned int
-SymmetricTensor<2,2>::component_to_unrolled_index (const TableIndices<2> &indices)
-{
-  Assert (indices[0] < dimension, ExcIndexRange(indices[0], 0, dimension));
-  Assert (indices[1] < dimension, ExcIndexRange(indices[1], 0, dimension));
 
-  static const unsigned int table[dimension][dimension] = {{0, 2},
-                                                          {2, 1}};
+  template <int dim, typename Number>
+  inline
+  Number
+  compute_norm (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data)
+  {
+    Number return_value;
+    const unsigned int n_independent_components = data.dimension;
+
+    switch (dim)
+      {
+      case 1:
+       return_value = std::fabs (data[0][0]);
+       break;
+      default:
+       return_value = 0;
+       for (unsigned int i=0; i<dim; ++i)
+         for (unsigned int j=0; j<dim; ++j)
+           return_value += data[i][j] * data[i][j];
+       for (unsigned int i=0; i<dim; ++i)
+         for (unsigned int j=dim; j<n_independent_components; ++j)
+           return_value += 2 * data[i][j] * data[i][j];
+       for (unsigned int i=dim; i<n_independent_components; ++i)
+         for (unsigned int j=0; j<dim; ++j)
+           return_value += 2 * data[i][j] * data[i][j];
+       for (unsigned int i=dim; i<n_independent_components; ++i)
+         for (unsigned int j=dim; j<n_independent_components; ++j)
+           return_value += 4 * data[i][j] * data[i][j];
+       return_value = std::sqrt(return_value);
+      }
+
+    return return_value;
+  }
 
-  return table[indices[0]][indices[1]];
-}
+} // end of namespace internal
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-unsigned int
-SymmetricTensor<2,3>::component_to_unrolled_index (const TableIndices<2> &indices)
+Number
+SymmetricTensor<rank,dim,Number>::norm () const
 {
-  const unsigned int dim = 3;
-  Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
-  Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
-
-  static const unsigned int table[dim][dim] = {{0, 3, 4},
-                                              {3, 1, 5},
-                                              {4, 5, 2}};
-
-  return table[indices[0]][indices[1]];
+  return internal::compute_norm<dim,Number> (data);
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
 unsigned int
-SymmetricTensor<2,4>::component_to_unrolled_index (const TableIndices<2> &indices)
+SymmetricTensor<rank,dim,Number>::component_to_unrolled_index
+  (const TableIndices<rank> &indices)
 {
-  const unsigned int dim = 4;
+  Assert (rank == 2, ExcNotImplemented());
   Assert (indices[0] < dim, ExcIndexRange(indices[0], 0, dim));
   Assert (indices[1] < dim, ExcIndexRange(indices[1], 0, dim));
 
-  static const unsigned int table[dim][dim] = {{0, 4, 5, 6},
-                                              {4, 1, 7, 8},
-                                              {5, 7, 2, 9},
-                                              {6, 8, 9, 3}};
-
-  return table[indices[0]][indices[1]];
-}
-
-
-
-template <>
-inline
-TableIndices<2>
-SymmetricTensor<2,1>::unrolled_to_component_indices (const unsigned int i)
-{
-  Assert (i < n_independent_components, ExcIndexRange(i, 0, n_independent_components));
-
-  return TableIndices<2>(0,0);
+  switch(dim)
+    {
+    case 1:
+      return 0;
+    case 2:
+      {
+       static const unsigned int table[2][2] = {{0, 2},
+                                                {2, 1}};
+       return table[indices[0]][indices[1]];
+      }
+    case 3:
+      {
+       static const unsigned int table[3][3] = {{0, 3, 4},
+                                                {3, 1, 5},
+                                                {4, 5, 2}};
+       return table[indices[0]][indices[1]];
+      }
+    case 4:
+      {
+       static const unsigned int table[4][4] = {{0, 4, 5, 6},
+                                                {4, 1, 7, 8},
+                                                {5, 7, 2, 9},
+                                                {6, 8, 9, 3}};
+       return table[indices[0]][indices[1]];
+      }
+    default:
+      Assert (false, ExcNotImplemented());
+      return 0;
+    }
 }
 
 
 
-template <>
+template <int rank, int dim, typename Number>
 inline
-TableIndices<2>
-SymmetricTensor<2,2>::unrolled_to_component_indices (const unsigned int i)
+TableIndices<rank>
+SymmetricTensor<rank,dim,Number>::unrolled_to_component_indices
+  (const unsigned int i)
 {
+  Assert (rank == 2, ExcNotImplemented());
   Assert (i < n_independent_components, ExcIndexRange(i, 0, n_independent_components));
-
-  static const TableIndices<2> table[n_independent_components] =
-    { TableIndices<2> (0,0),
-      TableIndices<2> (1,1),
-      TableIndices<2> (0,1) };
-
-  return table[i];
+  switch (dim)
+    {
+    case 1:
+      return TableIndices<2>(0,0);
+    case 2:
+      {
+       static const TableIndices<2> table[3] =
+         { TableIndices<2> (0,0),
+           TableIndices<2> (1,1),
+           TableIndices<2> (0,1) };
+       return table[i];
+      }
+    case 3:
+      {
+       static const TableIndices<2> table[6] =
+         { TableIndices<2> (0,0),
+           TableIndices<2> (1,1),
+           TableIndices<2> (2,2),
+           TableIndices<2> (0,1),
+           TableIndices<2> (0,2),
+           TableIndices<2> (1,2) };
+       return table[i];
+      }
+    default:
+      Assert (false, ExcNotImplemented());
+      return TableIndices<2>(0,0);
+    }
 }
 
 
 
-template <>
-inline
-TableIndices<2>
-SymmetricTensor<2,3>::unrolled_to_component_indices (const unsigned int i)
-{
-  Assert (i < n_independent_components, ExcIndexRange(i, 0, n_independent_components));
-
-  static const TableIndices<2> table[n_independent_components] =
-    { TableIndices<2> (0,0),
-      TableIndices<2> (1,1),
-      TableIndices<2> (2,2),
-      TableIndices<2> (0,1),
-      TableIndices<2> (0,2),
-      TableIndices<2> (1,2) };
-
-  return table[i];
-}
-
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 template <class Archive>
 inline
 void
-SymmetricTensor<rank,dim>::serialize(Archive & ar, const unsigned int)
+SymmetricTensor<rank,dim,Number>::serialize(Archive & ar, const unsigned int)
 {
   ar & data;
 }
@@ -2159,47 +2079,8 @@ SymmetricTensor<rank,dim>::serialize(Archive & ar, const unsigned int)
  * Compute the determinant of a tensor or rank 2. The determinant is
  * also commonly referred to as the third invariant of rank-2 tensors.
  *
- * For the present case of one-dimensional tensors, the determinant
- * equals the only element and is therefore equivalent to the trace.
- *
- * For greater notational simplicity, there is also a <tt>third_invariant</tt>
- * function that returns the determinant of a tensor.
- *
- * @relates SymmetricTensor
- * @author Wolfgang Bangerth, 2005
- */
-template <>
-inline
-double determinant (const SymmetricTensor<2,1> &t)
-{
-  return t.data[0];
-}
-
-
-
-/**
- * Compute the determinant of a tensor or rank 2. The determinant is
- * also commonly referred to as the third invariant of rank-2 tensors.
- *
- * For greater notational simplicity, there is also a <tt>third_invariant</tt>
- * function that returns the determinant of a tensor.
- *
- * @relates SymmetricTensor
- * @author Wolfgang Bangerth, 2005
- */
-template <>
-inline
-double determinant (const SymmetricTensor<2,2> &t)
-{
-  return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
-}
-
-
-
-
-/**
- * Compute the determinant of a tensor or rank 2. The determinant is
- * also commonly referred to as the third invariant of rank-2 tensors.
+ * For a one-dimensional tensor, the determinant equals the only element and
+ * is therefore equivalent to the trace.
  *
  * For greater notational simplicity, there is also a <tt>third_invariant</tt>
  * function that returns the determinant of a tensor.
@@ -2207,18 +2088,29 @@ double determinant (const SymmetricTensor<2,2> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <>
+template <int dim, typename Number>
 inline
-double determinant (const SymmetricTensor<2,3> &t)
+Number determinant (const SymmetricTensor<2,dim,Number> &t)
 {
+  switch (dim)
+    {
+    case 1:
+      return t.data[0];
+    case 2:
+      return (t.data[0] * t.data[1] - t.data[2]*t.data[2]);
+    case 3:
                                   // in analogy to general tensors, but
                                   // there's something to be simplified for
                                   // the present case
-  return ( t.data[0]*t.data[1]*t.data[2]
-          -t.data[0]*t.data[5]*t.data[5]
-          -t.data[1]*t.data[4]*t.data[4]
-          -t.data[2]*t.data[3]*t.data[3]
-          +2*t.data[3]*t.data[4]*t.data[5] );
+      return ( t.data[0]*t.data[1]*t.data[2]
+              -t.data[0]*t.data[5]*t.data[5]
+              -t.data[1]*t.data[4]*t.data[4]
+              -t.data[2]*t.data[3]*t.data[3]
+              +2*t.data[3]*t.data[4]*t.data[5] );
+    default:
+      Assert (false, ExcNotImplemented());
+      return 0;
+    }
 }
 
 
@@ -2232,9 +2124,9 @@ double determinant (const SymmetricTensor<2,3> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double third_invariant (const SymmetricTensor<2,dim> &t)
+double third_invariant (const SymmetricTensor<2,dim,Number> &t)
 {
   return determinant (t);
 }
@@ -2249,10 +2141,10 @@ double third_invariant (const SymmetricTensor<2,dim> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
-double trace (const SymmetricTensor<2,dim> &d)
+template <int dim, typename Number>
+Number trace (const SymmetricTensor<2,dim,Number> &d)
 {
-  double t=0;
+  Number t=0;
   for (unsigned int i=0; i<dim; ++i)
     t += d.data[i];
   return t;
@@ -2268,9 +2160,9 @@ double trace (const SymmetricTensor<2,dim> &d)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double first_invariant (const SymmetricTensor<2,dim> &t)
+Number first_invariant (const SymmetricTensor<2,dim,Number> &t)
 {
   return trace (t);
 }
@@ -2283,8 +2175,9 @@ double first_invariant (const SymmetricTensor<2,dim> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005, 2010
  */
+template <typename Number>
 inline
-double second_invariant (const SymmetricTensor<2,1> &)
+Number second_invariant (const SymmetricTensor<2,1,Number> &)
 {
   return 0;
 }
@@ -2298,8 +2191,9 @@ double second_invariant (const SymmetricTensor<2,1> &)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005, 2010
  */
+template <typename Number>
 inline
-double second_invariant (const SymmetricTensor<2,2> &t)
+Number second_invariant (const SymmetricTensor<2,2,Number> &t)
 {
   return t[0][0]*t[1][1] - t[0][1]*t[0][1];
 }
@@ -2313,8 +2207,9 @@ double second_invariant (const SymmetricTensor<2,2> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005, 2010
  */
+template <typename Number>
 inline
-double second_invariant (const SymmetricTensor<2,3> &t)
+Number second_invariant (const SymmetricTensor<2,3,Number> &t)
 {
   return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0]
          - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]);
@@ -2332,10 +2227,10 @@ double second_invariant (const SymmetricTensor<2,3> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
 inline
-SymmetricTensor<rank,dim>
-transpose (const SymmetricTensor<rank,dim> &t)
+SymmetricTensor<rank,dim,Number>
+transpose (const SymmetricTensor<rank,dim,Number> &t)
 {
   return t;
 }
@@ -2351,15 +2246,15 @@ transpose (const SymmetricTensor<rank,dim> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-SymmetricTensor<2,dim>
-deviator (const SymmetricTensor<2,dim> &t)
+SymmetricTensor<2,dim,Number>
+deviator (const SymmetricTensor<2,dim,Number> &t)
 {
-  SymmetricTensor<2,dim> tmp = t;
+  SymmetricTensor<2,dim,Number> tmp = t;
 
                                    // subtract scaled trace from the diagonal
-  const double tr = trace(t) / dim;
+  const Number tr = trace(t) / dim;
   for (unsigned int i=0; i<dim; ++i)
     tmp.data[i] -= tr;
 
@@ -2369,54 +2264,32 @@ deviator (const SymmetricTensor<2,dim> &t)
 
 
 /**
- * Return a unit symmetric tensor of rank 2 and dimension 1.
- *
- * @relates SymmetricTensor
- * @author Wolfgang Bangerth, 2005
- */
-template <>
-inline
-SymmetricTensor<2,1>
-unit_symmetric_tensor<1> ()
-{
-  SymmetricTensor<2,1> tmp;
-  tmp.data[0] = 1;
-  return tmp;
-}
-
-
-
-/**
- * Return a unit symmetric tensor of rank 2 and dimension 2.
- *
- * @relates SymmetricTensor
- * @author Wolfgang Bangerth, 2005
- */
-template <>
-inline
-SymmetricTensor<2,2>
-unit_symmetric_tensor<2> ()
-{
-  SymmetricTensor<2,2> tmp;
-  tmp.data[0] = tmp.data[1] = 1;
-  return tmp;
-}
-
-
-
-/**
- * Return a unit symmetric tensor of rank 2 and dimension 3.
+ * Return a unit symmetric tensor of rank 2.
  *
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <>
+template <int dim, typename Number>
 inline
-SymmetricTensor<2,3>
-unit_symmetric_tensor<3> ()
+SymmetricTensor<2,dim,Number>
+unit_symmetric_tensor ()
 {
-  SymmetricTensor<2,3> tmp;
-  tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
+  SymmetricTensor<2,dim,Number> tmp;
+  switch (dim)
+    {
+    case 1:
+      tmp.data[0] = 1;
+      break;
+    case 2:
+      tmp.data[0] = tmp.data[1] = 1;
+      break;
+    case 3:
+      tmp.data[0] = tmp.data[1] = tmp.data[2] = 1;
+      break;
+    default:
+      for (unsigned int d=0; d<dim; ++d)
+       tmp.data[d] = 1;
+    }
   return tmp;
 }
 
@@ -2437,12 +2310,12 @@ unit_symmetric_tensor<3> ()
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-SymmetricTensor<4,dim>
+SymmetricTensor<4,dim,Number>
 deviator_tensor ()
 {
-  SymmetricTensor<4,dim> tmp;
+  SymmetricTensor<4,dim,Number> tmp;
 
                                    // fill the elements treating the diagonal
   for (unsigned int i=0; i<dim; ++i)
@@ -2455,7 +2328,7 @@ deviator_tensor ()
                                    // off-diagonal elements twice, so simply
                                    // copying requires a weight of 1/2
   for (unsigned int i=dim;
-       i<internal::SymmetricTensorAccessors::StorageType<4,dim>::n_rank2_components;
+       i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
        ++i)
     tmp.data[i][i] = 0.5;
 
@@ -2485,12 +2358,12 @@ deviator_tensor ()
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-SymmetricTensor<4,dim>
+SymmetricTensor<4,dim,Number>
 identity_tensor ()
 {
-  SymmetricTensor<4,dim> tmp;
+  SymmetricTensor<4,dim,Number> tmp;
 
                                    // fill the elements treating the diagonal
   for (unsigned int i=0; i<dim; ++i)
@@ -2502,7 +2375,7 @@ identity_tensor ()
                                    // off-diagonal elements twice, so simply
                                    // copying requires a weight of 1/2
   for (unsigned int i=dim;
-       i<internal::SymmetricTensorAccessors::StorageType<4,dim>::n_rank2_components;
+       i<internal::SymmetricTensorAccessors::StorageType<4,dim,Number>::n_rank2_components;
        ++i)
     tmp.data[i][i] = 0.5;
 
@@ -2515,7 +2388,7 @@ identity_tensor ()
  * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are
  * mappings from and to symmetric rank-2 tensors, they can have an
  * inverse. This function computes it, if it exists, for the case that the
- * dimension equals 1.
+ * dimension equals either 1 or 2.
  *
  * If a tensor is not invertible, then the result is unspecified, but will
  * likely contain the results of a division by zero or a very small number at
@@ -2524,37 +2397,18 @@ identity_tensor ()
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <>
+template <int dim, typename Number>
 inline
-SymmetricTensor<4,1>
-invert (const SymmetricTensor<4,1> &t)
+SymmetricTensor<4,dim,Number>
+invert (const SymmetricTensor<4,dim,Number> &t)
 {
-  SymmetricTensor<4,1> tmp;
-  tmp.data[0][0] = 1./t.data[0][0];
-  return tmp;
-}
-
-
-
-/**
- * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are
- * mappings from and to symmetric rank-2 tensors, they can have an
- * inverse. This function computes it, if it exists, for the case that the
- * dimension equals 2.
- *
- * If a tensor is not invertible, then the result is unspecified, but will
- * likely contain the results of a division by zero or a very small number at
- * the very least.
- *
- * @relates SymmetricTensor
- * @author Wolfgang Bangerth, 2005
- */
-template <>
-inline
-SymmetricTensor<4,2>
-invert (const SymmetricTensor<4,2> &t)
-{
-  SymmetricTensor<4,2> tmp;
+  SymmetricTensor<4,dim,Number> tmp;
+  switch (dim)
+    {
+    case 1:
+      tmp.data[0][0] = 1./t.data[0][0];
+      break;
+    case 2:
 
                                    // inverting this tensor is a little more
                                    // complicated than necessary, since we
@@ -2588,32 +2442,38 @@ invert (const SymmetricTensor<4,2> &t)
                                    // column of the inverse by 1/2,
                                    // corresponding to the left and right
                                    // multiplication with mult^-1
-  const double t4 = t.data[0][0]*t.data[1][1],
-               t6 = t.data[0][0]*t.data[1][2],
-               t8 = t.data[0][1]*t.data[1][0],
-              t00 = t.data[0][2]*t.data[1][0],
-              t01 = t.data[0][1]*t.data[2][0],
-              t04 = t.data[0][2]*t.data[2][0],
-              t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-t8*t.data[2][2]+
-                         t00*t.data[2][1]+t01*t.data[1][2]-t04*t.data[1][1]);
-  tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
-  tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
-  tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
-  tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
-  tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
-  tmp.data[1][2] = -(t6-t00)*t07;
-  tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
-  tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
-  tmp.data[2][2] = (t4-t8)*t07;
+      {
+       const Number t4 = t.data[0][0]*t.data[1][1],
+                    t6 = t.data[0][0]*t.data[1][2],
+                    t8 = t.data[0][1]*t.data[1][0],
+                     t00 = t.data[0][2]*t.data[1][0],
+                     t01 = t.data[0][1]*t.data[2][0],
+                     t04 = t.data[0][2]*t.data[2][0],
+                     t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-
+                               t8*t.data[2][2]+t00*t.data[2][1]+
+                               t01*t.data[1][2]-t04*t.data[1][1]);
+       tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07;
+       tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07;
+       tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07;
+       tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07;
+       tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07;
+       tmp.data[1][2] = -(t6-t00)*t07;
+       tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07;
+       tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07;
+       tmp.data[2][2] = (t4-t8)*t07;
 
                                    // scale last row and column as mentioned
                                    // above
-  tmp.data[2][0] /= 2;
-  tmp.data[2][1] /= 2;
-  tmp.data[0][2] /= 2;
-  tmp.data[1][2] /= 2;
-  tmp.data[2][2] /= 4;
-
+       tmp.data[2][0] /= 2;
+       tmp.data[2][1] /= 2;
+       tmp.data[0][2] /= 2;
+       tmp.data[1][2] /= 2;
+       tmp.data[2][2] /= 4;
+      }
+      break;
+    default:
+      Assert (false, ExcNotImplemented());
+    }
   return tmp;
 }
 
@@ -2633,9 +2493,9 @@ invert (const SymmetricTensor<4,2> &t)
  * @author Wolfgang Bangerth, 2005
  */
 template <>
-SymmetricTensor<4,3>
-invert (const SymmetricTensor<4,3> &t);
-// this function is implemented in the .cc file
+SymmetricTensor<4,3,double>
+invert (const SymmetricTensor<4,3,double> &t);
+// this function is implemented in the .cc file for double data types
 
 
 /**
@@ -2644,21 +2504,21 @@ invert (const SymmetricTensor<4,3> &t);
  * <tt>T phi = t1 (t2 : phi)</tt> for all symmetric tensors <tt>phi</tt>.
  *
  * For example, the deviator tensor can be computed as
- * <tt>identity_tensor<dim>() -
- * 1/d*outer_product(unit_symmetric_tensor<dim>(),
- * unit_symmetric_tensor<dim>())</tt>, since the (double) contraction with the
+ * <tt>identity_tensor<dim,Number>() -
+ * 1/d*outer_product(unit_symmetric_tensor<dim,Number>(),
+ * unit_symmetric_tensor<dim,Number>())</tt>, since the (double) contraction with the
  * unit tensor yields the trace of a symmetric tensor.
  *
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-SymmetricTensor<4,dim>
-outer_product (const SymmetricTensor<2,dim> &t1,
-               const SymmetricTensor<2,dim> &t2)
+SymmetricTensor<4,dim,Number>
+outer_product (const SymmetricTensor<2,dim,Number> &t1,
+               const SymmetricTensor<2,dim,Number> &t2)
 {
-  SymmetricTensor<4,dim> tmp;
+  SymmetricTensor<4,dim,Number> tmp;
 
                                    // fill only the elements really needed
   for (unsigned int i=0; i<dim; ++i)
@@ -2680,13 +2540,14 @@ outer_product (const SymmetricTensor<2,dim> &t1,
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
-SymmetricTensor<2,1>
-symmetrize (const Tensor<2,1> &t)
+SymmetricTensor<2,1,Number>
+symmetrize (const Tensor<2,1,Number> &t)
 {
-  const double array[1]
+  const Number array[1]
     = { t[0][0] };
-  return SymmetricTensor<2,1>(array);
+  return SymmetricTensor<2,1,Number>(array);
 }
 
 
@@ -2699,13 +2560,14 @@ symmetrize (const Tensor<2,1> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
-SymmetricTensor<2,2>
-symmetrize (const Tensor<2,2> &t)
+SymmetricTensor<2,2,Number>
+symmetrize (const Tensor<2,2,Number> &t)
 {
-  const double array[3]
+  const Number array[3]
     = { t[0][0], t[1][1], (t[0][1] + t[1][0])/2 };
-  return SymmetricTensor<2,2>(array);
+  return SymmetricTensor<2,2,Number>(array);
 }
 
 
@@ -2718,21 +2580,79 @@ symmetrize (const Tensor<2,2> &t)
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
-SymmetricTensor<2,3>
-symmetrize (const Tensor<2,3> &t)
+SymmetricTensor<2,3,Number>
+symmetrize (const Tensor<2,3,Number> &t)
 {
-  const double array[6]
+  const Number array[6]
     = { t[0][0], t[1][1], t[2][2],
        (t[0][1] + t[1][0])/2,
        (t[0][2] + t[2][0])/2,
        (t[1][2] + t[2][1])/2 };
-  return SymmetricTensor<2,3>(array);
+  return SymmetricTensor<2,3,Number>(array);
 }
 
 
+
+/**
+ * Multiplication of a symmetric tensor of general rank with a scalar
+ * from the right.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim, typename Number>
+inline
+SymmetricTensor<rank,dim,Number>
+operator * (const SymmetricTensor<rank,dim,Number> &t,
+           const Number                            factor)
+{
+  SymmetricTensor<rank,dim,Number> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Multiplication of a symmetric tensor of general rank with a scalar
+ * from the left.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim, typename Number>
+inline
+SymmetricTensor<rank,dim,Number>
+operator * (const Number                            factor,
+           const SymmetricTensor<rank,dim,Number> &t)
+{
+  SymmetricTensor<rank,dim,Number> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Division of a symmetric tensor of general rank by a scalar.
+ *
+ * @relates SymmetricTensor
+ */
+template <int rank, int dim, typename Number>
+inline
+SymmetricTensor<rank,dim,Number>
+operator / (const SymmetricTensor<rank,dim,Number> &t,
+           const Number                            factor)
+{
+  SymmetricTensor<rank,dim,Number> tt = t;
+  tt /= factor;
+  return tt;
+}
+
+
+
 /**
- * Multiplication of a symmetric tensor of general rank with a scalar double
+ * Multiplication of a symmetric tensor of general rank with a scalar
  * from the right.
  *
  * @relates SymmetricTensor
@@ -2741,7 +2661,7 @@ template <int rank, int dim>
 inline
 SymmetricTensor<rank,dim>
 operator * (const SymmetricTensor<rank,dim> &t,
-           const double            factor)
+           const double                     factor)
 {
   SymmetricTensor<rank,dim> tt = t;
   tt *= factor;
@@ -2751,7 +2671,7 @@ operator * (const SymmetricTensor<rank,dim> &t,
 
 
 /**
- * Multiplication of a symmetric tensor of general rank with a scalar double
+ * Multiplication of a symmetric tensor of general rank with a scalar
  * from the left.
  *
  * @relates SymmetricTensor
@@ -2759,7 +2679,7 @@ operator * (const SymmetricTensor<rank,dim> &t,
 template <int rank, int dim>
 inline
 SymmetricTensor<rank,dim>
-operator * (const double            factor,
+operator * (const double                     factor,
            const SymmetricTensor<rank,dim> &t)
 {
   SymmetricTensor<rank,dim> tt = t;
@@ -2770,7 +2690,7 @@ operator * (const double            factor,
 
 
 /**
- * Division of a symmetric tensor of general rank by a scalar double.
+ * Division of a symmetric tensor of general rank by a scalar.
  *
  * @relates SymmetricTensor
  */
@@ -2778,7 +2698,7 @@ template <int rank, int dim>
 inline
 SymmetricTensor<rank,dim>
 operator / (const SymmetricTensor<rank,dim> &t,
-           const double            factor)
+           const double                     factor)
 {
   SymmetricTensor<rank,dim> tt = t;
   tt /= factor;
@@ -2803,11 +2723,12 @@ operator / (const SymmetricTensor<rank,dim> &t,
  * @related SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,1> &tmp,
-                const SymmetricTensor<4,1> &t,
-                const SymmetricTensor<2,1> &s)
+double_contract (SymmetricTensor<2,1,Number> &tmp,
+                const SymmetricTensor<4,1,Number> &t,
+                const SymmetricTensor<2,1,Number> &s)
 {
   tmp[0][0] = t[0][0][0][0] * s[0][0];
 }
@@ -2829,11 +2750,12 @@ double_contract (SymmetricTensor<2,1> &tmp,
  * @related SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,1> &tmp,
-                const SymmetricTensor<2,1> &s,
-                const SymmetricTensor<4,1> &t)
+double_contract (SymmetricTensor<2,1,Number> &tmp,
+                const SymmetricTensor<2,1,Number> &s,
+                const SymmetricTensor<4,1,Number> &t)
 {
   tmp[0][0] = t[0][0][0][0] * s[0][0];
 }
@@ -2854,11 +2776,12 @@ double_contract (SymmetricTensor<2,1> &tmp,
  *
  * @related SymmetricTensor @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,2> &tmp,
-                const SymmetricTensor<4,2> &t,
-                const SymmetricTensor<2,2> &s)
+double_contract (SymmetricTensor<2,2,Number> &tmp,
+                const SymmetricTensor<4,2,Number> &t,
+                const SymmetricTensor<2,2,Number> &s)
 {
   const unsigned int dim = 2;
 
@@ -2886,11 +2809,12 @@ double_contract (SymmetricTensor<2,2> &tmp,
  * @related SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,2> &tmp,
-                const SymmetricTensor<2,2> &s,
-                const SymmetricTensor<4,2> &t)
+double_contract (SymmetricTensor<2,2,Number> &tmp,
+                const SymmetricTensor<2,2,Number> &s,
+                const SymmetricTensor<4,2,Number> &t)
 {
   const unsigned int dim = 2;
 
@@ -2918,11 +2842,12 @@ double_contract (SymmetricTensor<2,2> &tmp,
  * @related SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,3> &tmp,
-                const SymmetricTensor<4,3> &t,
-                const SymmetricTensor<2,3> &s)
+double_contract (SymmetricTensor<2,3,Number> &tmp,
+                const SymmetricTensor<4,3,Number> &t,
+                const SymmetricTensor<2,3,Number> &s)
 {
   const unsigned int dim = 3;
 
@@ -2953,11 +2878,12 @@ double_contract (SymmetricTensor<2,3> &tmp,
  * @related SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
+template <typename Number>
 inline
 void
-double_contract (SymmetricTensor<2,3> &tmp,
-                const SymmetricTensor<2,3> &s,
-                const SymmetricTensor<4,3> &t)
+double_contract (SymmetricTensor<2,3,Number> &tmp,
+                const SymmetricTensor<2,3,Number> &s,
+                const SymmetricTensor<4,3,Number> &t)
 {
   const unsigned int dim = 3;
 
@@ -2989,12 +2915,12 @@ double_contract (SymmetricTensor<2,3> &tmp,
  * @relates SymmetricTensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
-Tensor<1,dim>
-operator * (const SymmetricTensor<2,dim> &src1,
-            const Tensor<1,dim> &src2)
+template <int dim, typename Number>
+Tensor<1,dim,Number>
+operator * (const SymmetricTensor<2,dim,Number> &src1,
+            const Tensor<1,dim,Number> &src2)
 {
-  Tensor<1,dim> dest;
+  Tensor<1,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       dest[i] += src1[i][j] * src2[j];
@@ -3011,15 +2937,15 @@ operator * (const SymmetricTensor<2,dim> &src1,
  *
  * @relates SymmetricTensor
  */
-template <int dim>
+template <int dim, typename Number>
 inline
 std::ostream & operator << (std::ostream &out,
-                           const SymmetricTensor<2,dim> &t)
+                           const SymmetricTensor<2,dim,Number> &t)
 {
                                   //make out lives a bit simpler by outputing
                                   //the tensor through the operator for the
                                   //general Tensor class
-  Tensor<2,dim> tt;
+  Tensor<2,dim,Number> tt;
 
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
@@ -3039,15 +2965,15 @@ std::ostream & operator << (std::ostream &out,
  *
  * @relates SymmetricTensor
  */
-template <int dim>
+template <int dim, typename Number>
 inline
 std::ostream & operator << (std::ostream &out,
-                           const SymmetricTensor<4,dim> &t)
+                           const SymmetricTensor<4,dim,Number> &t)
 {
                                   //make out lives a bit simpler by outputing
                                   //the tensor through the operator for the
                                   //general Tensor class
-  Tensor<4,dim> tt;
+  Tensor<4,dim,Number> tt;
 
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
index 0713fdb5d5c883d17283db5481b9eff44b70c191..8cfb1f71e6e7ded1aa97c01499d1facbf9c9141f 100644 (file)
@@ -17,8 +17,8 @@
 #include <deal.II/base/tensor_base.h>
 
 DEAL_II_NAMESPACE_OPEN
-template <int rank_, int dim> class Tensor;
-template <int dim> class Tensor<1,dim>;
+template <int rank_, int dim, typename Number> class Tensor;
+template <int dim, typename Number> class Tensor<1,dim,Number>;
 
 /**
  * Provide a general tensor class with an arbitrary rank, i.e. with
@@ -33,10 +33,15 @@ template <int dim> class Tensor<1,dim>;
  * produce far more efficient code than for matrices with
  * runtime-dependent dimension.
  *
+ * This class provides an optional template argument for the type of the
+ * underlying data. It defaults to @p double values. It can be used to base
+ * tensors on @p float or @p complex numbers or any other data type that
+ * implements basic arithmetic operations.
+ *
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 1998-2005
  */
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 class Tensor
 {
   public:
@@ -69,7 +74,7 @@ class Tensor
                                      * Type of stored objects. This
                                      * is a tensor of lower rank.
                                      */
-    typedef Tensor<rank_-1,dim> value_type;
+    typedef Tensor<rank_-1,dim,Number> value_type;
 
                                     /**
                                      * Declare an array type which
@@ -77,7 +82,7 @@ class Tensor
                                      * object of this type
                                      * statically.
                                      */
-    typedef typename Tensor<rank_-1,dim>::array_type array_type[dim];
+    typedef typename Tensor<rank_-1,dim,Number>::array_type array_type[dim];
 
                                     /**
                                      * Constructor. Initialize all entries
@@ -94,17 +99,17 @@ class Tensor
                                     /**
                                      * Read-Write access operator.
                                      */
-    Tensor<rank_-1,dim> &operator [] (const unsigned int i);
+    Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i);
 
                                     /**
                                      * Read-only access operator.
                                      */
-    const Tensor<rank_-1,dim> &operator [] (const unsigned int i) const;
+    const Tensor<rank_-1,dim,Number> &operator [] (const unsigned int i) const;
 
                                     /**
                                      *  Assignment operator.
                                      */
-    Tensor & operator = (const Tensor<rank_,dim> &);
+    Tensor & operator = (const Tensor<rank_,dim,Number> &);
 
                                     /**
                                      * This operator assigns a scalar
@@ -118,40 +123,40 @@ class Tensor
                                      * all elements of the tensor to
                                      * zero.
                                      */
-    Tensor<rank_,dim> & operator = (const double d);
+    Tensor<rank_,dim,Number> & operator = (const Number d);
 
                                     /**
                                      *  Test for equality of two tensors.
                                      */
-    bool operator == (const Tensor<rank_,dim> &) const;
+    bool operator == (const Tensor<rank_,dim,Number> &) const;
 
                                     /**
                                      *  Test for inequality of two tensors.
                                      */
-    bool operator != (const Tensor<rank_,dim> &) const;
+    bool operator != (const Tensor<rank_,dim,Number> &) const;
 
                                     /**
                                      *  Add another tensor.
                                      */
-    Tensor<rank_,dim> & operator += (const Tensor<rank_,dim> &);
+    Tensor<rank_,dim,Number> & operator += (const Tensor<rank_,dim,Number> &);
 
                                     /**
                                      *  Subtract another tensor.
                                      */
-    Tensor<rank_,dim> & operator -= (const Tensor<rank_,dim> &);
+    Tensor<rank_,dim,Number> & operator -= (const Tensor<rank_,dim,Number> &);
 
                                     /**
                                      *  Scale the tensor by <tt>factor</tt>,
                                      *  i.e. multiply all components by
                                      *  <tt>factor</tt>.
                                      */
-    Tensor<rank_,dim> & operator *= (const double factor);
+    Tensor<rank_,dim,Number> & operator *= (const Number factor);
 
                                     /**
                                      *  Scale the vector by
                                      *  <tt>1/factor</tt>.
                                      */
-    Tensor<rank_,dim> & operator /= (const double factor);
+    Tensor<rank_,dim,Number> & operator /= (const Number factor);
 
                                     /**
                                      *  Add two tensors. If possible, you
@@ -159,7 +164,7 @@ class Tensor
                                      *  instead since this does not need the
                                      *  creation of a temporary.
                                      */
-    Tensor<rank_,dim>   operator + (const Tensor<rank_,dim> &) const;
+    Tensor<rank_,dim,Number>   operator + (const Tensor<rank_,dim,Number> &) const;
 
                                     /**
                                      *  Subtract two tensors. If possible,
@@ -167,20 +172,20 @@ class Tensor
                                      *  instead since this does not need the
                                      *  creation of a temporary.
                                      */
-    Tensor<rank_,dim>   operator - (const Tensor<rank_,dim> &) const;
+    Tensor<rank_,dim,Number>   operator - (const Tensor<rank_,dim,Number> &) const;
 
                                     /**
                                      * Unary minus operator. Negate all
                                      * entries of a tensor.
                                      */
-    Tensor<rank_,dim>   operator - () const;
+    Tensor<rank_,dim,Number>   operator - () const;
 
                                      /**
                                       * Return the Frobenius-norm of a tensor,
                                       * i.e. the square root of the sum of
                                       * squares of all entries.
                                       */
-    double norm () const;
+    Number norm () const;
 
                                      /**
                                       * Return the square of the
@@ -194,7 +199,7 @@ class Tensor
                                      * may also be useful in other
                                      * contexts.
                                       */
-    double norm_square () const;
+    Number norm_square () const;
 
                                     /**
                                      * Fill a vector with all tensor elements.
@@ -205,7 +210,7 @@ class Tensor
                                      * usual in C++, the rightmost
                                      * index of the tensor marches fastest.
                                      */
-    void unroll (Vector<double> & result) const;
+    void unroll (Vector<Number> & result) const;
 
 
                                     /**
@@ -254,12 +259,12 @@ class Tensor
                                      * Array of tensors holding the
                                      * subelements.
                                      */
-    Tensor<rank_-1,dim> subtensor[dim];
+    Tensor<rank_-1,dim,Number> subtensor[dim];
 
                                     /**
                                      * Help function for unroll.
                                      */
-    void unroll_recursion(Vector<double> &result,
+    void unroll_recursion(Vector<Number> &result,
                          unsigned int   &start_index) const;
 
                                     // make the following class a
@@ -272,7 +277,7 @@ class Tensor
                                     // also, it would be sufficient to make
                                     // the function unroll_loops a friend,
                                     // but that seems to be impossible as well.
-    template <int, int> friend class Tensor;
+    template <int, int, typename> friend class Tensor;
 };
 
 
@@ -280,9 +285,9 @@ class Tensor
 
 #ifndef DOXYGEN
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim>::Tensor ()
+Tensor<rank_,dim,Number>::Tensor ()
 {
 // default constructor. not specifying an initializer list calls
 // the default constructor of the subobjects, which initialize them
@@ -290,19 +295,19 @@ Tensor<rank_,dim>::Tensor ()
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim>::Tensor (const array_type &initializer)
+Tensor<rank_,dim,Number>::Tensor (const array_type &initializer)
 {
   for (unsigned int i=0; i<dim; ++i)
-    subtensor[i] =  Tensor<rank_-1,dim>(initializer[i]);
+    subtensor[i] = Tensor<rank_-1,dim,Number>(initializer[i]);
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-typename Tensor<rank_,dim>::value_type&
-Tensor<rank_,dim>::operator[] (const unsigned int i)
+typename Tensor<rank_,dim,Number>::value_type&
+Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
 {
   Assert (i<dim, ExcIndexRange(i, 0, dim));
 
@@ -310,10 +315,10 @@ Tensor<rank_,dim>::operator[] (const unsigned int i)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-const typename Tensor<rank_,dim>::value_type&
-Tensor<rank_,dim>::operator[] (const unsigned int i) const
+const typename Tensor<rank_,dim,Number>::value_type&
+Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
 {
   Assert (i<dim, ExcIndexRange(i, 0, dim));
 
@@ -321,10 +326,10 @@ Tensor<rank_,dim>::operator[] (const unsigned int i) const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator = (const Tensor<rank_,dim> &t)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator = (const Tensor<rank_,dim,Number> &t)
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i] = t.subtensor[i];
@@ -332,10 +337,10 @@ Tensor<rank_,dim>::operator = (const Tensor<rank_,dim> &t)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator = (const double d)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator = (const Number d)
 {
   Assert (d==0, ExcMessage ("Only assignment with zero is allowed"));
 
@@ -345,10 +350,10 @@ Tensor<rank_,dim>::operator = (const double d)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
 bool
-Tensor<rank_,dim>::operator == (const Tensor<rank_,dim> &p) const
+Tensor<rank_,dim,Number>::operator == (const Tensor<rank_,dim,Number> &p) const
 {
   for (unsigned int i=0; i<dim; ++i)
     if (subtensor[i] != p.subtensor[i]) return false;
@@ -356,19 +361,19 @@ Tensor<rank_,dim>::operator == (const Tensor<rank_,dim> &p) const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
 bool
-Tensor<rank_,dim>::operator != (const Tensor<rank_,dim> &p) const
+Tensor<rank_,dim,Number>::operator != (const Tensor<rank_,dim,Number> &p) const
 {
   return !((*this) == p);
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator += (const Tensor<rank_,dim> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator += (const Tensor<rank_,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i] += p.subtensor[i];
@@ -376,10 +381,10 @@ Tensor<rank_,dim>::operator += (const Tensor<rank_,dim> &p)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator -= (const Tensor<rank_,dim> &p)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator -= (const Tensor<rank_,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i] -= p.subtensor[i];
@@ -387,10 +392,10 @@ Tensor<rank_,dim>::operator -= (const Tensor<rank_,dim> &p)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator *= (const double s)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator *= (const Number s)
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i] *= s;
@@ -398,10 +403,10 @@ Tensor<rank_,dim>::operator *= (const double s)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim> &
-Tensor<rank_,dim>::operator /= (const double s)
+Tensor<rank_,dim,Number> &
+Tensor<rank_,dim,Number>::operator /= (const Number s)
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i] /= s;
@@ -409,12 +414,12 @@ Tensor<rank_,dim>::operator /= (const double s)
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim>
-Tensor<rank_,dim>::operator + (const Tensor<rank_,dim> &t) const
+Tensor<rank_,dim,Number>
+Tensor<rank_,dim,Number>::operator + (const Tensor<rank_,dim,Number> &t) const
 {
-  Tensor<rank_,dim> tmp(*this);
+  Tensor<rank_,dim,Number> tmp(*this);
 
   for (unsigned int i=0; i<dim; ++i)
     tmp.subtensor[i] += t.subtensor[i];
@@ -423,12 +428,12 @@ Tensor<rank_,dim>::operator + (const Tensor<rank_,dim> &t) const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim>
-Tensor<rank_,dim>::operator - (const Tensor<rank_,dim> &t) const
+Tensor<rank_,dim,Number>
+Tensor<rank_,dim,Number>::operator - (const Tensor<rank_,dim,Number> &t) const
 {
-  Tensor<rank_,dim> tmp(*this);
+  Tensor<rank_,dim,Number> tmp(*this);
 
   for (unsigned int i=0; i<dim; ++i)
     tmp.subtensor[i] -= t.subtensor[i];
@@ -437,12 +442,12 @@ Tensor<rank_,dim>::operator - (const Tensor<rank_,dim> &t) const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim>
-Tensor<rank_,dim>::operator - () const
+Tensor<rank_,dim,Number>
+Tensor<rank_,dim,Number>::operator - () const
 {
-  Tensor<rank_,dim> tmp;
+  Tensor<rank_,dim,Number> tmp;
 
   for (unsigned int i=0; i<dim; ++i)
     tmp.subtensor[i] = -subtensor[i];
@@ -451,19 +456,19 @@ Tensor<rank_,dim>::operator - () const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-double Tensor<rank_,dim>::norm () const
+Number Tensor<rank_,dim,Number>::norm () const
 {
   return std::sqrt (norm_square());
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-double Tensor<rank_,dim>::norm_square () const
+Number Tensor<rank_,dim,Number>::norm_square () const
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     s += subtensor[i].norm_square();
 
@@ -471,29 +476,29 @@ double Tensor<rank_,dim>::norm_square () const
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-void Tensor<rank_,dim>::clear ()
+void Tensor<rank_,dim,Number>::clear ()
 {
   for (unsigned int i=0; i<dim; ++i)
     subtensor[i].clear();
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
 std::size_t
-Tensor<rank_,dim>::memory_consumption ()
+Tensor<rank_,dim,Number>::memory_consumption ()
 {
-  return sizeof(Tensor<rank_,dim>);
+  return sizeof(Tensor<rank_,dim,Number>);
 }
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 template <class Archive>
 inline
 void
-Tensor<rank_,dim>::serialize(Archive & ar, const unsigned int)
+Tensor<rank_,dim,Number>::serialize(Archive & ar, const unsigned int)
 {
   ar & subtensor;
 }
@@ -510,9 +515,9 @@ Tensor<rank_,dim>::serialize(Archive & ar, const unsigned int)
  *
  * @relates Tensor
  */
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 inline
-std::ostream & operator << (std::ostream &out, const Tensor<rank_,dim> &p)
+std::ostream & operator << (std::ostream &out, const Tensor<rank_,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim-1; ++i)
     out << p[i] << ' ';
@@ -545,12 +550,12 @@ std::ostream & operator << (std::ostream &out, const Tensor<rank_,1> &p)
  * @relates Tensor
  * @author Guido Kanschat, 2000
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double contract (const Tensor<1,dim> &src1,
-                const Tensor<1,dim> &src2)
+Number contract (const Tensor<1,dim,Number> &src1,
+                const Tensor<1,dim,Number> &src2)
 {
-  double res = 0.;
+  Number res = 0.;
   for (unsigned int i=0; i<dim; ++i)
     res += src1[i] * src2[i];
 
@@ -574,11 +579,11 @@ double contract (const Tensor<1,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double
-operator * (const Tensor<1,dim> &src1,
-            const Tensor<1,dim> &src2)
+Number
+operator * (const Tensor<1,dim,Number> &src1,
+            const Tensor<1,dim,Number> &src2)
 {
   return contract(src1, src2);
 }
@@ -591,11 +596,11 @@ operator * (const Tensor<1,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<1,dim>       &dest,
-              const Tensor<2,dim> &src1,
-              const Tensor<1,dim> &src2)
+void contract (Tensor<1,dim,Number>       &dest,
+              const Tensor<2,dim,Number> &src1,
+              const Tensor<1,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -620,12 +625,12 @@ void contract (Tensor<1,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
-Tensor<1,dim>
-operator * (const Tensor<2,dim> &src1,
-            const Tensor<1,dim> &src2)
+template <int dim, typename Number>
+Tensor<1,dim,Number>
+operator * (const Tensor<2,dim,Number> &src1,
+            const Tensor<1,dim,Number> &src2)
 {
-  Tensor<1,dim> dest;
+  Tensor<1,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       dest[i] += src1[i][j] * src2[j];
@@ -640,11 +645,11 @@ operator * (const Tensor<2,dim> &src1,
  * @relates Tensor
  * @author Guido Kanschat, 2001
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<1,dim>       &dest,
-              const Tensor<1,dim> &src1,
-              const Tensor<2,dim> &src2)
+void contract (Tensor<1,dim,Number>       &dest,
+              const Tensor<1,dim,Number> &src1,
+              const Tensor<2,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -669,13 +674,13 @@ void contract (Tensor<1,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim>
-operator * (const Tensor<1,dim> &src1,
-            const Tensor<2,dim> &src2)
+Tensor<1,dim,Number>
+operator * (const Tensor<1,dim,Number> &src1,
+            const Tensor<2,dim,Number> &src2)
 {
-  Tensor<1,dim> dest;
+  Tensor<1,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       dest[i] += src1[j] * src2[j][i];
@@ -690,11 +695,11 @@ operator * (const Tensor<1,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<2,dim>       &dest,
-              const Tensor<2,dim> &src1,
-              const Tensor<2,dim> &src2)
+void contract (Tensor<2,dim,Number>       &dest,
+              const Tensor<2,dim,Number> &src1,
+              const Tensor<2,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -720,13 +725,13 @@ void contract (Tensor<2,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<2,dim>
-operator * (const Tensor<2,dim> &src1,
-            const Tensor<2,dim> &src2)
+Tensor<2,dim,Number>
+operator * (const Tensor<2,dim,Number> &src1,
+            const Tensor<2,dim,Number> &src2)
 {
-  Tensor<2,dim> dest;
+  Tensor<2,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -749,11 +754,11 @@ operator * (const Tensor<2,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<2,dim>       &dest,
-              const Tensor<2,dim> &src1,   const unsigned int index1,
-              const Tensor<2,dim> &src2,   const unsigned int index2)
+void contract (Tensor<2,dim,Number>       &dest,
+              const Tensor<2,dim,Number> &src1,   const unsigned int index1,
+              const Tensor<2,dim,Number> &src2,   const unsigned int index2)
 {
   dest.clear ();
 
@@ -777,7 +782,7 @@ void contract (Tensor<2,dim>       &dest,
 
                default:
                      Assert (false,
-                              (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+                              (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
              };
            break;
       case 2:
@@ -798,12 +803,12 @@ void contract (Tensor<2,dim>       &dest,
 
                default:
                      Assert (false,
-                              (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+                              (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
              };
            break;
 
       default:
-           Assert (false, (typename Tensor<2,dim>::ExcInvalidTensorIndex (index1)));
+           Assert (false, (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
     };
 }
 
@@ -819,11 +824,11 @@ void contract (Tensor<2,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<2,dim>       &dest,
-              const Tensor<3,dim> &src1,   const unsigned int index1,
-              const Tensor<1,dim> &src2)
+void contract (Tensor<2,dim,Number>       &dest,
+              const Tensor<3,dim,Number> &src1,   const unsigned int index1,
+              const Tensor<1,dim,Number> &src2)
 {
   dest.clear ();
 
@@ -852,7 +857,7 @@ void contract (Tensor<2,dim>       &dest,
 
       default:
            Assert (false,
-                    (typename Tensor<2,dim>::ExcInvalidTensorIndex (index1)));
+                    (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
     };
 }
 
@@ -864,11 +869,11 @@ void contract (Tensor<2,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<3,dim>       &dest,
-              const Tensor<3,dim> &src1,
-              const Tensor<2,dim> &src2)
+void contract (Tensor<3,dim,Number>       &dest,
+              const Tensor<3,dim,Number> &src1,
+              const Tensor<2,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -892,11 +897,11 @@ void contract (Tensor<3,dim>       &dest,
  *
  * @relates Tensor
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<3,dim>       &dest,
-              const Tensor<3,dim> &src1, const unsigned int index1,
-              const Tensor<2,dim> &src2, const unsigned int index2)
+void contract (Tensor<3,dim,Number>       &dest,
+              const Tensor<3,dim,Number> &src1, const unsigned int index1,
+              const Tensor<2,dim,Number> &src2, const unsigned int index2)
 {
   dest.clear ();
 
@@ -921,7 +926,7 @@ void contract (Tensor<3,dim>       &dest,
                      break;
                default:
                      Assert (false,
-                             (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+                             (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
              }
 
            break;
@@ -944,7 +949,7 @@ void contract (Tensor<3,dim>       &dest,
                      break;
                default:
                      Assert (false,
-                             (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+                             (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
              }
 
            break;
@@ -967,13 +972,13 @@ void contract (Tensor<3,dim>       &dest,
                      break;
                default:
                      Assert (false,
-                             (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+                             (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
              }
 
            break;
       default:
            Assert (false,
-                   (typename Tensor<3,dim>::ExcInvalidTensorIndex (index1)));
+                   (typename Tensor<3,dim,Number>::ExcInvalidTensorIndex (index1)));
     }
 }
 
@@ -994,13 +999,13 @@ void contract (Tensor<3,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<3,dim>
-operator * (const Tensor<3,dim> &src1,
-            const Tensor<2,dim> &src2)
+Tensor<3,dim,Number>
+operator * (const Tensor<3,dim,Number> &src1,
+            const Tensor<2,dim,Number> &src2)
 {
-  Tensor<3,dim> dest;
+  Tensor<3,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1017,11 +1022,11 @@ operator * (const Tensor<3,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void contract (Tensor<3,dim>       &dest,
-              const Tensor<2,dim> &src1,
-              const Tensor<3,dim> &src2)
+void contract (Tensor<3,dim,Number>       &dest,
+              const Tensor<2,dim,Number> &src1,
+              const Tensor<3,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -1048,13 +1053,13 @@ void contract (Tensor<3,dim>       &dest,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<3,dim>
-operator * (const Tensor<2,dim> &src1,
-            const Tensor<3,dim> &src2)
+Tensor<3,dim,Number>
+operator * (const Tensor<2,dim,Number> &src1,
+            const Tensor<3,dim,Number> &src2)
 {
-  Tensor<3,dim> dest;
+  Tensor<3,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1071,13 +1076,13 @@ operator * (const Tensor<2,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<4,dim>
-operator * (const Tensor<3,dim> &src1,
-            const Tensor<3,dim> &src2)
+Tensor<4,dim,Number>
+operator * (const Tensor<3,dim,Number> &src1,
+            const Tensor<3,dim,Number> &src2)
 {
-  Tensor<4,dim> dest;
+  Tensor<4,dim,Number> dest;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1096,11 +1101,11 @@ operator * (const Tensor<3,dim> &src1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2005
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-void double_contract (Tensor<2,dim>       &dest,
-                     const Tensor<4,dim> &src1,
-                     const Tensor<2,dim> &src2)
+void double_contract (Tensor<2,dim,Number>       &dest,
+                     const Tensor<4,dim,Number> &src1,
+                     const Tensor<2,dim,Number> &src2)
 {
   dest.clear ();
   for (unsigned int i=0; i<dim; ++i)
@@ -1118,13 +1123,13 @@ void double_contract (Tensor<2,dim>       &dest,
  * @relates Tensor
  * @author Guido Kanschat, 2004
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double contract3 (const Tensor<1,dim>& u,
-                 const Tensor<2,dim>& A,
-                 const Tensor<1,dim>& v)
+Number contract3 (const Tensor<1,dim,Number>& u,
+                 const Tensor<2,dim,Number>& A,
+                 const Tensor<1,dim,Number>& v)
 {
-  double result = 0.;
+  Number result = 0.;
 
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
@@ -1140,14 +1145,14 @@ double contract3 (const Tensor<1,dim>& u,
  * @relates Tensor
  * @author Toby D. Young, 2011
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double
-contract3 (const Tensor<1,dim> &t1,
-          const Tensor<3,dim> &t2,
-          const Tensor<2,dim> &t3)
+Number
+contract3 (const Tensor<1,dim,Number> &t1,
+          const Tensor<3,dim,Number> &t2,
+          const Tensor<2,dim,Number> &t3)
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1163,14 +1168,14 @@ contract3 (const Tensor<1,dim> &t1,
  * @relates Tensor
  * @author Toby D. Young, 2011
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double
-contract3 (const Tensor<2,dim> &t1,
-          const Tensor<3,dim> &t2,
-          const Tensor<1,dim> &t3)
+Number
+contract3 (const Tensor<2,dim,Number> &t1,
+          const Tensor<3,dim,Number> &t2,
+          const Tensor<1,dim,Number> &t3)
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1186,14 +1191,14 @@ contract3 (const Tensor<2,dim> &t1,
  * @relates Tensor
  * @author Toby D. Young, 2011
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double
-contract3 (const Tensor<2,dim> &t1,
-          const Tensor<4,dim> &t2,
-          const Tensor<2,dim> &t3)
+Number
+contract3 (const Tensor<2,dim,Number> &t1,
+          const Tensor<4,dim,Number> &t2,
+          const Tensor<2,dim,Number> &t3)
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       for (unsigned int k=0; k<dim; ++k)
@@ -1210,10 +1215,10 @@ contract3 (const Tensor<2,dim> &t1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
-void outer_product (Tensor<2,dim>       &dst,
-                   const Tensor<1,dim> &src1,
-                   const Tensor<1,dim> &src2)
+template <int dim, typename Number>
+void outer_product (Tensor<2,dim,Number>       &dst,
+                   const Tensor<1,dim,Number> &src1,
+                   const Tensor<1,dim,Number> &src2)
 {
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
@@ -1228,10 +1233,10 @@ void outer_product (Tensor<2,dim>       &dst,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
-void outer_product (Tensor<3,dim>       &dst,
-                   const Tensor<1,dim> &src1,
-                   const Tensor<2,dim> &src2)
+template <int dim, typename Number>
+void outer_product (Tensor<3,dim,Number>       &dst,
+                   const Tensor<1,dim,Number> &src1,
+                   const Tensor<2,dim,Number> &src2)
 {
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
@@ -1247,10 +1252,10 @@ void outer_product (Tensor<3,dim>       &dst,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
-void outer_product (Tensor<3,dim>       &dst,
-                   const Tensor<2,dim> &src1,
-                   const Tensor<1,dim> &src2)
+template <int dim, typename Number>
+void outer_product (Tensor<3,dim,Number>       &dst,
+                   const Tensor<2,dim,Number> &src1,
+                   const Tensor<1,dim,Number> &src2)
 {
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
@@ -1270,10 +1275,10 @@ void outer_product (Tensor<3,dim>       &dst,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
-void outer_product (Tensor<1,dim>       &dst,
-                   const double         src1,
-                   const Tensor<1,dim> &src2)
+template <int dim, typename Number>
+void outer_product (Tensor<1,dim,Number>       &dst,
+                   const Number                src1,
+                   const Tensor<1,dim,Number> &src2)
 {
   for (unsigned int i=0; i<dim; ++i)
     dst[i] = src1 * src2[i];
@@ -1292,10 +1297,10 @@ void outer_product (Tensor<1,dim>       &dst,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
-void outer_product (Tensor<1,dim>       &dst,
-                   const Tensor<1,dim>  src1,
-                   const double         src2)
+template <int dim, typename Number>
+void outer_product (Tensor<1,dim,Number>       &dst,
+                   const Tensor<1,dim,Number>  src1,
+                   const Number         src2)
 {
   for (unsigned int i=0; i<dim; ++i)
     dst[i] = src1[i] * src2;
@@ -1314,13 +1319,13 @@ void outer_product (Tensor<1,dim>       &dst,
  * @relates Tensor
  * @author Guido Kanschat, 2001
  */
-template <int dim>
+template <int dim, typename Number>
 inline
 void
-cross_product (Tensor<1,dim>       &dst,
-              const Tensor<1,dim> &src)
+cross_product (Tensor<1,dim,Number>       &dst,
+              const Tensor<1,dim,Number> &src)
 {
-  Assert (dim==2, ExcImpossibleInDim(dim));
+  Assert (dim==2, ExcInternalError());
 
   dst[0] = src[1];
   dst[1] = -src[0];
@@ -1337,14 +1342,14 @@ cross_product (Tensor<1,dim>       &dst,
  * @relates Tensor
  * @author Guido Kanschat, 2001
  */
-template <int dim>
+template <int dim, typename Number>
 inline
 void
-cross_product (Tensor<1,dim>       &dst,
-              const Tensor<1,dim> &src1,
-              const Tensor<1,dim> &src2)
+cross_product (Tensor<1,dim,Number>       &dst,
+              const Tensor<1,dim,Number> &src1,
+              const Tensor<1,dim,Number> &src2)
 {
-  Assert (dim==3, ExcImpossibleInDim(dim));
+  Assert (dim==3, ExcInternalError());
 
   dst[0] = src1[1]*src2[2] - src1[2]*src2[1];
   dst[1] = src1[2]*src2[0] - src1[0]*src2[2];
@@ -1362,13 +1367,13 @@ cross_product (Tensor<1,dim>       &dst,
  * @relates Tensor
  * @author Wolfgang Bangerth, 2008
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double
-scalar_product (const Tensor<2,dim> &t1,
-               const Tensor<2,dim> &t2)
+Number
+scalar_product (const Tensor<2,dim,Number> &t1,
+               const Tensor<2,dim,Number> &t2)
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=0; j<dim; ++j)
       s += t1[i][j] * t2[i][j];
@@ -1384,9 +1389,9 @@ scalar_product (const Tensor<2,dim> &t1,
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
-template <int rank>
+template <int rank, typename Number>
 inline
-double determinant (const Tensor<rank,1> &t)
+Number determinant (const Tensor<rank,1,Number> &t)
 {
   return determinant(t[0]);
 }
@@ -1401,8 +1406,9 @@ double determinant (const Tensor<rank,1> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
+template <typename Number>
 inline
-double determinant (const Tensor<1,1> &t)
+Number determinant (const Tensor<1,1,Number> &t)
 {
   return t[0];
 }
@@ -1416,8 +1422,9 @@ double determinant (const Tensor<1,1> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
+template <typename Number>
 inline
-double determinant (const Tensor<2,1> &t)
+Number determinant (const Tensor<2,1,Number> &t)
 {
   return t[0][0];
 }
@@ -1430,8 +1437,9 @@ double determinant (const Tensor<2,1> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
+template <typename Number>
 inline
-double determinant (const Tensor<2,2> &t)
+Number determinant (const Tensor<2,2,Number> &t)
 {
   return ((t[0][0] * t[1][1]) -
          (t[1][0] * t[0][1]));
@@ -1444,8 +1452,9 @@ double determinant (const Tensor<2,2> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 1998
  */
+template <typename Number>
 inline
-double determinant (const Tensor<2,3> &t)
+Number determinant (const Tensor<2,3,Number> &t)
 {
                                   // get this using Maple:
                                   // with(linalg);
@@ -1469,9 +1478,9 @@ double determinant (const Tensor<2,3> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2009
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-double determinant (const Tensor<2,dim> &t)
+Number determinant (const Tensor<2,dim,Number> &t)
 {
                                   // compute the determinant using the
                                   // Laplace expansion of the
@@ -1481,7 +1490,7 @@ double determinant (const Tensor<2,dim> &t)
                                   //
                                   // for some algorithmic simplicity, we use
                                   // the expansion along the last row
-  double det = 0;
+  Number det = 0;
 
   for (unsigned int k=0; k<dim; ++k)
     {
@@ -1490,13 +1499,13 @@ double determinant (const Tensor<2,dim> &t)
        for (unsigned int j=0; j<dim-1; ++j)
          minor[i][j] = t[i][j<k ? j : j+1];
 
-      const double cofactor = std::pow (-1., static_cast<double>(k+1)) *
+      const Number cofactor = std::pow (-1., static_cast<Number>(k+1)) *
                              determinant (minor);
 
       det += t[dim-1][k] * cofactor;
     }
 
-  return std::pow (-1., static_cast<double>(dim)) * det;
+  return std::pow (-1., static_cast<Number>(dim)) * det;
 }
 
 
@@ -1508,10 +1517,10 @@ double determinant (const Tensor<2,dim> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2001
  */
-template <int dim>
-double trace (const Tensor<2,dim> &d)
+template <int dim, typename Number>
+Number trace (const Tensor<2,dim,Number> &d)
 {
-  double t=0;
+  Number t=0;
   for (unsigned int i=0; i<dim; ++i)
     t += d[i][i];
   return t;
@@ -1528,12 +1537,12 @@ double trace (const Tensor<2,dim> &d)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2000
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<2,dim>
-invert (const Tensor<2,dim> &t)
+Tensor<2,dim,Number>
+invert (const Tensor<2,dim,Number> &t)
 {
-  Tensor<2,dim> return_tensor;
+  Tensor<2,dim,Number> return_tensor;
   switch (dim)
     {
       case 1:
@@ -1544,7 +1553,7 @@ invert (const Tensor<2,dim> &t)
                                             // this is Maple output,
                                             // thus a bit unstructured
       {
-       const double t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]);
+       const Number t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]);
        return_tensor[0][0] = t[1][1]*t4;
        return_tensor[0][1] = -t[0][1]*t4;
        return_tensor[1][0] = -t[1][0]*t4;
@@ -1554,7 +1563,7 @@ invert (const Tensor<2,dim> &t)
 
       case 3:
       {
-       const double t4 = t[0][0]*t[1][1],
+       const Number t4 = t[0][0]*t[1][1],
                     t6 = t[0][0]*t[1][2],
                     t8 = t[0][1]*t[1][0],
                    t00 = t[0][2]*t[1][0],
@@ -1597,16 +1606,16 @@ invert (const Tensor<2,dim> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2002
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<2,dim>
-transpose (const Tensor<2,dim> &t)
+Tensor<2,dim,Number>
+transpose (const Tensor<2,dim,Number> &t)
 {
-  Tensor<2,dim> tt = t;
+  Tensor<2,dim,Number> tt = t;
   for (unsigned int i=0; i<dim; ++i)
     for (unsigned int j=i+1; j<dim; ++j)
       {
-        const double x = tt[i][j];
+        const Number x = tt[i][j];
         tt[i][j] = tt[j][i];
         tt[j][i] = x;
       };
@@ -1622,9 +1631,10 @@ transpose (const Tensor<2,dim> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2002
  */
+template <typename Number>
 inline
-Tensor<2,1>
-transpose (const Tensor<2,1> &t)
+Tensor<2,1,Number>
+transpose (const Tensor<2,1,Number> &t)
 {
   return t;
 }
@@ -1639,12 +1649,13 @@ transpose (const Tensor<2,1> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2002
  */
+template <typename Number>
 inline
-Tensor<2,2>
-transpose (const Tensor<2,2> &t)
+Tensor<2,2,Number>
+transpose (const Tensor<2,2,Number> &t)
 {
-  const double x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}};
-  return Tensor<2,2>(x);
+  const Number x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}};
+  return Tensor<2,2,Number>(x);
 }
 
 
@@ -1657,19 +1668,78 @@ transpose (const Tensor<2,2> &t)
  * @relates Tensor
  * @author Wolfgang Bangerth, 2002
  */
+template <typename Number>
 inline
-Tensor<2,3>
-transpose (const Tensor<2,3> &t)
+Tensor<2,3,Number>
+transpose (const Tensor<2,3,Number> &t)
 {
-  const double x[3][3] = {{t[0][0], t[1][0], t[2][0]},
+  const Number x[3][3] = {{t[0][0], t[1][0], t[2][0]},
                           {t[0][1], t[1][1], t[2][1]},
                           {t[0][2], t[1][2], t[2][2]}};
-  return Tensor<2,3>(x);
+  return Tensor<2,3,Number>(x);
 }
 
 #endif // DOXYGEN
 
 
+
+/**
+ * Multiplication of a tensor of general rank with a scalar Number
+ * from the right.
+ *
+ * @relates Tensor
+ */
+template <int rank, int dim, typename Number>
+inline
+Tensor<rank,dim,Number>
+operator * (const Tensor<rank,dim,Number> &t,
+           const Number                   factor)
+{
+  Tensor<rank,dim,Number> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Multiplication of a tensor of general rank with a scalar Number
+ * from the left.
+ *
+ * @relates Tensor
+ */
+template <int rank, int dim, typename Number>
+inline
+Tensor<rank,dim,Number>
+operator * (const Number                   factor,
+           const Tensor<rank,dim,Number> &t)
+{
+  Tensor<rank,dim,Number> tt = t;
+  tt *= factor;
+  return tt;
+}
+
+
+
+/**
+ * Division of a tensor of general rank by a scalar Number.
+ *
+ * @relates Tensor
+ */
+template <int rank, int dim, typename Number>
+inline
+Tensor<rank,dim,Number>
+operator / (const Tensor<rank,dim,Number> &t,
+           const Number                   factor)
+{
+  Tensor<rank,dim,Number> tt = t;
+  tt /= factor;
+  return tt;
+}
+
+
+
+
 /**
  * Multiplication of a tensor of general rank with a scalar double
  * from the right.
index 9a2e511bf891fec51c20d19437537613eb75d580..bbfc19eb9a71e2a25937bf144f033fdc7283c0f7 100644 (file)
 #define __deal2__tensor_base_h
 
 
-// single this file out from tensor.h, since we want to derive Point<dim>
-// from Tensor<1,dim>. However, the point class will not need all the
-// tensor stuff, so we don't want the whole tensor package to be included
-// everytime we use a point.
+// single this file out from tensor.h, since we want to derive
+// Point<dim,Number> from Tensor<1,dim,Number>. However, the point class will
+// not need all the tensor stuff, so we don't want the whole tensor package to
+// be included everytime we use a point.
 
 
 #include <deal.II/base/config.h>
@@ -35,33 +35,39 @@ DEAL_II_NAMESPACE_OPEN
 #endif
 
 template <typename number> class Vector;
-template <int dim> class Point;
+
+// 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 <int dim, typename Number=double> 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>;
+template <int rank_, int dim, typename Number=double> class Tensor;
+template <int dim, typename Number> class Tensor<0,dim,Number>;
+template <int dim, typename Number> class Tensor<1,dim,Number>;
 
 
 
 /**
- * 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 is a specialized version of the <tt>Tensor<rank,dim,Number></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
+ * objects of type Tensor@<spacedim-dim,dim,Number@>, 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
+ * class that acts as a scalar (i.e. @p Number) for all purposes but
  * is part of the Tensor template family.
  *
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 2009
  */
-template <int dim>
-class Tensor<0,dim>
+template <int dim, typename Number>
+class Tensor<0,dim,Number>
 {
   public:
                                     /**
@@ -91,11 +97,11 @@ class Tensor<0,dim>
 
                                     /**
                                      * Type of stored objects. This
-                                     * is a double for a rank 1 tensor.
+                                     * is a Number for a rank 1 tensor.
                                      */
 
-    typedef double value_type;
-    
+    typedef Number value_type;
+
                                     /**
                                      * Constructor. Set to zero.
                                      */
@@ -107,21 +113,21 @@ class Tensor<0,dim>
                                      * array.
                                      */
     Tensor (const value_type &initializer);
-    
+
                                     /**
                                      * Copy constructor.
                                      */
-    Tensor (const Tensor<0,dim> &);
+    Tensor (const Tensor<0,dim,Number> &);
 
                                     /**
-                                     * Conversion to double. Since
+                                     * Conversion to Number. Since
                                      * rank-0 tensors are scalars,
                                      * this is a natural operation.
                                      */
-    operator double () const;
+    operator Number () const;
 
                                     /**
-                                     * Conversion to double. Since
+                                     * Conversion to Number. Since
                                      * rank-0 tensors are scalars,
                                      * this is a natural operation.
                                      *
@@ -129,59 +135,59 @@ class Tensor<0,dim>
                                      * conversion operator that
                                      * returns a writable reference.
                                      */
-    operator double& ();
+    operator Number& ();
 
                                     /**
                                      * Assignment operator.
                                      */
-    Tensor<0,dim> & operator = (const Tensor<0,dim> &);
+    Tensor<0,dim,Number> & operator = (const Tensor<0,dim,Number> &);
 
                                     /**
                                      * Assignment operator.
                                      */
-    Tensor<0,dim> & operator = (const double d);
+    Tensor<0,dim,Number> & operator = (const Number d);
 
                                     /**
                                      * Test for equality of two
                                      * tensors.
                                      */
-    bool operator == (const Tensor<0,dim> &) const;
+    bool operator == (const Tensor<0,dim,Number> &) const;
 
                                     /**
                                      * Test for inequality of two
                                      * tensors.
                                      */
-    bool operator != (const Tensor<0,dim> &) const;
+    bool operator != (const Tensor<0,dim,Number> &) const;
 
                                     /**
                                      * Add another vector, i.e. move
                                      * this point by the given
                                      * offset.
                                      */
-    Tensor<0,dim> & operator += (const Tensor<0,dim> &);
-    
+    Tensor<0,dim,Number> & operator += (const Tensor<0,dim,Number> &);
+
                                     /**
                                      * Subtract another vector.
                                      */
-    Tensor<0,dim> & operator -= (const Tensor<0,dim> &);
+    Tensor<0,dim,Number> & operator -= (const Tensor<0,dim,Number> &);
 
                                     /**
                                      * Scale the vector by
                                      * <tt>factor</tt>, i.e. multiply all
                                      * coordinates by <tt>factor</tt>.
                                      */
-    Tensor<0,dim> & operator *= (const double factor);
+    Tensor<0,dim,Number> & operator *= (const Number factor);
 
                                     /**
                                      * Scale the vector by <tt>1/factor</tt>.
                                      */
-    Tensor<0,dim> & operator /= (const double factor);
+    Tensor<0,dim,Number> & operator /= (const Number factor);
 
                                     /**
                                      * Returns the scalar product of
                                      * two vectors.
                                      */
-    double          operator * (const Tensor<0,dim> &) const;
+    Number                 operator * (const Tensor<0,dim,Number> &) const;
 
                                     /**
                                      * Add two tensors. If possible,
@@ -189,7 +195,7 @@ class Tensor<0,dim>
                                      * since this does not need to
                                      * copy a point at least once.
                                      */
-    Tensor<0,dim>   operator + (const Tensor<0,dim> &) const;
+    Tensor<0,dim,Number>   operator + (const Tensor<0,dim,Number> &) const;
 
                                     /**
                                      * Subtract two tensors. If
@@ -198,13 +204,13 @@ class Tensor<0,dim>
                                      * need to copy a point at least
                                      * once.
                                      */
-    Tensor<0,dim>   operator - (const Tensor<0,dim> &) const;
+    Tensor<0,dim,Number>   operator - (const Tensor<0,dim,Number> &) const;
 
                                     /**
                                      * Tensor with inverted entries.
                                      */
-    Tensor<0,dim>   operator - () const;
-    
+    Tensor<0,dim,Number>   operator - () const;
+
                                      /**
                                       * Return the Frobenius-norm of a
                                       * tensor, i.e. the square root
@@ -215,7 +221,7 @@ class Tensor<0,dim>
                                       * <tt>l<sub>2</sub></tt> norm of
                                       * the vector.
                                       */
-    double norm () const;
+    Number norm () const;
 
                                      /**
                                       * Return the square of the
@@ -229,7 +235,7 @@ class Tensor<0,dim>
                                      * may also be useful in other
                                      * contexts.
                                       */
-    double norm_square () const;
+    Number norm_square () const;
 
                                     /**
                                      * Reset all values to zero.
@@ -266,9 +272,9 @@ class Tensor<0,dim>
                     << "dim must be positive, but was " << arg1);
 
                      /**
-                      * Read or write the data of this object to or 
+                      * Read or write the data of this object to or
                       * from a stream for the purpose of serialization
-                      */ 
+                      */
     template <class Archive>
     void serialize(Archive & ar, const unsigned int version);
 
@@ -276,33 +282,35 @@ class Tensor<0,dim>
                                     /**
                                      * The value of this scalar object.
                                      */
-    double value;
+    Number 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
- * provides the basis for the functionality needed for tensors of higher rank.
+ * This class is a specialized version of the <tt>Tensor<rank,dim,Number></tt>
+ * class.  It handles tensors with one index, i.e. vectors, of fixed dimension
+ * and provides the basis for the functionality needed for tensors of higher
+ * rank.
  *
  * Within deal.II, the distinction between this class and its derived class
  * <tt>Point</tt> is that we use the <tt>Point</tt> class mainly to denote the
  * points that make up geometric objects. As such, they have a small number of
  * additional operations over general tensors of rank 1 for which we use the
- * <tt>Tensor<1,dim></tt> class. In particular, there is a distance() function
- * to compute the Euclidian distance between two points in space.
+ * <tt>Tensor<1,dim,Number></tt> class. In particular, there is a distance()
+ * function to compute the Euclidian distance between two points in space.
  *
  * However, the <tt>Point</tt> class is really only used where the coordinates
  * of an object can be thought to possess the dimension of a length. For all
  * other uses, such as the gradient of a scalar function (which is a tensor of
  * rank 1, or vector, with as many elements as a point object, but with
- * different physical units), we use the <tt>Tensor<1,dim></tt> class.
+ * different physical units), we use the <tt>Tensor<1,dim,Number></tt> class.
  *
  * @ingroup geomprimitives
  * @author Wolfgang Bangerth, 1998-2005
  */
-template <int dim>
-class Tensor<1,dim>
+template <int dim,typename Number>
+class Tensor<1,dim,Number>
 {
   public:
                                     /**
@@ -332,11 +340,11 @@ class Tensor<1,dim>
 
                                     /**
                                      * Type of stored objects. This
-                                     * is a double for a rank 1 tensor.
+                                     * is a Number for a rank 1 tensor.
                                      */
 
-    typedef double value_type;
-    
+    typedef Number value_type;
+
                                     /**
                                      * Declare an array type which can
                                      * be used to initialize statically
@@ -348,8 +356,8 @@ class Tensor<1,dim>
                                      * likely to overflow memory
                                      * limits.
                                      */
-    typedef double array_type[(dim!=0) ? dim : 100000000];
-    
+    typedef Number array_type[(dim!=0) ? dim : 100000000];
+
                                     /**
                                      * Constructor. Initialize all entries
                                      * to zero if <tt>initialize==true</tt>; this
@@ -363,11 +371,11 @@ class Tensor<1,dim>
                                      * array.
                                      */
     Tensor (const array_type &initializer);
-    
+
                                     /**
                                      * Copy constructor.
                                      */
-    Tensor (const Tensor<1,dim> &);
+    Tensor (const Tensor<1,dim,Number> &);
 
                                     /**
                                      * Read access to the <tt>index</tt>th
@@ -379,7 +387,7 @@ class Tensor<1,dim>
                                      * operator for
                                      * backcompatibility.
                                      */
-    double   operator [] (const unsigned int index) const;
+    Number   operator [] (const unsigned int index) const;
 
                                     /**
                                      * Read and write access to the
@@ -391,12 +399,12 @@ class Tensor<1,dim>
                                      * operator for
                                      * backcompatibility.
                                      */
-    double & operator [] (const unsigned int index);
+    Number & operator [] (const unsigned int index);
 
                                     /**
                                      * Assignment operator.
                                      */
-    Tensor<1,dim> & operator = (const Tensor<1,dim> &);
+    Tensor<1,dim,Number> & operator = (const Tensor<1,dim,Number> &);
 
                                     /**
                                      * This operator assigns a scalar
@@ -410,49 +418,49 @@ class Tensor<1,dim>
                                      * all elements of the tensor to
                                      * zero.
                                      */
-    Tensor<1,dim> & operator = (const double d);
+    Tensor<1,dim,Number> & operator = (const Number d);
 
                                     /**
                                      * Test for equality of two
                                      * tensors.
                                      */
-    bool operator == (const Tensor<1,dim> &) const;
+    bool operator == (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      * Test for inequality of two
                                      * tensors.
                                      */
-    bool operator != (const Tensor<1,dim> &) const;
+    bool operator != (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      * Add another vector, i.e. move
                                      * this point by the given
                                      * offset.
                                      */
-    Tensor<1,dim> & operator += (const Tensor<1,dim> &);
-    
+    Tensor<1,dim,Number> & operator += (const Tensor<1,dim,Number> &);
+
                                     /**
                                      * Subtract another vector.
                                      */
-    Tensor<1,dim> & operator -= (const Tensor<1,dim> &);
+    Tensor<1,dim,Number> & operator -= (const Tensor<1,dim,Number> &);
 
                                     /**
                                      * Scale the vector by
                                      * <tt>factor</tt>, i.e. multiply all
                                      * coordinates by <tt>factor</tt>.
                                      */
-    Tensor<1,dim> & operator *= (const double factor);
+    Tensor<1,dim,Number> & operator *= (const Number factor);
 
                                     /**
                                      * Scale the vector by <tt>1/factor</tt>.
                                      */
-    Tensor<1,dim> & operator /= (const double factor);
+    Tensor<1,dim,Number> & operator /= (const Number factor);
 
                                     /**
                                      * Returns the scalar product of
                                      * two vectors.
                                      */
-    double          operator * (const Tensor<1,dim> &) const;
+    Number                 operator * (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      * Add two tensors. If possible,
@@ -460,7 +468,7 @@ class Tensor<1,dim>
                                      * since this does not need to
                                      * copy a point at least once.
                                      */
-    Tensor<1,dim>   operator + (const Tensor<1,dim> &) const;
+    Tensor<1,dim,Number>   operator + (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      * Subtract two tensors. If
@@ -469,13 +477,13 @@ class Tensor<1,dim>
                                      * need to copy a point at least
                                      * once.
                                      */
-    Tensor<1,dim>   operator - (const Tensor<1,dim> &) const;
+    Tensor<1,dim,Number>   operator - (const Tensor<1,dim,Number> &) const;
 
                                     /**
                                      * Tensor with inverted entries.
                                      */
-    Tensor<1,dim>   operator - () const;
-    
+    Tensor<1,dim,Number>   operator - () const;
+
                                      /**
                                       * Return the Frobenius-norm of a
                                       * tensor, i.e. the square root
@@ -486,7 +494,7 @@ class Tensor<1,dim>
                                       * <tt>l<sub>2</sub></tt> norm of
                                       * the vector.
                                       */
-    double norm () const;
+    Number norm () const;
 
                                      /**
                                       * Return the square of the
@@ -500,7 +508,7 @@ class Tensor<1,dim>
                                      * may also be useful in other
                                      * contexts.
                                       */
-    double norm_square () const;
+    Number norm_square () const;
 
                                     /**
                                      * Reset all values to zero.
@@ -531,7 +539,7 @@ class Tensor<1,dim>
                                      * usual in C++, the rightmost
                                      * index marches fastest.
                                      */
-    void unroll (Vector<double> &result) const;
+    void unroll (Vector<Number> &result) const;
 
                                     /**
                                      * Determine an estimate for
@@ -554,14 +562,14 @@ class Tensor<1,dim>
     DeclException1 (ExcDimTooSmall,
                     int,
                     << "dim must be positive, but was " << arg1);
-                    
+
                      /**
-                      * Read or write the data of this object to or 
+                      * Read or write the data of this object to or
                       * from a stream for the purpose of serialization
-                      */ 
+                      */
     template <class Archive>
     void serialize(Archive & ar, const unsigned int version);
-    
+
   private:
                                     /**
                                      * Store the values in a simple
@@ -572,7 +580,7 @@ class Tensor<1,dim>
                                      * constructor to disallow the
                                      * creation of such an object.
                                      */
-    double values[(dim!=0) ? (dim) : 1];
+    Number values[(dim!=0) ? (dim) : 1];
 
 #ifdef DEAL_II_TEMPLATE_SPEC_ACCESS_WORKAROUND
   public:
@@ -588,9 +596,9 @@ class Tensor<1,dim>
                                      * case, even if it should be
                                      * public for your compiler.
                                      */
-    void unroll_recursion (Vector<double> &result,
+    void unroll_recursion (Vector<Number> &result,
                           unsigned int   &start_index) const;
-    
+
   private:
                                     /**
                                      * Make the following classes
@@ -605,7 +613,7 @@ class Tensor<1,dim>
                                      * that seems to be impossible as
                                      * well.
                                      */
-    template <int otherrank, int otherdim>  friend class dealii::Tensor;
+    template <int otherrank, int otherdim, typename OtherNumber>  friend class dealii::Tensor;
 
                                     /**
                                      * Point is allowed access to
@@ -619,24 +627,24 @@ class Tensor<1,dim>
                                 /**
                                  *  Prints the value of this scalar.
                                  */
-template <int dim>
-std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p);
+template <int dim,typename Number>
+std::ostream & operator << (std::ostream &out, const Tensor<0,dim,Number> &p);
 
                                 /**
                                  *  Prints the values of this tensor in the
                                  *  form <tt>x1 x2 x3 etc</tt>.
                                  */
-template <int dim>
-std::ostream & operator << (std::ostream &out, const Tensor<1,dim> &p);
+template <int dim,typename Number>
+std::ostream & operator << (std::ostream &out, const Tensor<1,dim,Number> &p);
 
 
 #ifndef DOXYGEN
 
 /*---------------------------- Inline functions: Tensor<0,dim> ------------------------*/
 
-template <int dim>
+template <int dim,typename Number>
 inline
-Tensor<0,dim>::Tensor ()
+Tensor<0,dim,Number>::Tensor ()
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
@@ -645,9 +653,9 @@ Tensor<0,dim>::Tensor ()
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim>::Tensor (const value_type &initializer)
+Tensor<0,dim,Number>::Tensor (const value_type &initializer)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
@@ -656,39 +664,39 @@ Tensor<0,dim>::Tensor (const value_type &initializer)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim>::Tensor (const Tensor<0,dim> &p)
+Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,Number> &p)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
-  
+
   value = p.value;
 }
 
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim>::operator double () const
+Tensor<0,dim,Number>::operator Number () const
 {
   return value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim>::operator double & ()
+Tensor<0,dim,Number>::operator Number & ()
 {
   return value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p)
 {
   value = p.value;
   return *this;
@@ -696,9 +704,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator = (const double d)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Number d)
 {
   value = d;
   return *this;
@@ -706,27 +714,27 @@ Tensor<0,dim> & Tensor<0,dim>::operator = (const double d)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-bool Tensor<0,dim>::operator == (const Tensor<0,dim> &p) const
+bool Tensor<0,dim,Number>::operator == (const Tensor<0,dim,Number> &p) const
 {
   return (value == p.value);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-bool Tensor<0,dim>::operator != (const Tensor<0,dim> &p) const
+bool Tensor<0,dim,Number>::operator != (const Tensor<0,dim,Number> &p) const
 {
   return !((*this) == p);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator += (const Tensor<0,dim,Number> &p)
 {
   value += p.value;
   return *this;
@@ -734,9 +742,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Number> &p)
 {
   value -= p.value;
   return *this;
@@ -744,9 +752,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator *= (const Number s)
 {
   value *= s;
   return *this;
@@ -754,9 +762,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s)
+Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator /= (const Number s)
 {
   value /= s;
   return *this;
@@ -764,81 +772,81 @@ Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<0,dim>::operator * (const Tensor<0,dim> &p) const
+Number Tensor<0,dim,Number>::operator * (const Tensor<0,dim,Number> &p) const
 {
   return value*p.value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> Tensor<0,dim>::operator + (const Tensor<0,dim> &p) const
+Tensor<0,dim,Number> Tensor<0,dim,Number>::operator + (const Tensor<0,dim,Number> &p) const
 {
   return value+p.value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> Tensor<0,dim>::operator - (const Tensor<0,dim> &p) const
+Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - (const Tensor<0,dim,Number> &p) const
 {
   return value-p.value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<0,dim> Tensor<0,dim>::operator - () const
+Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const
 {
   return -value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<0,dim>::norm () const
+Number Tensor<0,dim,Number>::norm () const
 {
   return std::abs (value);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<0,dim>::norm_square () const
+Number Tensor<0,dim,Number>::norm_square () const
 {
   return value*value;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-void Tensor<0,dim>::clear ()
+void Tensor<0,dim,Number>::clear ()
 {
   value = 0;
 }
 
-template <int dim>
+template <int dim, typename Number>
 template <class Archive>
 inline
-void Tensor<0,dim>::serialize(Archive & ar, const unsigned int)
+void Tensor<0,dim,Number>::serialize(Archive & ar, const unsigned int)
 {
   ar & value;
 }
 
-/*---------------------------- Inline functions: Tensor<1,dim> ------------------------*/
+/*---------------------------- Inline functions: Tensor<1,dim,Number> ------------------------*/
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim>::Tensor (const bool initialize)
+Tensor<1,dim,Number>::Tensor (const bool initialize)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
@@ -849,9 +857,9 @@ Tensor<1,dim>::Tensor (const bool initialize)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim>::Tensor (const array_type &initializer)
+Tensor<1,dim,Number>::Tensor (const array_type &initializer)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
 
@@ -861,12 +869,12 @@ Tensor<1,dim>::Tensor (const array_type &initializer)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim>::Tensor (const Tensor<1,dim> &p)
+Tensor<1,dim,Number>::Tensor (const Tensor<1,dim,Number> &p)
 {
   Assert (dim>0, ExcDimTooSmall(dim));
-  
+
   for (unsigned int i=0; i<dim; ++i)
     values[i] = p.values[i];
 }
@@ -875,7 +883,7 @@ Tensor<1,dim>::Tensor (const Tensor<1,dim> &p)
 
 template <>
 inline
-Tensor<1,0>::Tensor (const Tensor<1,0> &)
+Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &)
 {
                                   // at some places in the library,
                                   // we have Point<0> for formal
@@ -891,9 +899,9 @@ Tensor<1,0>::Tensor (const Tensor<1,0> &)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<1,dim>::operator [] (const unsigned int index) const
+Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const
 {
   Assert (index<dim, ExcIndexRange (index, 0, dim));
   return values[index];
@@ -901,9 +909,9 @@ double Tensor<1,dim>::operator [] (const unsigned int index) const
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double & Tensor<1,dim>::operator [] (const unsigned int index)
+Number & Tensor<1,dim,Number>::operator [] (const unsigned int index)
 {
   Assert (index<dim, ExcIndexRange (index, 0, dim));
   return values[index];
@@ -913,7 +921,7 @@ double & Tensor<1,dim>::operator [] (const unsigned int index)
 
 template <>
 inline
-Tensor<1,0> & Tensor<1,0>::operator = (const Tensor<1,0> &)
+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
@@ -930,82 +938,56 @@ Tensor<1,0> & Tensor<1,0>::operator = (const Tensor<1,0> &)
 
 
 
-template <>
+template <int dim, typename Number>
 inline
-Tensor<1,1> & Tensor<1,1>::operator = (const Tensor<1,1> &p)
+Tensor<1,dim,Number> &
+Tensor<1,dim,Number>::operator = (const Tensor<1,dim,Number> &p)
 {
                                   // unroll by hand since this is a
                                   // frequently called function and
                                   // some compilers don't want to
                                   // always unroll the loop in the
                                   // general template
-  values[0] = p.values[0];
+  switch (dim)
+    {
+    case 1:
+      values[0] = p.values[0];
+      break;
+    case 2:
+      values[0] = p.values[0];
+      values[1] = p.values[1];
+      break;
+    case 3:
+      values[0] = p.values[0];
+      values[1] = p.values[1];
+      values[2] = p.values[2];
+      break;
+    default:
+      for (unsigned int i=0; i<dim; ++i)
+       values[i] = p.values[i];
+    }
   return *this;
 }
 
 
 
-template <>
+template <int dim, typename Number>
 inline
-Tensor<1,2> & Tensor<1,2>::operator = (const Tensor<1,2> &p)
+Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator = (const Number d)
 {
-                                  // unroll by hand since this is a
-                                  // frequently called function and
-                                  // some compilers don't want to
-                                  // always unroll the loop in the
-                                  // general template
-  values[0] = p.values[0];
-  values[1] = p.values[1];
-  return *this;
-}
-
-
+  Assert (d==Number(0), ExcMessage ("Only assignment with zero is allowed"));
 
-template <>
-inline
-Tensor<1,3> & Tensor<1,3>::operator = (const Tensor<1,3> &p)
-{
-                                  // unroll by hand since this is a
-                                  // frequently called function and
-                                  // some compilers don't want to
-                                  // always unroll the loop in the
-                                  // general template
-  values[0] = p.values[0];
-  values[1] = p.values[1];
-  values[2] = p.values[2];
-  return *this;
-}
-
-
-
-template <int dim>
-inline
-Tensor<1,dim> & Tensor<1,dim>::operator = (const Tensor<1,dim> &p)
-{
-  for (unsigned int i=0; i<dim; ++i)
-    values[i] = p.values[i];
-  return *this;
-}
-
-
-
-template <int dim>
-inline
-Tensor<1,dim> & Tensor<1,dim>::operator = (const double d)
-{
-  Assert (d==0, ExcMessage ("Only assignment with zero is allowed"));
-  
   for (unsigned int i=0; i<dim; ++i)
     values[i] = 0;
-  
+
   return *this;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-bool Tensor<1,dim>::operator == (const Tensor<1,dim> &p) const
+bool Tensor<1,dim,Number>::operator == (const Tensor<1,dim,Number> &p) const
 {
   for (unsigned int i=0; i<dim; ++i)
     if (values[i] != p.values[i])
@@ -1017,25 +999,25 @@ bool Tensor<1,dim>::operator == (const Tensor<1,dim> &p) const
 
 template <>
 inline
-bool Tensor<1,0>::operator == (const Tensor<1,0> &) const
+bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const
 {
   return true;
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-bool Tensor<1,dim>::operator != (const Tensor<1,dim> &p) const
+bool Tensor<1,dim,Number>::operator != (const Tensor<1,dim,Number> &p) const
 {
   return !((*this) == p);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim> & Tensor<1,dim>::operator += (const Tensor<1,dim> &p)
+Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator += (const Tensor<1,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] += p.values[i];
@@ -1044,9 +1026,9 @@ Tensor<1,dim> & Tensor<1,dim>::operator += (const Tensor<1,dim> &p)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim> & Tensor<1,dim>::operator -= (const Tensor<1,dim> &p)
+Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator -= (const Tensor<1,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] -= p.values[i];
@@ -1055,9 +1037,9 @@ Tensor<1,dim> & Tensor<1,dim>::operator -= (const Tensor<1,dim> &p)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim> & Tensor<1,dim>::operator *= (const double s)
+Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator *= (const Number s)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] *= s;
@@ -1066,9 +1048,9 @@ Tensor<1,dim> & Tensor<1,dim>::operator *= (const double s)
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-Tensor<1,dim> & Tensor<1,dim>::operator /= (const double s)
+Tensor<1,dim,Number> & Tensor<1,dim,Number>::operator /= (const Number s)
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] /= s;
@@ -1077,86 +1059,63 @@ Tensor<1,dim> & Tensor<1,dim>::operator /= (const double s)
 
 
 
-template <>
+template <int dim, typename Number>
 inline
-double Tensor<1,1>::operator * (const Tensor<1,1> &p) const
+Number
+Tensor<1,dim,Number>::operator * (const Tensor<1,dim,Number> &p) const
 {
                                   // unroll by hand since this is a
                                   // frequently called function and
                                   // some compilers don't want to
                                   // always unroll the loop in the
                                   // general template
-  return (values[0] * p.values[0]);
+  switch (dim)
+    {
+    case 1:
+      return (values[0] * p.values[0]);
+      break;
+    case 2:
+      return (values[0] * p.values[0] +
+             values[1] * p.values[1]);
+      break;
+    case 3:
+      return (values[0] * p.values[0] +
+             values[1] * p.values[1] +
+             values[2] * p.values[2]);
+      break;
+    default:
+      Number q=0;
+      for (unsigned int i=0; i<dim; ++i)
+       q += values[i] * p.values[i];
+      return q;
+    }
 }
 
 
 
-template <>
+template <int dim, typename Number>
 inline
-double Tensor<1,2>::operator * (const Tensor<1,2> &p) const
+Tensor<1,dim,Number> Tensor<1,dim,Number>::operator + (const Tensor<1,dim,Number> &p) const
 {
-                                  // unroll by hand since this is a
-                                  // frequently called function and
-                                  // some compilers don't want to
-                                  // always unroll the loop in the
-                                  // general template
-  return (values[0] * p.values[0] +
-         values[1] * p.values[1]);
+  return (Tensor<1,dim,Number>(*this) += p);
 }
 
 
 
-template <>
+template <int dim, typename Number>
 inline
-double Tensor<1,3>::operator * (const Tensor<1,3> &p) const
+Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - (const Tensor<1,dim,Number> &p) const
 {
-                                  // unroll by hand since this is a
-                                  // frequently called function and
-                                  // some compilers don't want to
-                                  // always unroll the loop in the
-                                  // general template
-  return (values[0] * p.values[0] +
-         values[1] * p.values[1] +
-         values[2] * p.values[2]);
+  return (Tensor<1,dim,Number>(*this) -= p);
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<1,dim>::operator * (const Tensor<1,dim> &p) const
+Tensor<1,dim,Number> Tensor<1,dim,Number>::operator - () const
 {
-  double q=0;
-  for (unsigned int i=0; i<dim; ++i)
-    q += values[i] * p.values[i];
-  return q;
-}
-
-
-
-template <int dim>
-inline
-Tensor<1,dim> Tensor<1,dim>::operator + (const Tensor<1,dim> &p) const
-{
-  return (Tensor<1,dim>(*this) += p);
-}
-
-
-
-template <int dim>
-inline
-Tensor<1,dim> Tensor<1,dim>::operator - (const Tensor<1,dim> &p) const
-{
-  return (Tensor<1,dim>(*this) -= p);
-}
-
-
-
-template <int dim>
-inline
-Tensor<1,dim> Tensor<1,dim>::operator - () const
-{
-  Tensor<1,dim> result;
+  Tensor<1,dim,Number> result;
   for (unsigned int i=0; i<dim; ++i)
     result.values[i] = -values[i];
   return result;
@@ -1164,20 +1123,20 @@ Tensor<1,dim> Tensor<1,dim>::operator - () const
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<1,dim>::norm () const
+Number Tensor<1,dim,Number>::norm () const
 {
   return std::sqrt (norm_square());
 }
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-double Tensor<1,dim>::norm_square () const
+Number Tensor<1,dim,Number>::norm_square () const
 {
-  double s = 0;
+  Number s = 0;
   for (unsigned int i=0; i<dim; ++i)
     s += values[i] * values[i];
 
@@ -1186,9 +1145,9 @@ double Tensor<1,dim>::norm_square () const
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
-void Tensor<1,dim>::clear ()
+void Tensor<1,dim,Number>::clear ()
 {
   for (unsigned int i=0; i<dim; ++i)
     values[i] = 0;
@@ -1196,19 +1155,19 @@ void Tensor<1,dim>::clear ()
 
 
 
-template <int dim>
+template <int dim, typename Number>
 inline
 std::size_t
-Tensor<1,dim>::memory_consumption ()
+Tensor<1,dim,Number>::memory_consumption ()
 {
-  return sizeof(Tensor<1,dim>);
+  return sizeof(Tensor<1,dim,Number>);
 }
 
 
-template <int dim>
+template <int dim, typename Number>
 template <class Archive>
 inline
-void Tensor<1,dim>::serialize(Archive & ar, const unsigned int)
+void Tensor<1,dim,Number>::serialize(Archive & ar, const unsigned int)
 {
   ar & values;
 }
@@ -1219,13 +1178,13 @@ void Tensor<1,dim>::serialize(Archive & ar, const unsigned int)
  * Output operator for tensors of rank 0. Since such tensors are
  * scalars, we simply print this one value.
  *
- * @relates Tensor<0,dim>
+ * @relates Tensor<0,dim,Number>
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p)
+std::ostream & operator << (std::ostream &out, const Tensor<0,dim,Number> &p)
 {
-  out << static_cast<double>(p);
+  out << static_cast<Number>(p);
   return out;
 }
 
@@ -1235,11 +1194,11 @@ std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p)
  * Output operator for tensors of rank 1. Print the elements
  * consecutively, with a space in between.
  *
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
  */
-template <int dim>
+template <int dim, typename Number>
 inline
-std::ostream & operator << (std::ostream &out, const Tensor<1,dim> &p)
+std::ostream & operator << (std::ostream &out, const Tensor<1,dim,Number> &p)
 {
   for (unsigned int i=0; i<dim-1; ++i)
     out << p[i] << ' ';
@@ -1255,10 +1214,10 @@ std::ostream & operator << (std::ostream &out, const Tensor<1,dim> &p)
  * implemented specialized from the general template in order to avoid
  * a compiler warning that the loop is empty.
  *
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
  */
 inline
-std::ostream & operator << (std::ostream &out, const Tensor<1,1> &p)
+std::ostream & operator << (std::ostream &out, const Tensor<1,1,double> &p)
 {
   out << p[0];
 
@@ -1267,10 +1226,67 @@ std::ostream & operator << (std::ostream &out, const Tensor<1,1> &p)
 
 
 
+/**
+ * Multiplication of a tensor of rank 1 with a scalar Number from the right.
+ *
+ * @relates Tensor<1,dim,Number>
+ */
+template <int dim, typename Number>
+inline
+Tensor<1,dim,Number>
+operator * (const Tensor<1,dim,Number> &t,
+           const Number                factor)
+{
+  Tensor<1,dim,Number> tt;
+  for (unsigned int d=0; d<dim; ++d)
+    tt[d] = t[d] * factor;
+  return tt;
+}
+
+
+
+/**
+ * Multiplication of a tensor of rank 1 with a scalar Number from the left.
+ *
+ * @relates Tensor<1,dim,Number>
+ */
+template <int dim, typename Number>
+inline
+Tensor<1,dim,Number>
+operator * (const Number                factor,
+           const Tensor<1,dim,Number> &t)
+{
+  Tensor<1,dim,Number> tt;
+  for (unsigned int d=0; d<dim; ++d)
+    tt[d] = t[d] * factor;
+  return tt;
+}
+
+
+
+/**
+ * Division of a tensor of rank 1 by a scalar Number.
+ *
+ * @relates Tensor<1,dim,Number>
+ */
+template <int dim, typename Number>
+inline
+Tensor<1,dim,Number>
+operator / (const Tensor<1,dim,Number> &t,
+           const Number                factor)
+{
+  Tensor<1,dim,Number> tt;
+  for (unsigned int d=0; d<dim; ++d)
+    tt[d] = t[d] / factor;
+  return tt;
+}
+
+
+
 /**
  * Multiplication of a tensor of rank 1 with a scalar double from the right.
  *
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
  */
 template <int dim>
 inline
@@ -1289,7 +1305,7 @@ operator * (const Tensor<1,dim> &t,
 /**
  * Multiplication of a tensor of rank 1 with a scalar double from the left.
  *
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
  */
 template <int dim>
 inline
@@ -1308,7 +1324,7 @@ operator * (const double         factor,
 /**
  * Division of a tensor of rank 1 by a scalar double.
  *
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
  */
 template <int dim>
 inline
@@ -1321,6 +1337,8 @@ operator / (const Tensor<1,dim> &t,
     tt[d] = t[d] / factor;
   return tt;
 }
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 2f9b90f2b3a496a771a09e388b69a11b7ccc69b7..61e71f093c1232e897dd75b7fa536a36c1a853f0 100644 (file)
@@ -17,6 +17,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/table.h>
 #include <deal.II/base/index_set.h>
+#include <deal.II/base/point.h>
 #include <deal.II/lac/constraint_matrix.h>
 #include <deal.II/dofs/function_map.h>
 #include <deal.II/dofs/dof_handler.h>
@@ -32,7 +33,6 @@ template<int dim, class T> class Table;
 class SparsityPattern;
 template <typename number> class Vector;
 template <int dim> class Function;
-template <int dim> class Point;
 template <int dim, int spacedim> class FiniteElement;
 template <int dim, int spacedim> class DoFHandler;
 namespace hp
index 9fc27fb6ac53b4bd082e5646f6a7f0f9722fcc45..e93b1ec0daa0818a6efdfade7ddd7e95ed4a551c 100644 (file)
@@ -17,6 +17,8 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/geometry_info.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/symmetric_tensor.h>
 
 #include <vector>
 
@@ -24,8 +26,6 @@ DEAL_II_NAMESPACE_OPEN
 
 template <typename number> class FullMatrix;
 template <typename number> class Vector;
-template <int rank, int dim> class SymmetricTensor;
-template <int rank, int dim> class Tensor;
 template <int dim> class Quadrature;
 template <int dim, int spacedim> class FiniteElement;
 template <int dim, int spacedim> class DoFHandler;
index 9c64767dda06c0c4ca64ee2acdafd73b106cf707..5f34bb8943ca46111a912a1e510129e19259e0fb 100644 (file)
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/smartpointer.h>
+#include <deal.II/base/point.h>
 #include <iostream>
 #include <vector>
 #include <string>
 
 DEAL_II_NAMESPACE_OPEN
 
-template <int dim> class Point;
 template <int dim, int space_dim> class Triangulation;
 template <int dim> class CellData;
 class SubCellData;
index 5061bbadeccb002dea3d2f3368b95e9bd4807299..426ed34b517480a6f22183d9719d1e5c3a48fda0 100644 (file)
@@ -16,6 +16,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/geometry_info.h>
+#include <deal.II/base/point.h>
 #include <deal.II/grid/tria_iterator_base.h>
 #include <deal.II/grid/tria_iterator_selector.h>
 
@@ -28,8 +29,6 @@ namespace std
 
 DEAL_II_NAMESPACE_OPEN
 
-template <int dim> class Point;
-
 template <int dim, int spacedim> class Triangulation;
 template <typename Accessor> class TriaRawIterator;
 template <typename Accessor> class TriaIterator;
index 94f3f70e19b9cff43e6bd1dbbb9d32ed49bcc950..1f61a511a8f8eb1d753e40675af4256b28700602 100644 (file)
@@ -18,6 +18,7 @@
 #include <deal.II/base/table.h>
 #include <deal.II/lac/exceptions.h>
 #include <deal.II/lac/identity_matrix.h>
+#include <deal.II/base/tensor.h>
 
 #include <vector>
 #include <iomanip>
@@ -29,7 +30,6 @@ DEAL_II_NAMESPACE_OPEN
 // forward declarations
 template <typename number> class Vector;
 template <typename number> class LAPACKFullMatrix;
-template <int rank, int dim> class Tensor;
 
 
 /*! @addtogroup Matrix1
index f1b4532938cd0955c7ed25e5f4b306e74411e516..3b74811600e977982812057f02ad39eaf66a41ae 100644 (file)
@@ -16,6 +16,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/point.h>
 #include <deal.II/dofs/function_map.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/hp/mapping_collection.h>
@@ -26,7 +27,6 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-template <int dim> class Point;
 template <int dim> class Function;
 template <int dim> class FunctionMap;
 template <int dim> class Quadrature;
@@ -1173,8 +1173,8 @@ class VectorTools
      const unsigned char boundary_component,
      ConstraintMatrix& constraints,
      const hp::MappingCollection<dim, dim>& mapping_collection = hp::StaticMappingQ1<dim>::mapping_collection);
-   
-   
+
+
                                     /**
                                      * Compute the constraints that
                                      * correspond to boundary conditions of
index e142d3be6b4c495415bd88706c915eee7ae2d80b..2a8f0043b4a56caa8afbde59deb15bbb2460fa63 100644 (file)
@@ -14,8 +14,8 @@
 
 #include <deal.II/base/config.h>
 #include <cmath>
-#include <limits>
 #include <complex>
+#include <limits>
 
 DEAL_II_NAMESPACE_OPEN
 
index c4a895f6100e9b33f502760bb8262baddd7f1764..3b08e510bec5a535b15bcc6244a30a5e26a2beb8 100644 (file)
@@ -13,7 +13,6 @@
 
 #include <deal.II/base/function_parser.h>
 #include <deal.II/base/utilities.h>
-#include <deal.II/base/point.h>
 #include <deal.II/lac/vector.h>
 
 #ifndef DEAL_II_DISABLE_PARSER
@@ -38,7 +37,7 @@ FunctionParser<dim>::FunctionParser(const unsigned int n_components,
                 :
                 Function<dim>(n_components, initial_time),
                 fp (0)
-{ 
+{
   fp = new fparser::FunctionParser[n_components];
 }
 
@@ -64,7 +63,7 @@ void FunctionParser<dim>::initialize (const std::string                   &varia
                 constants,
                 std::map< std::string, double >(),
                 time_dependent,
-                use_degrees);  
+                use_degrees);
 }
 
 
@@ -104,7 +103,7 @@ void FunctionParser<dim>::initialize (const std::string   &variables,
         }
 
 
-                       
+
                                       // Add the various constants to
                                       // the parser.
       std::map< std::string, double >::const_iterator
@@ -115,21 +114,21 @@ void FunctionParser<dim>::initialize (const std::string   &variables,
          const bool success = fp[i].AddConstant(constant->first, constant->second);
          AssertThrow (success, ExcMessage("Invalid Constant Name"));
        }
-      
-      const int ret_value = fp[i].Parse(expressions[i], 
-                                       variables, 
+
+      const int ret_value = fp[i].Parse(expressions[i],
+                                       variables,
                                        use_degrees);
-      AssertThrow (ret_value == -1, 
+      AssertThrow (ret_value == -1,
                   ExcParseError(ret_value, fp[i].ErrorMsg()));
-      
+
                                       // The fact that the parser did
                                       // not throw an error does not
                                       // mean that everything went
                                       // ok... we can still have
                                       // problems with the number of
                                       // variables...
-    } 
-    
+    }
+
                                   // Now we define how many variables
                                   // we expect to read in.  We
                                   // distinguish between two cases:
@@ -143,16 +142,16 @@ void FunctionParser<dim>::initialize (const std::string   &variables,
                                   // parsed the variables string, if
                                   // none of this is the case, then
                                   // an exception is thrown.
-  if (time_dependent) 
+  if (time_dependent)
     n_vars = dim+1;
-  else 
+  else
     n_vars = dim;
-    
+
                                   // Let's check if the number of
                                   // variables is correct...
-  AssertThrow (n_vars == fp[0].NVars(), 
+  AssertThrow (n_vars == fp[0].NVars(),
               ExcDimensionMismatch(n_vars,fp[0].NVars()));
-    
+
                                   // Now set the initialization bit.
   initialized = true;
 }
@@ -172,7 +171,7 @@ void FunctionParser<dim>::initialize (const std::string &variables,
               Utilities::split_string_list (expression, ';'),
               constants,
               time_dependent,
-              use_degrees);  
+              use_degrees);
 }
 
 
@@ -199,19 +198,19 @@ void FunctionParser<dim>::initialize (const std::string &variables,
 
 template <int dim>
 double FunctionParser<dim>::value (const Point<dim>  &p,
-                                  const unsigned int component) const 
+                                  const unsigned int component) const
 {
   Assert (initialized==true, ExcNotInitialized());
   Assert (component < this->n_components,
          ExcIndexRange(component, 0, this->n_components));
-  
+
                                   // Statically allocate dim+1
                                   // double variables.
   double vars[dim+1];
-  
+
   for (unsigned int i=0; i<dim; ++i)
     vars[i] = p(i);
-  
+
                                   // We need the time variable only
                                   // if the number of variables is
                                   // different from the dimension. In
@@ -220,9 +219,9 @@ double FunctionParser<dim>::value (const Point<dim>  &p,
                                   // have already been thrown
   if (dim != n_vars)
     vars[dim] = this->get_time();
-  
+
   double my_value = fp[component].Eval((double*)vars);
-  
+
   AssertThrow (fp[component].EvalError() == 0,
               ExcMessage(fp[component].ErrorMsg()));
   return my_value;
@@ -232,19 +231,19 @@ double FunctionParser<dim>::value (const Point<dim>  &p,
 
 template <int dim>
 void FunctionParser<dim>::vector_value (const Point<dim> &p,
-                                       Vector<double>   &values) const 
+                                       Vector<double>   &values) const
 {
   Assert (initialized==true, ExcNotInitialized());
   Assert (values.size() == this->n_components,
          ExcDimensionMismatch (values.size(), this->n_components));
-  
+
                                   // Statically allocates dim+1
                                   // double variables.
   double vars[dim+1];
-    
+
   for(unsigned int i=0; i<dim; ++i)
     vars[i] = p(i);
-    
+
                                   // We need the time variable only
                                   // if the number of variables is
                                   // different from the dimension. In
@@ -253,14 +252,14 @@ void FunctionParser<dim>::vector_value (const Point<dim> &p,
                                   // have already been thrown
   if(dim != n_vars)
     vars[dim] = this->get_time();
-    
+
   for(unsigned int component = 0; component < this->n_components;
       ++component)
     {
       values(component) = fp[component].Eval((double*)vars);
       AssertThrow(fp[component].EvalError() == 0,
                  ExcMessage(fp[component].ErrorMsg()));
-    } 
+    }
 }
 
 #else
@@ -292,7 +291,7 @@ FunctionParser<dim>::initialize(const std::string &,
 
 template <int dim>
 double FunctionParser<dim>::value (
-  const Point<dim> &, unsigned int) const 
+  const Point<dim> &, unsigned int) const
 {
   Assert(false, ExcDisabled("parser"));
   return 0.;
@@ -301,7 +300,7 @@ double FunctionParser<dim>::value (
 
 template <int dim>
 void FunctionParser<dim>::vector_value (
-  const Point<dim>&, Vector<double>&) const 
+  const Point<dim>&, Vector<double>&) const
 {
   Assert(false, ExcDisabled("parser"));
 }
index c9b009af88967a0274205fa31d43e15adcd3b2f7..30ef5f61bbaa08053ccc9f3cf4cd1f15c90b2345 100644 (file)
@@ -18,10 +18,10 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template <>
-SymmetricTensor<4,3>
-invert (const SymmetricTensor<4,3> &t)
+SymmetricTensor<4,3,double>
+invert<3,double> (const SymmetricTensor<4,3,double> &t)
 {
-  SymmetricTensor<4,3> tmp = t;
+  SymmetricTensor<4,3,double> tmp = t;
 
                                    // this function follows the exact same
                                    // scheme as the 2d case, except that
@@ -45,7 +45,7 @@ invert (const SymmetricTensor<4,3> &t)
   for (unsigned int i=0; i<N; ++i)
     diagonal_sum += std::fabs(tmp.data[i][i]);
   const double typical_diagonal_element = diagonal_sum/N;
-  
+
   unsigned int p[N];
   for (unsigned int i=0; i<N; ++i)
     p[i] = i;
@@ -68,7 +68,7 @@ invert (const SymmetricTensor<4,3> &t)
                                       // too small
       Assert(max > 1.e-16*typical_diagonal_element,
             ExcMessage("This tensor seems to be noninvertible"));
-      
+
                                       // row interchange
       if (r>j)
        {
@@ -120,7 +120,7 @@ invert (const SymmetricTensor<4,3> &t)
   for (unsigned int i=3; i<6; ++i)
     for (unsigned int j=3; j<6; ++j)
       tmp.data[i][j] /= 4;
-  
+
   return tmp;
 }
 
index 16cbecd6efc13cb02fa71c1d5c64a0c6ab664c45..5ade4a855413cdc5613cc5f54804fd7b1c5339be 100644 (file)
 //---------------------------------------------------------------------------
 
 
-#include <deal.II/base/tensor.h>
+#include <base/tensor.h>
 #include <cmath>
-#include <deal.II/lac/vector.h>
+#include <lac/vector.h>
 
 DEAL_II_NAMESPACE_OPEN
 
 
 // storage for static variables
-template <int dim>
-const unsigned int Tensor<1,dim>::dimension;
+template <int dim, typename Number>
+const unsigned int Tensor<1,dim,Number>::dimension;
 
-template <int rank, int dim>
-const unsigned int Tensor<rank,dim>::dimension;
+template <int rank, int dim, typename Number>
+const unsigned int Tensor<rank,dim,Number>::dimension;
 
 
-template <int dim>
+template <int dim, typename Number>
 void
-Tensor<1,dim>::unroll (Vector<double> &result) const
+Tensor<1,dim,Number>::unroll (Vector<Number> &result) const
 {
   Assert (result.size()==dim,
           ExcDimensionMismatch(dim, result.size()));
@@ -40,12 +40,12 @@ Tensor<1,dim>::unroll (Vector<double> &result) const
 
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 void
-Tensor<rank_, dim>::unroll (Vector<double> &result) const
+Tensor<rank_, dim, Number>::unroll (Vector<Number> &result) const
 {
-  Assert(result.size()==std::pow(static_cast<double>(dim),rank_),
-        ExcDimensionMismatch(static_cast<unsigned int>(std::pow(static_cast<double>(dim),rank_)),
+  Assert(result.size()==std::pow(static_cast<Number>(dim),rank_),
+        ExcDimensionMismatch(static_cast<unsigned int>(std::pow(static_cast<Number>(dim),rank_)),
                               result.size()));
 
   unsigned index = 0;
@@ -54,40 +54,40 @@ Tensor<rank_, dim>::unroll (Vector<double> &result) const
 
 
 
-template <int rank_, int dim>
+template <int rank_, int dim, typename Number>
 void
-Tensor<rank_, dim>::unroll_recursion (Vector<double> &result,
-                                      unsigned int   &index) const
+Tensor<rank_, dim, Number>::unroll_recursion (Vector<Number> &result,
+                                             unsigned int   &index) const
 {
   for (unsigned i=0; i<dim; ++i)
     {
       operator[](i).unroll_recursion(result, index);
-    }    
+    }
 }
 
 
 
-template<int dim>
+template<int dim, typename Number>
 void
-Tensor<1,dim>::unroll_recursion (Vector<double> &result,
-                                 unsigned int   &index) const
+Tensor<1,dim,Number>::unroll_recursion (Vector<Number> &result,
+                                       unsigned int   &index) const
 {
   for (unsigned i=0; i<dim; ++i)
-    result(index++) = operator[](i);  
+    result(index++) = operator[](i);
 }
 
 
-template class Tensor<1, 1>;
-template class Tensor<1, 2>;
-template class Tensor<1, 3>;
-template class Tensor<2, 1>;
-template class Tensor<2, 2>;
-template class Tensor<2, 3>;
-template class Tensor<3, 1>;
-template class Tensor<3, 2>;
-template class Tensor<3, 3>;
-template class Tensor<4, 1>;
-template class Tensor<4, 2>;
-template class Tensor<4, 3>;
+template class Tensor<1, 1, double>;
+template class Tensor<1, 2, double>;
+template class Tensor<1, 3, double>;
+template class Tensor<2, 1, double>;
+template class Tensor<2, 2, double>;
+template class Tensor<2, 3, double>;
+template class Tensor<3, 1, double>;
+template class Tensor<3, 2, double>;
+template class Tensor<3, 3, double>;
+template class Tensor<4, 1, double>;
+template class Tensor<4, 2, double>;
+template class Tensor<4, 3, double>;
 
 DEAL_II_NAMESPACE_CLOSE
index fc2fa9901f267380e5f54f2db8f7d534d46ea953..9b09cc5103aad1deceb8461802f02f733c8211c8 100644 (file)
@@ -83,6 +83,27 @@ namespace FETools
   {
     return new FE_Q<3>(quad);
   }
+
+
+                               // Specializations for FE_DGQArbitraryNodes.
+  template <>
+  FiniteElement<1, 1>*
+  FEFactory<FE_DGQ<1> >::get (const Quadrature<1> &quad) const
+  {
+    return new FE_DGQArbitraryNodes<1>(quad);
+  }
+  template <>
+  FiniteElement<2, 2>*
+  FEFactory<FE_DGQ<2> >::get (const Quadrature<1> &quad) const
+  {
+    return new FE_DGQArbitraryNodes<2>(quad);
+  }
+  template <>
+  FiniteElement<3, 3>*
+  FEFactory<FE_DGQ<3> >::get (const Quadrature<1> &quad) const
+  {
+    return new FE_DGQArbitraryNodes<3>(quad);
+  }
 }
 
 namespace
index f612f60f3120e815a872302abaf4d5ed81578f9b..ad814bd73be1124418ffe54ab29aa81e8e3c0e01 100644 (file)
@@ -1118,7 +1118,7 @@ namespace internal
                                                       // FEFaceValues in 1d
                      Assert (false, ExcNotImplemented());
                      cross_product (boundary_forms[point],
-                                    (face_no == 0 ? 1 : -1) * cell_normal);
+                                    (face_no == 0 ? 1. : -1.) * cell_normal);
                    }
                  else if ( (dim==2) && (spacedim==3) )
                    cross_product (boundary_forms[point],
index c3175e6d4c60194b3997b59064a3142597330aca..3b9e8959f556f39cfda143317a4e5ed79bb7b8fa 100644 (file)
@@ -1607,7 +1607,7 @@ namespace hp
                        "this class does not currently support this."));
 
     create_active_fe_table ();
-    
+
     tria_listeners.push_back
     (tria.signals.pre_refinement
      .connect (std_cxx1x::bind (&DoFHandler<dim,spacedim>::pre_refinement_action,
@@ -1631,7 +1631,7 @@ namespace hp
     for (unsigned int i=0; i<tria_listeners.size(); ++i)
       tria_listeners[i].disconnect ();
     tria_listeners.clear ();
-    
+
                                      // ...and release allocated memory
     clear ();
   }
index 0db68343a901d91a33204e5fc37c99df899bfb99..60e34866f8e84df30007f27a25db7ab737f75d9d 100644 (file)
@@ -13,8 +13,6 @@
 
 
 #include <deal.II/base/logstream.h>
-#include <deal.II/base/tensor.h>
-#include <deal.II/base/tensor_base.h>
 #include <deal.II/lac/full_matrix.templates.h>
 
 DEAL_II_NAMESPACE_OPEN

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.