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>
</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)
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
<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)
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
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
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:
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
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,
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
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>
#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;
/**
* template programming.
*/
static const unsigned int dimension = dim;
-
+
/**
* Number of vector components.
*/
*/
Function (const unsigned int n_components = 1,
const double initial_time = 0.0);
-
+
/**
* Virtual destructor; absolutely
* necessary in this case.
* to.
*/
Function & operator= (const Function &f);
-
+
/**
* Return the value of the
* function at the given
*/
virtual void vector_value (const Point<dim> &p,
Vector<double> &values) const;
-
+
/**
* Set <tt>values</tt> to the point
* values of the specified
*/
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
*/
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
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
*/
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.
*
*/
virtual void vector_laplacian (const Point<dim> &p,
Vector<double> &values) const;
-
+
/**
* Compute the Laplacian of one
* component at a set of points.
* 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
*/
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
*/
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>,
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.
*
*/
ConstantFunction (const double value,
const unsigned int n_components = 1);
-
+
/**
* Virtual destructor; absolutely
* necessary in this case.
* 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
#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>
DEAL_II_NAMESPACE_OPEN
-template <int> class Point;
-template <int, int> class Tensor;
-template <int dim> class Tensor<1,dim>;
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.
* 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
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,
"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
*
* 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
* 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
* 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()
* 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
* 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
const ConstMap &units,
const bool time_dependent = false,
const bool use_degrees = false);
-
+
/**
* Initialize the function. Same as
* above, but accepts a string rather
/** @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
* in the initialize() methods.
*/
bool initialized;
-
+
/**
* Number of variables. If this
* is also a function of time,
return "";
}
-
+
DEAL_II_NAMESPACE_CLOSE
#endif
-
+
* @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:
/**
* an origin.
*/
Point ();
-
+
/**
* Constructor. Initialize all
* entries to zero if
* 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
* 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
* 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
* 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
* 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
* 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
* 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
* 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);
};
#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;
-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;
-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;
}
-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];
-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];
-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;
-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
/*------------------------------- 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] << ' ';
* 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];
#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];
* <i>Q<sub>k</sub></i> contained.
*/
PolynomialsRaviartThomas (const unsigned int k);
-
+
/**
* Computes the value and the
* first and second derivatives
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
* which is <tt>RaviartThomas</tt>.
*/
std::string name () const;
-
+
/**
* Return the number of
* polynomials in the space
* 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
};
+
template <int dim>
inline unsigned int
PolynomialsRaviartThomas<dim>::n() const
}
+
template <int dim>
inline unsigned int
PolynomialsRaviartThomas<dim>::degree() const
}
+
template <int dim>
inline std::string
PolynomialsRaviartThomas<dim>::name() const
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
* 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;
};
* 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;
};
* 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
* 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;
};
* 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
*/
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
* 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;
};
* 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;
/**
*
* 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;
};
/**
* 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;
};
*
* @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:
* 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:
/**
/**
* 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:
/**
// 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
};
*
* @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:
/**
// 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
};
}
* @ingroup geomprimitives
* @author Wolfgang Bangerth, 2005
*/
-template <int rank, int dim>
+template <int rank, int dim, typename Number>
class SymmetricTensor
{
public:
* 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;
/**
* 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
* 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.
* all elements of the tensor to
* zero.
*/
- SymmetricTensor & operator = (const double d);
+ SymmetricTensor & operator = (const Number d);
/**
* Convert the present symmetric tensor
* 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.
* 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
* 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
* 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
* 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;
/**
* 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);
/**
* specify the entire set of
* indices.
*/
- double
+ Number
operator [] (const TableIndices<rank> &indices) const;
/**
* specify the entire set of
* indices.
*/
- double &
+ Number &
operator [] (const TableIndices<rank> &indices);
/**
* just one of them, although they are
* equal for symmetric tensors).
*/
- double norm () const;
+ Number norm () const;
/**
* Tensors can be unrolled by
* properties of the base tensor.
*/
typedef
- internal::SymmetricTensorAccessors::StorageType<rank,dim>
+ internal::SymmetricTensorAccessors::StorageType<rank,dim,Number>
base_tensor_descriptor;
/**
/**
* 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> &);
};
{
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)
:
- 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)
:
- 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));
}
-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;
-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"));
-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;
-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;
-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;
-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;
-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;
-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;
-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;
-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
// 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
// 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
// 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
// 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;
}
* 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.
* @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;
+ }
}
* @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);
}
* @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;
* @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);
}
* @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;
}
* @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];
}
* @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]);
* @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;
}
* @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;
/**
- * 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;
}
* @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)
// 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;
* @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)
// 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;
* 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
* @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
// 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;
}
* @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
/**
* <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)
* @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);
}
* @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);
}
* @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
inline
SymmetricTensor<rank,dim>
operator * (const SymmetricTensor<rank,dim> &t,
- const double factor)
+ const double factor)
{
SymmetricTensor<rank,dim> tt = t;
tt *= factor;
/**
- * 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
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;
/**
- * 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
*/
inline
SymmetricTensor<rank,dim>
operator / (const SymmetricTensor<rank,dim> &t,
- const double factor)
+ const double factor)
{
SymmetricTensor<rank,dim> tt = t;
tt /= factor;
* @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];
}
* @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];
}
*
* @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;
* @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;
* @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;
* @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;
* @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];
*
* @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)
*
* @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)
#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
* 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:
* 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
* 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
/**
* 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
* 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
* 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,
* 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
* may also be useful in other
* contexts.
*/
- double norm_square () const;
+ Number norm_square () const;
/**
* Fill a vector with all tensor elements.
* usual in C++, the rightmost
* index of the tensor marches fastest.
*/
- void unroll (Vector<double> & result) const;
+ void unroll (Vector<Number> & result) const;
/**
* 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
// 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;
};
#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
}
-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));
}
-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));
}
-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];
}
-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"));
}
-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;
}
-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];
}
-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];
}
-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;
}
-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;
}
-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];
}
-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];
}
-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];
}
-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();
}
-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;
}
*
* @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] << ' ';
* @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];
* @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);
}
* @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)
* @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];
* @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)
* @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];
* @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)
* @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)
* @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 ();
default:
Assert (false,
- (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+ (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
};
break;
case 2:
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)));
};
}
* @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 ();
default:
Assert (false,
- (typename Tensor<2,dim>::ExcInvalidTensorIndex (index1)));
+ (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1)));
};
}
* @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)
*
* @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 ();
break;
default:
Assert (false,
- (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+ (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
}
break;
break;
default:
Assert (false,
- (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2)));
+ (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2)));
}
break;
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)));
}
}
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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)
* @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];
* @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;
* @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];
* @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];
* @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];
* @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]);
}
* @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];
}
* @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];
}
* @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]));
* @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);
* @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
//
// 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)
{
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;
}
* @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;
* @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:
// 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;
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],
* @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;
};
* @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;
}
* @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);
}
* @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.
#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>
#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:
/**
/**
* 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.
*/
* 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.
*
* 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,
* 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
* 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
* <tt>l<sub>2</sub></tt> norm of
* the vector.
*/
- double norm () const;
+ Number norm () const;
/**
* Return the square of the
* may also be useful in other
* contexts.
*/
- double norm_square () const;
+ Number norm_square () const;
/**
* Reset all values to zero.
<< "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);
/**
* 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:
/**
/**
* 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
* 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
* 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
* operator for
* backcompatibility.
*/
- double operator [] (const unsigned int index) const;
+ Number operator [] (const unsigned int index) const;
/**
* Read and write access to the
* 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
* 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,
* 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
* 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
* <tt>l<sub>2</sub></tt> norm of
* the vector.
*/
- double norm () const;
+ Number norm () const;
/**
* Return the square of the
* may also be useful in other
* contexts.
*/
- double norm_square () const;
+ Number norm_square () const;
/**
* Reset all values to zero.
* usual in C++, the rightmost
* index marches fastest.
*/
- void unroll (Vector<double> &result) const;
+ void unroll (Vector<Number> &result) const;
/**
* Determine an estimate for
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
* 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:
* 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
* 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
/**
* 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));
-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));
-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;
-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;
-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;
-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;
-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;
-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;
-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));
-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));
-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];
}
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
-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];
-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];
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
-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])
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];
-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];
-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;
-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;
-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;
-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];
-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;
-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;
}
* 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;
}
* 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] << ' ';
* 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];
+/**
+ * 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
/**
* 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
/**
* Division of a tensor of rank 1 by a scalar double.
*
- * @relates Tensor<1,dim>
+ * @relates Tensor<1,dim,Number>
*/
template <int dim>
inline
tt[d] = t[d] / factor;
return tt;
}
+
+
DEAL_II_NAMESPACE_CLOSE
#endif
#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>
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
#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>
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;
#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;
#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>
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;
#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>
// forward declarations
template <typename number> class Vector;
template <typename number> class LAPACKFullMatrix;
-template <int rank, int dim> class Tensor;
/*! @addtogroup Matrix1
#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>
DEAL_II_NAMESPACE_OPEN
-template <int dim> class Point;
template <int dim> class Function;
template <int dim> class FunctionMap;
template <int dim> class Quadrature;
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
#include <deal.II/base/config.h>
#include <cmath>
-#include <limits>
#include <complex>
+#include <limits>
DEAL_II_NAMESPACE_OPEN
#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
:
Function<dim>(n_components, initial_time),
fp (0)
-{
+{
fp = new fparser::FunctionParser[n_components];
}
constants,
std::map< std::string, double >(),
time_dependent,
- use_degrees);
+ use_degrees);
}
}
-
+
// Add the various constants to
// the parser.
std::map< std::string, double >::const_iterator
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:
// 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;
}
Utilities::split_string_list (expression, ';'),
constants,
time_dependent,
- use_degrees);
+ use_degrees);
}
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
// 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;
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
// 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
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.;
template <int dim>
void FunctionParser<dim>::vector_value (
- const Point<dim>&, Vector<double>&) const
+ const Point<dim>&, Vector<double>&) const
{
Assert(false, ExcDisabled("parser"));
}
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
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;
// too small
Assert(max > 1.e-16*typical_diagonal_element,
ExcMessage("This tensor seems to be noninvertible"));
-
+
// row interchange
if (r>j)
{
for (unsigned int i=3; i<6; ++i)
for (unsigned int j=3; j<6; ++j)
tmp.data[i][j] /= 4;
-
+
return tmp;
}
//---------------------------------------------------------------------------
-#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()));
-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;
-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
{
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
// 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],
"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,
for (unsigned int i=0; i<tria_listeners.size(); ++i)
tria_listeners[i].disconnect ();
tria_listeners.clear ();
-
+
// ...and release allocated memory
clear ();
}
#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
projects / dealii-svn.git / commitdiff