From: kronbichler Date: Tue, 2 Aug 2011 09:47:19 +0000 (+0000) Subject: Extend classes Tensor, SymmetricTensor and Point to work with templates for the numbe... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=03d659809da981ef5bfa4ab981b230debb2574ef;p=dealii-svn.git Extend classes Tensor, SymmetricTensor and Point to work with templates for the number type. Default template value is double, so all old code remains valid, except for forward declarations without the number type. git-svn-id: https://svn.dealii.org/trunk@23989 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/doxygen/headers/coding_conventions.h b/deal.II/doc/doxygen/headers/coding_conventions.h index 7265d69103..9136b101b5 100644 --- a/deal.II/doc/doxygen/headers/coding_conventions.h +++ b/deal.II/doc/doxygen/headers/coding_conventions.h @@ -102,6 +102,10 @@ code is written, without having to look up the exact definition of something. and DoFHandler (named cell_iterator, active_line_iterator, etc) to make the connection to the STL classes clear. +
  • For classes with multiple template arguments, the dimension is usually + put before the data type specifier, i.e., we use Point and not + Point. +
  • Each class has to have at least 200 pages of documentation ;-)
    • diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index be4b64a31e..08d6c50b3e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -21,6 +21,19 @@ inconvenience this causes.

      +
    1. 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 template + class Point that might be present in user-defined header +files. Now forward declarations need to specify the type as well, i.e., +template class Point. However, +nothing changes if the full declarations in deal.II/base/tensor.h, +deal.II/base/symmetric_tensor.h and deal.II/base/point.h +are included. +
      +(Martin Kronbichler, 2011/08/02) +
    2. Removed: deal.II no longer supports Trilinos versions prior to 10.0.
      (Wolfgang Bangerth, 2011/06/29) @@ -50,7 +63,7 @@ corresponding elements of namespace std::placeholders into the global namespace if your compiler supported this part of the C++ 1x standard, or otherwise using the BOOST counterparts which are already in the global namespace. However, this leads to a conflict if one has a C++ 1x enabled compiler (e.g. GCC 4.6) -and #includes certain BOOST headers, since the importation of symbols +and includes certain BOOST headers, since the importation of symbols into the global namespace now leads to ambiguous names. The only solution to the problem is to not import names into the global namespace, but rather import the names from either BOOST or namespace std into the deal.II namespace @@ -65,7 +78,18 @@ changed to use std_cxx1x::_1, std_cxx1x::_2, etc from now on.

      General

        -
      1. Fixed: The function VectorTools::create_right_hand_side now also works for objects of type hp::DoFHandler with different finite elements. +
      2. 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 Tensor<1,dim,std::complex >. deal.II +uses a default template argument double that ensures that all +code using e.g. Tensor<1,dim> remains valid. +
        +(Martin Kronbichler, 2011/08/02) + + +
      3. Fixed: The function VectorTools::create_right_hand_side now also works +for objects of type hp::DoFHandler with different finite elements.
        (Daniel Gerecht, 2011/07/20) @@ -383,12 +407,10 @@ more than rhs.size() CG iterations if the number of degrees of free is very small.
        (Jichao Yin, WB, 2011/04/06) -
      4. New: There is now a new function ConditionalOStream::get_stream().
        (WB, 2011/03/09) -
      5. Fixed: FESystem::get_unit_face_support_points would refuse to return anything if one of the base elements did not have support points. This @@ -397,7 +419,6 @@ base element has no support points and also has degrees of freedom on the face.
        (WB, 2011/03/07) -
      6. Fixed: Objects of type FE_Nothing could be generated with multiple vector components by passing an argument to the constructor. Yet, the FE_Nothing::get_name() function @@ -405,18 +426,15 @@ always just returned the string FE_Nothing@() independently o number of components. This is now fixed.
        (WB, 2011/03/07) -
      7. 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.
        (Timo Heister, 2011/02/23) -
      8. 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.
        (Wolfgang Bangerth, 2011/02/18) -
      9. New: There is now a function DataOutBase::write_visit_record that does the equivalent for VisIt that DataOutBase::write_pvtu_record does for ParaView: @@ -424,25 +442,21 @@ generate a file that contains a list of all other VTK or VTU files of which the current parallel simulation consists.
        (Wolfgang Bangerth, 2011/02/16) -
      10. New: There is now a function TrilinosWrappers::VectorBase::minimal_value.
        (Wolfgang Bangerth, 2011/02/16) -
      11. Fixed: TableBase::operator= could not be compiled if the type of the elements of the table was bool. This is now fixed.
        (Wolfgang Bangerth, 2011/02/16) -
      12. 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.
        (Wolfgang Bangerth, 2011/02/12) -
      13. Fixed: There are systems where the libz library is installed but the zlib.h header file is not available. Since the latter @@ -450,7 +464,6 @@ condition was not tested, this would result in compiler errors. This is now fixed.
        (Wolfgang Bangerth, 2011/02/09) -
      14. Fixed: Prolongation and restriction matrices were not computed at all for elements of type FE_DGQ if dim@. Consequently, @@ -459,7 +472,6 @@ the DoFCellAccess::set_dof_values_by_interpolation function did not work either and simply returned zero results. This is now fixed.
        (M. Sebastian Pauletti, Wolfgang Bangerth, 2011/02/09) -
      15. Fixed: When refining a mesh with codimension one, edges were refined using the same manifold description as adjacent cells, but this ignored that a @@ -467,63 +479,58 @@ boundary indicator might have been purposefully set for edges that are truly at the boundary of the mesh. For such edges, the boundary indicator is now honored.
        (M. Sebastian Pauletti, Wolfgang Bangerth, 2011/02/09) -
      16. Fixed: The functions VectorTools::compute_mean_value and VectorTools::integrate_difference now also work for distributed triangulations of type parallel::distributed::Triangulation.
        (Wolfgang Bangerth, 2011/02/07) -
      17. Changed: If the libz 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).
        (Wolfgang Bangerth, 2011/01/28) -
      18. New: Trilinos and PETSc vectors now have a function has_ghost_elements().
        (Timo Heister, 2011/01/26) -
      19. 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.
        (Wolfgang Bangerth, 2011/01/21) -
      20. Fixed: When calling DataOut::build_patches with a mapping, requesting more than one subdivision, and when dim@, then some cells were not properly mapped. This is now fixed.
        (Wolfgang Bangerth, 2011/01/18) -
      21. 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().
        (Timo Heister, 2011/01/17) -
      22. 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.
        +matrix columns before applying new boundary values. +
        (Martin Stoll, Martin Kronbichler, 2011/01/14) -
      23. Extended: base/tensor.h has an additional collection of contractions between three tensors (ie. contract3). This can be useful for writing matrix/vector assembly in a more compact -form than before.
        +form than before. +
        (Toby D. Young, 2011/01/12) +
      diff --git a/deal.II/include/deal.II/base/function.h b/deal.II/include/deal.II/base/function.h index bd2096a9c7..c9daa892be 100644 --- a/deal.II/include/deal.II/base/function.h +++ b/deal.II/include/deal.II/base/function.h @@ -17,12 +17,11 @@ #include #include #include +#include +#include #include DEAL_II_NAMESPACE_OPEN -template class Point; -template class Tensor; -template class Tensor<1,dim>; template class Vector; /** @@ -109,7 +108,7 @@ class Function : public FunctionTime, * template programming. */ static const unsigned int dimension = dim; - + /** * Number of vector components. */ @@ -126,7 +125,7 @@ class Function : public FunctionTime, */ Function (const unsigned int n_components = 1, const double initial_time = 0.0); - + /** * Virtual destructor; absolutely * necessary in this case. @@ -182,7 +181,7 @@ class Function : public FunctionTime, * to. */ Function & operator= (const Function &f); - + /** * Return the value of the * function at the given @@ -208,7 +207,7 @@ class Function : public FunctionTime, */ virtual void vector_value (const Point &p, Vector &values) const; - + /** * Set values to the point * values of the specified @@ -245,7 +244,7 @@ class Function : public FunctionTime, */ virtual void vector_value_list (const std::vector > &points, std::vector > &values) const; - + /** * For each component of the * function, fill a vector of @@ -277,7 +276,7 @@ class Function : public FunctionTime, */ virtual void vector_gradient (const Point &p, std::vector > &gradients) const; - + /** * Set gradients to the * gradients of the specified @@ -290,7 +289,7 @@ class Function : public FunctionTime, virtual void gradient_list (const std::vector > &points, std::vector > &gradients, const unsigned int component = 0) const; - + /** * For each component of the * function, fill a vector of @@ -306,12 +305,12 @@ class Function : public FunctionTime, */ virtual void vector_gradients (const std::vector > &points, std::vector > > &gradients) const; - + /** * Set gradients to the gradients of * the function at the points, * for all components. - * It is assumed that gradients + * It is assumed that gradients * already has the right size, i.e. * the same size as the points array. * @@ -338,7 +337,7 @@ class Function : public FunctionTime, */ virtual void vector_laplacian (const Point &p, Vector &values) const; - + /** * Compute the Laplacian of one * component at a set of points. @@ -396,13 +395,13 @@ class ZeroFunction : public Function * components is preset to one. */ ZeroFunction (const unsigned int n_components = 1); - + /** * Virtual destructor; absolutely * necessary in this case. */ virtual ~ZeroFunction (); - + /** * Return the value of the function * at the given point for one @@ -443,7 +442,7 @@ class ZeroFunction : public Function */ virtual void vector_value_list (const std::vector > &points, std::vector > &values) const; - + /** * Return the gradient of the function * at the given point, for the @@ -460,7 +459,7 @@ class ZeroFunction : public Function */ virtual void vector_gradient (const Point &p, std::vector > &gradients) const; - + /** * Set gradients to the gradients of * the function at the points, @@ -472,12 +471,12 @@ class ZeroFunction : public Function virtual void gradient_list (const std::vector > &points, std::vector > &gradients, const unsigned int component = 0) const; - + /** * Set gradients to the gradients of * the function at the points, * for all components. - * It is assumed that gradients + * It is assumed that gradients * already has the right size, i.e. * the same size as the points array. * @@ -527,7 +526,7 @@ class ConstantFunction : public ZeroFunction */ ConstantFunction (const double value, const unsigned int n_components = 1); - + /** * Virtual destructor; absolutely * necessary in this case. @@ -656,8 +655,8 @@ class ComponentSelectFunction : public ConstantFunction * number of vector components. */ ComponentSelectFunction (const std::pair &selected, - const unsigned int n_components); - + const unsigned int n_components); + /** * Return the value of the function * at the given point for all diff --git a/deal.II/include/deal.II/base/function_parser.h b/deal.II/include/deal.II/base/function_parser.h index 745f9b347f..40317b10f1 100644 --- a/deal.II/include/deal.II/base/function_parser.h +++ b/deal.II/include/deal.II/base/function_parser.h @@ -20,6 +20,8 @@ #include #include #include +#include +#include #include #include @@ -31,9 +33,6 @@ namespace fparser DEAL_II_NAMESPACE_OPEN -template class Point; -template class Tensor; -template class Tensor<1,dim>; template class Vector; @@ -59,23 +58,23 @@ template class Vector; // Define some constants that will be used by the function parser std::map 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 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 units in expressions. @@ -135,12 +134,12 @@ template class Vector; * For more information on this feature, please see * contrib/functionparser/fparser.txt - - * + + * * See http://warp.povusers.org/FunctionParser/ for an - * explanation on how the underlying library works. + * explanation on how the underlying library works. * - + From the fparser.txt file: @verbatim @@ -182,9 +181,9 @@ template class Vector; example the following expression is valid: x*-y Note that the '=' comparison can be inaccurate due to floating point precision problems (eg. "sqrt(100)=10" probably returns 0, not 1). - + The class supports these functions: - + abs(A) : Absolute value of A. If A is negative, returns -A otherwise returns A. acos(A) : Arc-cosine of A. Returns the angle, measured in radians, @@ -241,9 +240,9 @@ template class Vector; "x-1" "-sin(sqrt(x^2+y^2))" "sqrt(XCoord*XCoord + YCoord*YCoord)" - + An example of a recursive function is the factorial function: - + "if(n>1, n*eval(n-1), 1)" Note that a recursive call has some overhead, which makes it a bit slower @@ -263,53 +262,53 @@ template class Vector; * * An example of time dependent scalar function is the following: @verbatim - + // Empty constants object std::map 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 constants; - + // Variables that will be used inside the expressions std::string variables = "x,y"; - + // Define the expression of the vector valued function. std::string expression = "cos(2*pi*x)*y^2; sin(2*pi*x)*exp(y)"; - + // Generate an empty vector valued function FunctionParser<2> function(2); - + // And populate it with the newly created objects. function.initialize(variables, expression, - constants); - + constants); + @endverbatim - * + * * * @ingroup functions * @author Luca Heltai, 2005 @@ -333,8 +332,8 @@ class FunctionParser : public Function * exception is thrown. */ FunctionParser (const unsigned int n_components = 1, - const double initial_time = 0.0); - + const double initial_time = 0.0); + /** * Destructor. Explicitly delete * the FunctionParser objects @@ -342,14 +341,14 @@ class FunctionParser : public Function * component of the function). */ ~FunctionParser(); - + /** * Type for the constant * map. Used by the initialize() * method. */ typedef std::map ConstMap; - + /** * Iterator for the constants * map. Used by the initialize() @@ -380,7 +379,7 @@ class FunctionParser : public Function * is different from dim (if this * function is not time-dependent) or * from dim+1 (if it is time-dependent). - * + * * expressions: a list of strings * containing the expressions that will * be byte compiled by the internal @@ -404,7 +403,7 @@ class FunctionParser : public Function * name is not valid (eg: * constants["sin"] = 1.5;) an * exception is thrown. - * + * * time_dependent. If this is a * time dependent function, then the last * variable declared in vars is @@ -448,7 +447,7 @@ class FunctionParser : public Function const ConstMap &units, const bool time_dependent = false, const bool use_degrees = false); - + /** * Initialize the function. Same as * above, but accepts a string rather @@ -521,24 +520,24 @@ class FunctionParser : public Function /** @addtogroup Exceptions * @{ */ DeclException2 (ExcParseError, - int, char*, + int, char*, << "Parsing Error at Column " << arg1 << ". The parser said: " << arg2); - + DeclException2 (ExcInvalidExpressionSize, - int, int, + int, int, << "The number of components (" << arg1 - << ") is not equal to the number of expressions (" + << ") is not equal to the number of expressions (" << arg2 << ")."); - + //@} - private: + private: /** * A pointer to the actual * function parsers. */ fparser::FunctionParser * fp; - + /** * State of usability. This * variable is checked every time @@ -547,7 +546,7 @@ class FunctionParser : public Function * in the initialize() methods. */ bool initialized; - + /** * Number of variables. If this * is also a function of time, @@ -584,10 +583,10 @@ FunctionParser::default_variable_names () return ""; } - + DEAL_II_NAMESPACE_CLOSE #endif - + diff --git a/deal.II/include/deal.II/base/point.h b/deal.II/include/deal.II/base/point.h index 600b18275c..72cfc9e22e 100644 --- a/deal.II/include/deal.II/base/point.h +++ b/deal.II/include/deal.II/base/point.h @@ -47,8 +47,8 @@ DEAL_II_NAMESPACE_OPEN * @ingroup geomprimitives * @author Wolfgang Bangerth, 1997 */ -template -class Point : public Tensor<1,dim> +template +class Point : public Tensor<1,dim,Number> { public: /** @@ -56,7 +56,7 @@ class Point : public Tensor<1,dim> * an origin. */ Point (); - + /** * Constructor. Initialize all * entries to zero if @@ -69,8 +69,8 @@ class Point : public Tensor<1,dim> * additional data is inside a point, * this is ok. */ - Point (const Tensor<1,dim> &); - + Point (const Tensor<1,dim,Number> &); + /** * Constructor for one dimensional * points. This function is only @@ -78,7 +78,7 @@ class Point : public Tensor<1,dim> * the usage is considered unsafe for * points with dim!=1. */ - explicit Point (const double x); + explicit Point (const Number x); /** * Constructor for two dimensional @@ -87,8 +87,8 @@ class Point : public Tensor<1,dim> * the usage is considered unsafe for * points with dim!=2. */ - Point (const double x, const double y); - + Point (const Number x, const Number y); + /** * Constructor for three dimensional * points. This function is only @@ -96,38 +96,38 @@ class Point : public Tensor<1,dim> * the usage is considered unsafe for * points with dim!=3. */ - 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 i. */ - static Point unit_vector(const unsigned int i); - + static Point unit_vector(const unsigned int i); + /** * Read access to the indexth * coordinate. */ - double operator () (const unsigned int index) const; + Number operator () (const unsigned int index) const; /** * Read and write access to the * indexth 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 operator += instead * since this does not need to copy a * point at least once. */ - Point operator + (const Tensor<1,dim>&) const; + Point operator + (const Tensor<1,dim,Number>&) const; /** * Subtract two point vectors. If @@ -135,13 +135,13 @@ class Point : public Tensor<1,dim> * instead since this does not need to * copy a point at least once. */ - Point operator - (const Tensor<1,dim>&) const; + Point operator - (const Tensor<1,dim,Number>&) const; /** * The opposite vector. */ - Point operator - () const; - + Point operator - () const; + /** * Multiply by a factor. If possible, * use operator *= instead @@ -151,13 +151,13 @@ class Point : public Tensor<1,dim> * There is a commutative complement to this * function also */ - Point operator * (const double) const; + Point operator * (const Number) const; /** * Returns the scalar product of two * vectors. */ - double operator * (const Tensor<1,dim> &) const; + Number operator * (const Tensor<1,dim,Number> &) const; /** * Divide by a factor. If possible, use @@ -165,15 +165,15 @@ class Point : public Tensor<1,dim> * this does not need to copy a point at * least once. */ - Point operator / (const double) const; + Point 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 * this point to the point @@ -181,12 +181,12 @@ class Point : public Tensor<1,dim> * of the difference between the vectors * representing the two points. */ - double distance (const Point &p) const; + Number distance (const Point &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 void serialize(Archive & ar, const unsigned int version); }; @@ -195,34 +195,34 @@ class Point : public Tensor<1,dim> #ifndef DOXYGEN -template +template inline -Point::Point () +Point::Point () {} -template +template inline -Point::Point (const bool initialize) - : - Tensor<1,dim>(initialize) +Point::Point (const bool initialize) + : + Tensor<1,dim,Number>(initialize) {} -template +template inline -Point::Point (const Tensor<1,dim> &t) - : - Tensor<1,dim>(t) +Point::Point (const Tensor<1,dim,Number> &t) + : + Tensor<1,dim,Number>(t) {} -template +template inline -Point::Point (const double x) +Point::Point (const Number x) { Assert (dim==1, StandardExceptions::ExcInvalidConstructorCall()); this->values[0] = x; @@ -230,9 +230,9 @@ Point::Point (const double x) -template +template inline -Point::Point (const double x, const double y) +Point::Point (const Number x, const Number y) { Assert (dim==2, StandardExceptions::ExcInvalidConstructorCall()); this->values[0] = x; @@ -241,9 +241,9 @@ Point::Point (const double x, const double y) -template +template inline -Point::Point (const double x, const double y, const double z) +Point::Point (const Number x, const Number y, const Number z) { Assert (dim==3, StandardExceptions::ExcInvalidConstructorCall()); this->values[0] = x; @@ -252,20 +252,21 @@ Point::Point (const double x, const double y, const double z) } -template +template inline -Point -Point::unit_vector(unsigned int i) +Point +Point::unit_vector(unsigned int i) { - Point p; + Point p; p[i] = 1.; return p; } -template +template inline -double Point::operator () (const unsigned int index) const +Number +Point::operator () (const unsigned int index) const { AssertIndexRange((int) index, dim); return this->values[index]; @@ -273,9 +274,10 @@ double Point::operator () (const unsigned int index) const -template +template inline -double & Point::operator () (const unsigned int index) +Number & +Point::operator () (const unsigned int index) { AssertIndexRange((int) index, dim); return this->values[index]; @@ -283,29 +285,32 @@ double & Point::operator () (const unsigned int index) -template +template inline -Point Point::operator + (const Tensor<1,dim> &p) const +Point +Point::operator + (const Tensor<1,dim,Number> &p) const { - return (Point(*this) += p); + return (Point(*this) += p); } -template +template inline -Point Point::operator - (const Tensor<1,dim> &p) const +Point +Point::operator - (const Tensor<1,dim,Number> &p) const { - return (Point(*this) -= p); + return (Point(*this) -= p); } -template +template inline -Point Point::operator - () const +Point +Point::operator - () const { - Point result; + Point result; for (unsigned int i=0; ivalues[i]; return result; @@ -313,67 +318,74 @@ Point Point::operator - () const -template +template inline -Point Point::operator * (const double factor) const +Point +Point::operator * (const Number factor) const { - return (Point(*this) *= factor); + return (Point(*this) *= factor); } -template +template inline -double Point::operator * (const Tensor<1,dim> &p) const +Number +Point::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 +template inline -double Point::square () const +Number +Point::square () const { - double q=0; + Number q=0; for (unsigned int i=0; ivalues[i] * this->values[i]; return q; } -template + +template inline -double Point::distance (const Point &p) const +Number +Point::distance (const Point &p) const { - double sum=0; + Number sum=0; for (unsigned int i=0; ivalues[i]-p(i); sum += diff*diff; } - + return std::sqrt(sum); } -template + +template inline -Point Point::operator / (const double factor) const +Point Point::operator / (const Number factor) const { - return (Point(*this) /= factor); + return (Point(*this) /= factor); } -template + +template template inline -void -Point::serialize(Archive & ar, const unsigned int) +void +Point::serialize(Archive & ar, const unsigned int) { // forward to serialization // function in the base class - ar & static_cast &>(*this); + ar & static_cast &>(*this); } #endif // DOXYGEN @@ -382,27 +394,43 @@ Point::serialize(Archive & ar, const unsigned int) /*------------------------------- Global functions: Point ---------------------------*/ +/** + * Global operator scaling a point vector by a scalar. + * @relates Point + */ +template +inline +Point operator * (const Number factor, + const Point &p) +{ + return p*factor; +} + + + /** * Global operator scaling a point vector by a scalar. * @relates Point */ template inline -Point operator * (const double factor, const Point &p) +Point operator * (const double factor, + const Point &p) { return p*factor; } + /** * Output operator for points. Print the elements consecutively, * with a space in between. * @relates Point */ -template +template inline -std::ostream & operator << (std::ostream &out, - const Point &p) +std::ostream & operator << (std::ostream &out, + const Point &p) { for (unsigned int i=0; i +template inline -std::istream & operator >> (std::istream &in, - Point &p) +std::istream & operator >> (std::istream &in, + Point &p) { for (unsigned int i=0; i> p[i]; @@ -432,13 +460,15 @@ std::istream & operator >> (std::istream &in, #ifndef DOXYGEN -/** +/** * Output operator for points of dimension 1. This is implemented * specialized from the general template in order to avoid a compiler - * warning that the loop is empty. + * warning that the loop is empty. */ +template 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]; diff --git a/deal.II/include/deal.II/base/polynomials_raviart_thomas.h b/deal.II/include/deal.II/base/polynomials_raviart_thomas.h index f6425982b1..5e00ef9cb3 100644 --- a/deal.II/include/deal.II/base/polynomials_raviart_thomas.h +++ b/deal.II/include/deal.II/base/polynomials_raviart_thomas.h @@ -59,7 +59,7 @@ class PolynomialsRaviartThomas * Qk contained. */ PolynomialsRaviartThomas (const unsigned int k); - + /** * Computes the value and the * first and second derivatives @@ -87,13 +87,13 @@ class PolynomialsRaviartThomas void compute (const Point &unit_point, std::vector > &values, std::vector > &grads, - std::vector > &grad_grads) const; - + std::vector > &grad_grads) const; + /** * Returns the number of Raviart-Thomas polynomials. */ unsigned int n () const; - + /** * Returns the degree of the Raviart-Thomas * space, which is one less than @@ -106,7 +106,7 @@ class PolynomialsRaviartThomas * which is RaviartThomas. */ std::string name () const; - + /** * Return the number of * polynomials in the space @@ -117,14 +117,14 @@ class PolynomialsRaviartThomas * classes. */ static unsigned int compute_n_pols(unsigned int degree); - + private: /** * The degree of this object as * given to the constructor. */ const unsigned int my_degree; - + /** * An object representing the * polynomial space for a single @@ -153,6 +153,7 @@ class PolynomialsRaviartThomas }; + template inline unsigned int PolynomialsRaviartThomas::n() const @@ -161,6 +162,7 @@ PolynomialsRaviartThomas::n() const } + template inline unsigned int PolynomialsRaviartThomas::degree() const @@ -169,6 +171,7 @@ PolynomialsRaviartThomas::degree() const } + template inline std::string PolynomialsRaviartThomas::name() const diff --git a/deal.II/include/deal.II/base/symmetric_tensor.h b/deal.II/include/deal.II/base/symmetric_tensor.h index 0ab2a688fc..4029965e0b 100644 --- a/deal.II/include/deal.II/base/symmetric_tensor.h +++ b/deal.II/include/deal.II/base/symmetric_tensor.h @@ -18,18 +18,23 @@ DEAL_II_NAMESPACE_OPEN -template class SymmetricTensor; +template class SymmetricTensor; -template SymmetricTensor<2,dim> unit_symmetric_tensor (); -template SymmetricTensor<4,dim> deviator_tensor (); -template SymmetricTensor<4,dim> identity_tensor (); -template SymmetricTensor<4,dim> invert (const SymmetricTensor<4,dim> &); -template double trace (const SymmetricTensor<2,dim2> &); +template SymmetricTensor<2,dim,Number> +unit_symmetric_tensor (); +template SymmetricTensor<4,dim,Number> +deviator_tensor (); +template SymmetricTensor<4,dim,Number> +identity_tensor (); +template SymmetricTensor<4,dim,Number> +invert (const SymmetricTensor<4,dim,Number> &); +template Number +trace (const SymmetricTensor<2,dim2,Number> &); -template SymmetricTensor<2,dim> -deviator (const SymmetricTensor<2,dim> &); -template double -determinant (const SymmetricTensor<2,dim> &); +template SymmetricTensor<2,dim,Number> +deviator (const SymmetricTensor<2,dim,Number> &); +template Number +determinant (const SymmetricTensor<2,dim,Number> &); namespace internal @@ -116,15 +121,15 @@ namespace internal * tensor of rank * rank1+rank2-4, 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 + template struct double_contraction_result { - typedef ::dealii::SymmetricTensor type; + typedef ::dealii::SymmetricTensor type; }; @@ -137,15 +142,15 @@ namespace internal * tensor of rank * rank1+rank2-4, 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 - struct double_contraction_result<2,2,dim> + template + struct double_contraction_result<2,2,dim,Number> { - typedef double type; + typedef Number type; }; @@ -169,15 +174,15 @@ namespace internal * about the storage format in your * application programs. */ - template + template struct StorageType; /** * Specialization of StorageType for * rank-2 tensors. */ - template - struct StorageType<2,dim> + template + struct StorageType<2,dim,Number> { /** * Number of independent components of a @@ -191,7 +196,7 @@ namespace internal * Declare the type in which we actually * store the data. */ - typedef Tensor<1,n_independent_components> base_tensor_type; + typedef Tensor<1,n_independent_components,Number> base_tensor_type; }; @@ -200,8 +205,8 @@ namespace internal * Specialization of StorageType for * rank-4 tensors. */ - template - struct StorageType<4,dim> + template + struct StorageType<4,dim,Number> { /** * Number of independent components @@ -219,7 +224,7 @@ namespace internal */ static const unsigned int n_independent_components = (n_rank2_components * - StorageType<2,dim>::n_independent_components); + StorageType<2,dim,Number>::n_independent_components); /** * Declare the type in which we @@ -230,7 +235,7 @@ namespace internal * matrix if we represent the rank-2 * tensors as vectors. */ - typedef Tensor<2,n_rank2_components> base_tensor_type; + typedef Tensor<2,n_rank2_components,Number> base_tensor_type; }; @@ -241,7 +246,7 @@ namespace internal * switching on whether the tensor * should be constant or not. */ - template + template struct AccessorTypes; /** @@ -252,12 +257,12 @@ namespace internal * * Specialization for constant tensors. */ - template - struct AccessorTypes + template + struct AccessorTypes { - typedef const ::dealii::SymmetricTensor tensor_type; + typedef const ::dealii::SymmetricTensor tensor_type; - typedef double reference; + typedef Number reference; }; /** @@ -269,12 +274,12 @@ namespace internal * Specialization for non-constant * tensors. */ - template - struct AccessorTypes + template + struct AccessorTypes { - typedef ::dealii::SymmetricTensor tensor_type; + typedef ::dealii::SymmetricTensor tensor_type; - typedef double &reference; + typedef Number &reference; }; @@ -313,7 +318,7 @@ namespace internal * * @author Wolfgang Bangerth, 2002, 2005 */ - template + template class Accessor { public: @@ -321,8 +326,8 @@ namespace internal * Import two typedefs from the * switch class above. */ - typedef typename AccessorTypes::reference reference; - typedef typename AccessorTypes::tensor_type tensor_type; + typedef typename AccessorTypes::reference reference; + typedef typename AccessorTypes::tensor_type tensor_type; private: /** @@ -387,7 +392,7 @@ namespace internal /** * Index operator. */ - Accessor operator [] (const unsigned int i); + Accessor operator [] (const unsigned int i); private: /** @@ -402,16 +407,16 @@ namespace internal // work around bugs in some // compilers #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG - template friend class SymmetricTensor; - template + template friend class SymmetricTensor; + template friend class Accessor; # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG - friend class ::dealii::SymmetricTensor; - friend class Accessor; + friend class ::dealii::SymmetricTensor; + friend class Accessor; # endif #else - friend class SymmetricTensor; - friend class Accessor; + friend class SymmetricTensor; + friend class Accessor; #endif }; @@ -427,16 +432,16 @@ namespace internal * * @author Wolfgang Bangerth, 2002, 2005 */ - template - class Accessor + template + class Accessor { public: /** * Import two typedefs from the * switch class above. */ - typedef typename AccessorTypes::reference reference; - typedef typename AccessorTypes::tensor_type tensor_type; + typedef typename AccessorTypes::reference reference; + typedef typename AccessorTypes::tensor_type tensor_type; private: /** @@ -521,16 +526,16 @@ namespace internal // work around bugs in some // compilers #ifndef DEAL_II_NAMESP_TEMPL_FRIEND_BUG - template friend class SymmetricTensor; - template + template friend class SymmetricTensor; + template friend class SymmetricTensorAccessors::Accessor; # ifndef DEAL_II_TEMPL_SPEC_FRIEND_BUG - friend class ::dealii::SymmetricTensor; - friend class SymmetricTensorAccessors::Accessor; + friend class ::dealii::SymmetricTensor; + friend class SymmetricTensorAccessors::Accessor; # endif #else - friend class SymmetricTensor; - friend class Accessor; + friend class SymmetricTensor; + friend class Accessor; #endif }; } @@ -604,7 +609,7 @@ namespace internal * @ingroup geomprimitives * @author Wolfgang Bangerth, 2005 */ -template +template class SymmetricTensor { public: @@ -636,7 +641,7 @@ class SymmetricTensor * tensors of rank 2. */ static const unsigned int n_independent_components - = internal::SymmetricTensorAccessors::StorageType:: + = internal::SymmetricTensorAccessors::StorageType:: n_independent_components; /** @@ -663,7 +668,7 @@ class SymmetricTensor * practice to check before calling * symmetrize. */ - SymmetricTensor (const Tensor<2,dim> &t); + SymmetricTensor (const Tensor<2,dim,Number> &t); /** * A constructor that creates a @@ -692,7 +697,7 @@ class SymmetricTensor * namespace is to work around * bugs in some older compilers. */ - SymmetricTensor (const double (&array) [internal::SymmetricTensorAccessors::StorageType::n_independent_components]); + SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType::n_independent_components]); /** * Assignment operator. @@ -711,7 +716,7 @@ class SymmetricTensor * all elements of the tensor to * zero. */ - SymmetricTensor & operator = (const double d); + SymmetricTensor & operator = (const Number d); /** * Convert the present symmetric tensor @@ -719,7 +724,7 @@ class SymmetricTensor * elements, but using the different * storage scheme of full tensors. */ - operator Tensor () const; + operator Tensor () const; /** * Test for equality of two tensors. @@ -746,13 +751,13 @@ class SymmetricTensor * i.e. multiply all components by * factor. */ - SymmetricTensor & operator *= (const double factor); + SymmetricTensor & operator *= (const Number factor); /** * Scale the vector by * 1/factor. */ - SymmetricTensor & operator /= (const double factor); + SymmetricTensor & operator /= (const Number factor); /** * Add two tensors. If possible, you @@ -827,8 +832,8 @@ class SymmetricTensor * into the first argument to the * function. */ - typename internal::SymmetricTensorAccessors::double_contraction_result::type - operator * (const SymmetricTensor<2,dim> &s) const; + typename internal::SymmetricTensorAccessors::double_contraction_result::type + operator * (const SymmetricTensor<2,dim,Number> &s) const; /** * Contraction over two indices @@ -836,14 +841,14 @@ class SymmetricTensor * rank-4 symmetric tensor given * as argument. */ - typename internal::SymmetricTensorAccessors::double_contraction_result::type - operator * (const SymmetricTensor<4,dim> &s) const; + typename internal::SymmetricTensorAccessors::double_contraction_result::type + operator * (const SymmetricTensor<4,dim,Number> &s) const; /** * Return a read-write reference * to the indicated element. */ - double & operator() (const TableIndices &indices); + Number & operator() (const TableIndices &indices); /** * Return the value of the @@ -858,14 +863,14 @@ class SymmetricTensor * don't know here whether * copying is expensive or not. */ - double operator() (const TableIndices &indices) const; + Number operator() (const TableIndices &indices) const; /** * Access the elements of a row of this * symmetric tensor. This function is * called for constant tensors. */ - internal::SymmetricTensorAccessors::Accessor + internal::SymmetricTensorAccessors::Accessor operator [] (const unsigned int row) const; /** @@ -873,7 +878,7 @@ class SymmetricTensor * symmetric tensor. This function is * called for non-constant tensors. */ - internal::SymmetricTensorAccessors::Accessor + internal::SymmetricTensorAccessors::Accessor operator [] (const unsigned int row); /** @@ -881,7 +886,7 @@ class SymmetricTensor * specify the entire set of * indices. */ - double + Number operator [] (const TableIndices &indices) const; /** @@ -889,7 +894,7 @@ class SymmetricTensor * specify the entire set of * indices. */ - double & + Number & operator [] (const TableIndices &indices); /** @@ -907,7 +912,7 @@ class SymmetricTensor * just one of them, although they are * equal for symmetric tensors). */ - double norm () const; + Number norm () const; /** * Tensors can be unrolled by @@ -980,7 +985,7 @@ class SymmetricTensor * properties of the base tensor. */ typedef - internal::SymmetricTensorAccessors::StorageType + internal::SymmetricTensorAccessors::StorageType base_tensor_descriptor; /** @@ -998,32 +1003,32 @@ class SymmetricTensor /** * Make all other symmetric tensors friends. */ - template friend class SymmetricTensor; + template friend class SymmetricTensor; /** * Make a few more functions friends. */ - template - friend double trace (const SymmetricTensor<2,dim2> &d); + template + friend Number2 trace (const SymmetricTensor<2,dim2,Number2> &d); - template - friend double determinant (const SymmetricTensor<2,dim2> &t); + template + friend Number2 determinant (const SymmetricTensor<2,dim2,Number2> &t); - template - friend SymmetricTensor<2,dim2> - deviator (const SymmetricTensor<2,dim2> &t); + template + friend SymmetricTensor<2,dim2,Number2> + deviator (const SymmetricTensor<2,dim2,Number2> &t); - template - friend SymmetricTensor<2,dim2> unit_symmetric_tensor (); + template + friend SymmetricTensor<2,dim2,Number2> unit_symmetric_tensor (); - template - friend SymmetricTensor<4,dim2> deviator_tensor (); + template + friend SymmetricTensor<4,dim2,Number2> deviator_tensor (); - template - friend SymmetricTensor<4,dim2> identity_tensor (); + template + friend SymmetricTensor<4,dim2,Number2> identity_tensor (); - template - friend SymmetricTensor<4,dim2> invert (const SymmetricTensor<4,dim2> &); + template + friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &); }; @@ -1036,8 +1041,8 @@ namespace internal { namespace SymmetricTensorAccessors { - template - Accessor:: + template + Accessor:: Accessor (tensor_type &tensor, const TableIndices &previous_indices) : @@ -1047,18 +1052,18 @@ namespace internal - template - Accessor - Accessor::operator[] (const unsigned int i) + template + Accessor + Accessor::operator[] (const unsigned int i) { - return Accessor (tensor, - merge (previous_indices, i, rank-P)); + return Accessor (tensor, + merge (previous_indices, i, rank-P)); } - template - Accessor:: + template + Accessor:: Accessor (tensor_type &tensor, const TableIndices &previous_indices) : @@ -1068,9 +1073,9 @@ namespace internal - template - typename Accessor::reference - Accessor::operator[] (const unsigned int i) + template + typename Accessor::reference + Accessor::operator[] (const unsigned int i) { return tensor(merge (previous_indices, i, rank-1)); } @@ -1081,55 +1086,61 @@ namespace internal -template +template inline -SymmetricTensor::SymmetricTensor () +SymmetricTensor::SymmetricTensor () {} -template <> +template inline -SymmetricTensor<2,2>::SymmetricTensor (const Tensor<2,2> &t) +SymmetricTensor::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 +template inline -SymmetricTensor::SymmetricTensor (const double (&array) [internal::SymmetricTensorAccessors::StorageType::n_independent_components]) +SymmetricTensor::SymmetricTensor (const Number (&array) [internal::SymmetricTensorAccessors::StorageType::n_independent_components]) : data (array) {} -template + +template inline -SymmetricTensor & -SymmetricTensor::operator = (const SymmetricTensor &t) +SymmetricTensor & +SymmetricTensor::operator = (const SymmetricTensor &t) { data = t.data; return *this; @@ -1137,10 +1148,10 @@ SymmetricTensor::operator = (const SymmetricTensor &t) -template +template inline -SymmetricTensor & -SymmetricTensor::operator = (const double d) +SymmetricTensor & +SymmetricTensor::operator = (const Number d) { Assert (d==0, ExcMessage ("Only assignment with zero is allowed")); @@ -1151,66 +1162,81 @@ SymmetricTensor::operator = (const double d) -template <> -inline -SymmetricTensor<2,1>:: -operator Tensor<2,1> () const + // helper function to convert symmetric tensor + // to generic tensor +namespace internal { - const double t[1][1] = {{data[0]}}; - return Tensor<2,1>(t); -} - + template + 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 + 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 + 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 inline -SymmetricTensor<2,3>:: -operator Tensor<2,3> () const +SymmetricTensor:: +operator Tensor () 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 +template inline bool -SymmetricTensor::operator == (const SymmetricTensor &t) const +SymmetricTensor::operator == +(const SymmetricTensor &t) const { return data == t.data; } -template +template inline bool -SymmetricTensor::operator != (const SymmetricTensor &t) const +SymmetricTensor::operator != +(const SymmetricTensor &t) const { return data != t.data; } -template +template inline -SymmetricTensor & -SymmetricTensor::operator += (const SymmetricTensor &t) +SymmetricTensor & +SymmetricTensor::operator += +(const SymmetricTensor &t) { data += t.data; return *this; @@ -1218,10 +1244,11 @@ SymmetricTensor::operator += (const SymmetricTensor &t) -template +template inline -SymmetricTensor & -SymmetricTensor::operator -= (const SymmetricTensor &t) +SymmetricTensor & +SymmetricTensor::operator -= +(const SymmetricTensor &t) { data -= t.data; return *this; @@ -1229,10 +1256,10 @@ SymmetricTensor::operator -= (const SymmetricTensor &t) -template +template inline -SymmetricTensor & -SymmetricTensor::operator *= (const double d) +SymmetricTensor & +SymmetricTensor::operator *= (const Number d) { data *= d; return *this; @@ -1240,10 +1267,10 @@ SymmetricTensor::operator *= (const double d) -template +template inline -SymmetricTensor & -SymmetricTensor::operator /= (const double d) +SymmetricTensor & +SymmetricTensor::operator /= (const Number d) { data /= d; return *this; @@ -1251,10 +1278,10 @@ SymmetricTensor::operator /= (const double d) -template +template inline -SymmetricTensor -SymmetricTensor::operator + (const SymmetricTensor &t) const +SymmetricTensor +SymmetricTensor::operator + (const SymmetricTensor &t) const { SymmetricTensor tmp = *this; tmp.data += t.data; @@ -1263,10 +1290,10 @@ SymmetricTensor::operator + (const SymmetricTensor &t) const -template +template inline -SymmetricTensor -SymmetricTensor::operator - (const SymmetricTensor &t) const +SymmetricTensor +SymmetricTensor::operator - (const SymmetricTensor &t) const { SymmetricTensor tmp = *this; tmp.data -= t.data; @@ -1275,10 +1302,10 @@ SymmetricTensor::operator - (const SymmetricTensor &t) const -template +template inline -SymmetricTensor -SymmetricTensor::operator - () const +SymmetricTensor +SymmetricTensor::operator - () const { SymmetricTensor tmp = *this; tmp.data = -tmp.data; @@ -1287,417 +1314,335 @@ SymmetricTensor::operator - () const -template +template inline void -SymmetricTensor::clear () +SymmetricTensor::clear () { data.clear (); } -template +template inline std::size_t -SymmetricTensor::memory_consumption () +SymmetricTensor::memory_consumption () { return - internal::SymmetricTensorAccessors::StorageType::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::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 + 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 + 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 + 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 + 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 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::type +SymmetricTensor::operator * (const SymmetricTensor<2,dim,Number> &s) const { - const unsigned int dim = 2; - SymmetricTensor<4,dim> tmp; - for (unsigned int i=0; i (data, s.data); } -template <> +template 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::type +SymmetricTensor::operator * (const SymmetricTensor<4,dim,Number> &s) const { - const unsigned int dim = 3; - SymmetricTensor<4,dim> tmp; - for (unsigned int i=0; i::type tmp; + tmp.data = internal::perform_double_contraction (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 -inline -double -SymmetricTensor<2,1>::operator () (const TableIndices<2> &indices) const -{ - const unsigned int rank = 2; - for (unsigned int r=0; r -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 + 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 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 + 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 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 -inline -double -SymmetricTensor<4,1>::operator () (const TableIndices<4> &indices) const -{ - const unsigned int rank = 4; - for (unsigned int r=0; r 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 + inline + Number & + symmetric_tensor_access (const TableIndices<4> &indices, + typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data) + { + switch (dim) + { + case 1: + return data[0][0]; + case 2: // each entry of the tensor can be // thought of as an entry in a // matrix that maps the rolled-out @@ -1707,35 +1652,26 @@ SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices) // tensors. determine which // position the present entry is // stored in - unsigned int base_index[2] ; - if ((indices[0] == 0) && (indices[1] == 0)) - base_index[0] = 0; - else if ((indices[0] == 1) && (indices[1] == 1)) - base_index[0] = 1; - else - base_index[0] = 2; - - if ((indices[2] == 0) && (indices[3] == 0)) - base_index[1] = 0; - else if ((indices[2] == 1) && (indices[3] == 1)) - base_index[1] = 1; - else - base_index[1] = 2; - - return data[base_index[0]][base_index[1]]; -} - - - -template <> -inline -double -SymmetricTensor<4,2>::operator () (const TableIndices<4> &indices) const -{ - const unsigned int rank = 4; - for (unsigned int r=0; r::operator () (const TableIndices<4> &indices) const // tensors. determine which // position the present entry is // stored in - unsigned int base_index[2] ; - if ((indices[0] == 0) && (indices[1] == 0)) - base_index[0] = 0; - else if ((indices[0] == 1) && (indices[1] == 1)) - base_index[0] = 1; - else - base_index[0] = 2; - - if ((indices[2] == 0) && (indices[3] == 0)) - base_index[1] = 0; - else if ((indices[2] == 1) && (indices[3] == 1)) - base_index[1] = 1; - else - base_index[1] = 2; - - return data[base_index[0]][base_index[1]]; -} - + { + unsigned int base_index[2] ; + if ((indices[0] == 0) && (indices[1] == 0)) + base_index[0] = 0; + else if ((indices[0] == 1) && (indices[1] == 1)) + base_index[0] = 1; + else if ((indices[0] == 2) && (indices[1] == 2)) + base_index[0] = 2; + else if (((indices[0] == 0) && (indices[1] == 1)) || + ((indices[0] == 1) && (indices[1] == 0))) + base_index[0] = 3; + else if (((indices[0] == 0) && (indices[1] == 2)) || + ((indices[0] == 2) && (indices[1] == 0))) + base_index[0] = 4; + else + { + Assert (((indices[0] == 1) && (indices[1] == 2)) || + ((indices[0] == 2) && (indices[1] == 1)), + ExcInternalError()); + base_index[0] = 5; + } + + if ((indices[2] == 0) && (indices[3] == 0)) + base_index[1] = 0; + else if ((indices[2] == 1) && (indices[3] == 1)) + base_index[1] = 1; + else if ((indices[2] == 2) && (indices[3] == 2)) + base_index[1] = 2; + else if (((indices[2] == 0) && (indices[3] == 1)) || + ((indices[2] == 1) && (indices[3] == 0))) + base_index[1] = 3; + else if (((indices[2] == 0) && (indices[3] == 2)) || + ((indices[2] == 2) && (indices[3] == 0))) + base_index[1] = 4; + else + { + Assert (((indices[2] == 1) && (indices[3] == 2)) || + ((indices[2] == 2) && (indices[3] == 1)), + ExcInternalError()); + base_index[1] = 5; + } + + return data[base_index[0]][base_index[1]]; + } + } + } -template <> -inline -double & -SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices) -{ - const unsigned int rank = 4; - for (unsigned int r=0; r + inline + Number + symmetric_tensor_access (const TableIndices<4> &indices, + const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data) + { + switch (dim) + { + case 1: + return data[0][0]; + case 2: // each entry of the tensor can be // thought of as an entry in a // matrix that maps the rolled-out @@ -1783,61 +1750,26 @@ SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices) // tensors. determine which // position the present entry is // stored in - unsigned int base_index[2] ; - if ((indices[0] == 0) && (indices[1] == 0)) - base_index[0] = 0; - else if ((indices[0] == 1) && (indices[1] == 1)) - base_index[0] = 1; - else if ((indices[0] == 2) && (indices[1] == 2)) - base_index[0] = 2; - else if (((indices[0] == 0) && (indices[1] == 1)) || - ((indices[0] == 1) && (indices[1] == 0))) - base_index[0] = 3; - else if (((indices[0] == 0) && (indices[1] == 2)) || - ((indices[0] == 2) && (indices[1] == 0))) - base_index[0] = 4; - else - { - Assert (((indices[0] == 1) && (indices[1] == 2)) || - ((indices[0] == 2) && (indices[1] == 1)), - ExcInternalError()); - base_index[0] = 5; - } - - if ((indices[2] == 0) && (indices[3] == 0)) - base_index[1] = 0; - else if ((indices[2] == 1) && (indices[3] == 1)) - base_index[1] = 1; - else if ((indices[2] == 2) && (indices[3] == 2)) - base_index[1] = 2; - else if (((indices[2] == 0) && (indices[3] == 1)) || - ((indices[2] == 1) && (indices[3] == 0))) - base_index[1] = 3; - else if (((indices[2] == 0) && (indices[3] == 2)) || - ((indices[2] == 2) && (indices[3] == 0))) - base_index[1] = 4; - else - { - Assert (((indices[2] == 1) && (indices[3] == 2)) || - ((indices[2] == 2) && (indices[3] == 1)), - ExcInternalError()); - base_index[1] = 5; - } - - return data[base_index[0]][base_index[1]]; -} - - - -template <> -inline -double -SymmetricTensor<4,3>::operator () (const TableIndices<4> &indices) const -{ - const unsigned int rank = 4; - for (unsigned int r=0; r::operator () (const TableIndices<4> &indices) const // tensors. determine which // position the present entry is // stored in - unsigned int base_index[2] ; - if ((indices[0] == 0) && (indices[1] == 0)) - base_index[0] = 0; - else if ((indices[0] == 1) && (indices[1] == 1)) - base_index[0] = 1; - else if ((indices[0] == 2) && (indices[1] == 2)) - base_index[0] = 2; - else if (((indices[0] == 0) && (indices[1] == 1)) || - ((indices[0] == 1) && (indices[1] == 0))) - base_index[0] = 3; - else if (((indices[0] == 0) && (indices[1] == 2)) || - ((indices[0] == 2) && (indices[1] == 0))) - base_index[0] = 4; - else - { - Assert (((indices[0] == 1) && (indices[1] == 2)) || - ((indices[0] == 2) && (indices[1] == 1)), - ExcInternalError()); - base_index[0] = 5; - } - - if ((indices[2] == 0) && (indices[3] == 0)) - base_index[1] = 0; - else if ((indices[2] == 1) && (indices[3] == 1)) - base_index[1] = 1; - else if ((indices[2] == 2) && (indices[3] == 2)) - base_index[1] = 2; - else if (((indices[2] == 0) && (indices[3] == 1)) || - ((indices[2] == 1) && (indices[3] == 0))) - base_index[1] = 3; - else if (((indices[2] == 0) && (indices[3] == 2)) || - ((indices[2] == 2) && (indices[3] == 0))) - base_index[1] = 4; - else - { - Assert (((indices[2] == 1) && (indices[3] == 2)) || - ((indices[2] == 2) && (indices[3] == 1)), - ExcInternalError()); - base_index[1] = 5; - } - - return data[base_index[0]][base_index[1]]; -} - - - -template -internal::SymmetricTensorAccessors::Accessor -SymmetricTensor::operator [] (const unsigned int row) const -{ - return - internal::SymmetricTensorAccessors:: - Accessor (*this, TableIndices (row)); -} - - - -template -internal::SymmetricTensorAccessors::Accessor -SymmetricTensor::operator [] (const unsigned int row) -{ - return - internal::SymmetricTensorAccessors:: - Accessor (*this, TableIndices (row)); -} - - - -template -inline -double -SymmetricTensor::operator [] (const TableIndices &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 -inline -double & -SymmetricTensor::operator [] (const TableIndices &indices) -{ - return data[component_to_unrolled_index(indices)]; -} +} // end of namespace internal -template <> +template inline -double -SymmetricTensor<2,1>::norm () const +Number & +SymmetricTensor::operator () (const TableIndices &indices) { - return std::fabs(data[0]); + for (unsigned int r=0; r (indices, data); } -template <> +template inline -double -SymmetricTensor<2,2>::norm () const +Number +SymmetricTensor::operator () + (const TableIndices &indices) const { - return std::sqrt(data[0]*data[0] + data[1]*data[1] + 2*data[2]*data[2]); + for (unsigned int r=0; r (indices, data); } -template <> -inline -double -SymmetricTensor<2,3>::norm () const +template +internal::SymmetricTensorAccessors::Accessor +SymmetricTensor::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 (*this, TableIndices (row)); } -template <> -inline -double -SymmetricTensor<4,1>::norm () const +template +internal::SymmetricTensorAccessors::Accessor +SymmetricTensor::operator [] (const unsigned int row) { - return std::fabs(data[0][0]); + return + internal::SymmetricTensorAccessors:: + Accessor (*this, TableIndices (row)); } -template <> +template inline -double -SymmetricTensor<4,2>::norm () const +Number +SymmetricTensor::operator [] (const TableIndices &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 +template inline -double -SymmetricTensor<4,3>::norm () const +Number & +SymmetricTensor::operator [] (const TableIndices &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 -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 + 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 -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 + 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 +template inline -unsigned int -SymmetricTensor<2,3>::component_to_unrolled_index (const TableIndices<2> &indices) +Number +SymmetricTensor::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 (data); } -template <> +template inline unsigned int -SymmetricTensor<2,4>::component_to_unrolled_index (const TableIndices<2> &indices) +SymmetricTensor::component_to_unrolled_index + (const TableIndices &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 inline -TableIndices<2> -SymmetricTensor<2,2>::unrolled_to_component_indices (const unsigned int i) +TableIndices +SymmetricTensor::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 +template template inline void -SymmetricTensor::serialize(Archive & ar, const unsigned int) +SymmetricTensor::serialize(Archive & ar, const unsigned int) { ar & data; } @@ -2159,47 +2079,8 @@ SymmetricTensor::serialize(Archive & ar, const unsigned int) * Compute the determinant of a tensor or rank 2. The determinant is * also commonly referred to as the third invariant of rank-2 tensors. * - * For the present case of one-dimensional tensors, the determinant - * equals the only element and is therefore equivalent to the trace. - * - * For greater notational simplicity, there is also a third_invariant - * 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 third_invariant - * 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 third_invariant * function that returns the determinant of a tensor. @@ -2207,18 +2088,29 @@ double determinant (const SymmetricTensor<2,2> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template <> +template inline -double determinant (const SymmetricTensor<2,3> &t) +Number determinant (const SymmetricTensor<2,dim,Number> &t) { + switch (dim) + { + case 1: + return t.data[0]; + case 2: + return (t.data[0] * t.data[1] - t.data[2]*t.data[2]); + case 3: // in analogy to general tensors, but // there's something to be simplified for // the present case - return ( t.data[0]*t.data[1]*t.data[2] - -t.data[0]*t.data[5]*t.data[5] - -t.data[1]*t.data[4]*t.data[4] - -t.data[2]*t.data[3]*t.data[3] - +2*t.data[3]*t.data[4]*t.data[5] ); + return ( t.data[0]*t.data[1]*t.data[2] + -t.data[0]*t.data[5]*t.data[5] + -t.data[1]*t.data[4]*t.data[4] + -t.data[2]*t.data[3]*t.data[3] + +2*t.data[3]*t.data[4]*t.data[5] ); + default: + Assert (false, ExcNotImplemented()); + return 0; + } } @@ -2232,9 +2124,9 @@ double determinant (const SymmetricTensor<2,3> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template inline -double third_invariant (const SymmetricTensor<2,dim> &t) +double third_invariant (const SymmetricTensor<2,dim,Number> &t) { return determinant (t); } @@ -2249,10 +2141,10 @@ double third_invariant (const SymmetricTensor<2,dim> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template -double trace (const SymmetricTensor<2,dim> &d) +template +Number trace (const SymmetricTensor<2,dim,Number> &d) { - double t=0; + Number t=0; for (unsigned int i=0; i &d) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template inline -double first_invariant (const SymmetricTensor<2,dim> &t) +Number first_invariant (const SymmetricTensor<2,dim,Number> &t) { return trace (t); } @@ -2283,8 +2175,9 @@ double first_invariant (const SymmetricTensor<2,dim> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005, 2010 */ +template inline -double second_invariant (const SymmetricTensor<2,1> &) +Number second_invariant (const SymmetricTensor<2,1,Number> &) { return 0; } @@ -2298,8 +2191,9 @@ double second_invariant (const SymmetricTensor<2,1> &) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005, 2010 */ +template inline -double second_invariant (const SymmetricTensor<2,2> &t) +Number second_invariant (const SymmetricTensor<2,2,Number> &t) { return t[0][0]*t[1][1] - t[0][1]*t[0][1]; } @@ -2313,8 +2207,9 @@ double second_invariant (const SymmetricTensor<2,2> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005, 2010 */ +template inline -double second_invariant (const SymmetricTensor<2,3> &t) +Number second_invariant (const SymmetricTensor<2,3,Number> &t) { return (t[0][0]*t[1][1] + t[1][1]*t[2][2] + t[2][2]*t[0][0] - t[0][1]*t[0][1] - t[0][2]*t[0][2] - t[1][2]*t[1][2]); @@ -2332,10 +2227,10 @@ double second_invariant (const SymmetricTensor<2,3> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template inline -SymmetricTensor -transpose (const SymmetricTensor &t) +SymmetricTensor +transpose (const SymmetricTensor &t) { return t; } @@ -2351,15 +2246,15 @@ transpose (const SymmetricTensor &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &t) /** - * Return a unit symmetric tensor of rank 2 and dimension 1. - * - * @relates SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -template <> -inline -SymmetricTensor<2,1> -unit_symmetric_tensor<1> () -{ - SymmetricTensor<2,1> tmp; - tmp.data[0] = 1; - return tmp; -} - - - -/** - * Return a unit symmetric tensor of rank 2 and dimension 2. - * - * @relates SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -template <> -inline -SymmetricTensor<2,2> -unit_symmetric_tensor<2> () -{ - SymmetricTensor<2,2> tmp; - tmp.data[0] = tmp.data[1] = 1; - return tmp; -} - - - -/** - * Return a unit symmetric tensor of rank 2 and dimension 3. + * Return a unit symmetric tensor of rank 2. * * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template <> +template 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 () * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template 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::n_rank2_components; + i::n_rank2_components; ++i) tmp.data[i][i] = 0.5; @@ -2485,12 +2358,12 @@ deviator_tensor () * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template 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::n_rank2_components; + i::n_rank2_components; ++i) tmp.data[i][i] = 0.5; @@ -2515,7 +2388,7 @@ identity_tensor () * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are * mappings from and to symmetric rank-2 tensors, they can have an * inverse. This function computes it, if it exists, for the case that the - * dimension equals 1. + * dimension equals either 1 or 2. * * If a tensor is not invertible, then the result is unspecified, but will * likely contain the results of a division by zero or a very small number at @@ -2524,37 +2397,18 @@ identity_tensor () * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template <> +template inline -SymmetricTensor<4,1> -invert (const SymmetricTensor<4,1> &t) +SymmetricTensor<4,dim,Number> +invert (const SymmetricTensor<4,dim,Number> &t) { - SymmetricTensor<4,1> tmp; - tmp.data[0][0] = 1./t.data[0][0]; - return tmp; -} - - - -/** - * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are - * mappings from and to symmetric rank-2 tensors, they can have an - * inverse. This function computes it, if it exists, for the case that the - * dimension equals 2. - * - * If a tensor is not invertible, then the result is unspecified, but will - * likely contain the results of a division by zero or a very small number at - * the very least. - * - * @relates SymmetricTensor - * @author Wolfgang Bangerth, 2005 - */ -template <> -inline -SymmetricTensor<4,2> -invert (const SymmetricTensor<4,2> &t) -{ - SymmetricTensor<4,2> tmp; + SymmetricTensor<4,dim,Number> tmp; + switch (dim) + { + case 1: + tmp.data[0][0] = 1./t.data[0][0]; + break; + case 2: // inverting this tensor is a little more // complicated than necessary, since we @@ -2588,32 +2442,38 @@ invert (const SymmetricTensor<4,2> &t) // column of the inverse by 1/2, // corresponding to the left and right // multiplication with mult^-1 - const double t4 = t.data[0][0]*t.data[1][1], - t6 = t.data[0][0]*t.data[1][2], - t8 = t.data[0][1]*t.data[1][0], - t00 = t.data[0][2]*t.data[1][0], - t01 = t.data[0][1]*t.data[2][0], - t04 = t.data[0][2]*t.data[2][0], - t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]-t8*t.data[2][2]+ - t00*t.data[2][1]+t01*t.data[1][2]-t04*t.data[1][1]); - tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07; - tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07; - tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07; - tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07; - tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07; - tmp.data[1][2] = -(t6-t00)*t07; - tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07; - tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07; - tmp.data[2][2] = (t4-t8)*t07; + { + const Number t4 = t.data[0][0]*t.data[1][1], + t6 = t.data[0][0]*t.data[1][2], + t8 = t.data[0][1]*t.data[1][0], + t00 = t.data[0][2]*t.data[1][0], + t01 = t.data[0][1]*t.data[2][0], + t04 = t.data[0][2]*t.data[2][0], + t07 = 1.0/(t4*t.data[2][2]-t6*t.data[2][1]- + t8*t.data[2][2]+t00*t.data[2][1]+ + t01*t.data[1][2]-t04*t.data[1][1]); + tmp.data[0][0] = (t.data[1][1]*t.data[2][2]-t.data[1][2]*t.data[2][1])*t07; + tmp.data[0][1] = -(t.data[0][1]*t.data[2][2]-t.data[0][2]*t.data[2][1])*t07; + tmp.data[0][2] = -(-t.data[0][1]*t.data[1][2]+t.data[0][2]*t.data[1][1])*t07; + tmp.data[1][0] = -(t.data[1][0]*t.data[2][2]-t.data[1][2]*t.data[2][0])*t07; + tmp.data[1][1] = (t.data[0][0]*t.data[2][2]-t04)*t07; + tmp.data[1][2] = -(t6-t00)*t07; + tmp.data[2][0] = -(-t.data[1][0]*t.data[2][1]+t.data[1][1]*t.data[2][0])*t07; + tmp.data[2][1] = -(t.data[0][0]*t.data[2][1]-t01)*t07; + tmp.data[2][2] = (t4-t8)*t07; // scale last row and column as mentioned // above - tmp.data[2][0] /= 2; - tmp.data[2][1] /= 2; - tmp.data[0][2] /= 2; - tmp.data[1][2] /= 2; - tmp.data[2][2] /= 4; - + tmp.data[2][0] /= 2; + tmp.data[2][1] /= 2; + tmp.data[0][2] /= 2; + tmp.data[1][2] /= 2; + tmp.data[2][2] /= 4; + } + break; + default: + Assert (false, ExcNotImplemented()); + } return tmp; } @@ -2633,9 +2493,9 @@ invert (const SymmetricTensor<4,2> &t) * @author Wolfgang Bangerth, 2005 */ template <> -SymmetricTensor<4,3> -invert (const SymmetricTensor<4,3> &t); -// this function is implemented in the .cc file +SymmetricTensor<4,3,double> +invert (const SymmetricTensor<4,3,double> &t); +// this function is implemented in the .cc file for double data types /** @@ -2644,21 +2504,21 @@ invert (const SymmetricTensor<4,3> &t); * T phi = t1 (t2 : phi) for all symmetric tensors phi. * * For example, the deviator tensor can be computed as - * identity_tensor() - - * 1/d*outer_product(unit_symmetric_tensor(), - * unit_symmetric_tensor()), since the (double) contraction with the + * identity_tensor() - + * 1/d*outer_product(unit_symmetric_tensor(), + * unit_symmetric_tensor()), since the (double) contraction with the * unit tensor yields the trace of a symmetric tensor. * * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &t1, * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline -SymmetricTensor<2,1> -symmetrize (const Tensor<2,1> &t) +SymmetricTensor<2,1,Number> +symmetrize (const Tensor<2,1,Number> &t) { - const double array[1] + const Number array[1] = { t[0][0] }; - return SymmetricTensor<2,1>(array); + return SymmetricTensor<2,1,Number>(array); } @@ -2699,13 +2560,14 @@ symmetrize (const Tensor<2,1> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline -SymmetricTensor<2,2> -symmetrize (const Tensor<2,2> &t) +SymmetricTensor<2,2,Number> +symmetrize (const Tensor<2,2,Number> &t) { - const double array[3] + const Number array[3] = { t[0][0], t[1][1], (t[0][1] + t[1][0])/2 }; - return SymmetricTensor<2,2>(array); + return SymmetricTensor<2,2,Number>(array); } @@ -2718,21 +2580,79 @@ symmetrize (const Tensor<2,2> &t) * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template 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 +inline +SymmetricTensor +operator * (const SymmetricTensor &t, + const Number factor) +{ + SymmetricTensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Multiplication of a symmetric tensor of general rank with a scalar + * from the left. + * + * @relates SymmetricTensor + */ +template +inline +SymmetricTensor +operator * (const Number factor, + const SymmetricTensor &t) +{ + SymmetricTensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Division of a symmetric tensor of general rank by a scalar. + * + * @relates SymmetricTensor + */ +template +inline +SymmetricTensor +operator / (const SymmetricTensor &t, + const Number factor) +{ + SymmetricTensor tt = t; + tt /= factor; + return tt; +} + + + /** - * Multiplication of a symmetric tensor of general rank with a scalar double + * Multiplication of a symmetric tensor of general rank with a scalar * from the right. * * @relates SymmetricTensor @@ -2741,7 +2661,7 @@ template inline SymmetricTensor operator * (const SymmetricTensor &t, - const double factor) + const double factor) { SymmetricTensor tt = t; tt *= factor; @@ -2751,7 +2671,7 @@ operator * (const SymmetricTensor &t, /** - * Multiplication of a symmetric tensor of general rank with a scalar double + * Multiplication of a symmetric tensor of general rank with a scalar * from the left. * * @relates SymmetricTensor @@ -2759,7 +2679,7 @@ operator * (const SymmetricTensor &t, template inline SymmetricTensor -operator * (const double factor, +operator * (const double factor, const SymmetricTensor &t) { SymmetricTensor tt = t; @@ -2770,7 +2690,7 @@ operator * (const double factor, /** - * Division of a symmetric tensor of general rank by a scalar double. + * Division of a symmetric tensor of general rank by a scalar. * * @relates SymmetricTensor */ @@ -2778,7 +2698,7 @@ template inline SymmetricTensor operator / (const SymmetricTensor &t, - const double factor) + const double factor) { SymmetricTensor tt = t; tt /= factor; @@ -2803,11 +2723,12 @@ operator / (const SymmetricTensor &t, * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,1> &tmp, - const SymmetricTensor<4,1> &t, - const SymmetricTensor<2,1> &s) +double_contract (SymmetricTensor<2,1,Number> &tmp, + const SymmetricTensor<4,1,Number> &t, + const SymmetricTensor<2,1,Number> &s) { tmp[0][0] = t[0][0][0][0] * s[0][0]; } @@ -2829,11 +2750,12 @@ double_contract (SymmetricTensor<2,1> &tmp, * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,1> &tmp, - const SymmetricTensor<2,1> &s, - const SymmetricTensor<4,1> &t) +double_contract (SymmetricTensor<2,1,Number> &tmp, + const SymmetricTensor<2,1,Number> &s, + const SymmetricTensor<4,1,Number> &t) { tmp[0][0] = t[0][0][0][0] * s[0][0]; } @@ -2854,11 +2776,12 @@ double_contract (SymmetricTensor<2,1> &tmp, * * @related SymmetricTensor @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,2> &tmp, - const SymmetricTensor<4,2> &t, - const SymmetricTensor<2,2> &s) +double_contract (SymmetricTensor<2,2,Number> &tmp, + const SymmetricTensor<4,2,Number> &t, + const SymmetricTensor<2,2,Number> &s) { const unsigned int dim = 2; @@ -2886,11 +2809,12 @@ double_contract (SymmetricTensor<2,2> &tmp, * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,2> &tmp, - const SymmetricTensor<2,2> &s, - const SymmetricTensor<4,2> &t) +double_contract (SymmetricTensor<2,2,Number> &tmp, + const SymmetricTensor<2,2,Number> &s, + const SymmetricTensor<4,2,Number> &t) { const unsigned int dim = 2; @@ -2918,11 +2842,12 @@ double_contract (SymmetricTensor<2,2> &tmp, * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,3> &tmp, - const SymmetricTensor<4,3> &t, - const SymmetricTensor<2,3> &s) +double_contract (SymmetricTensor<2,3,Number> &tmp, + const SymmetricTensor<4,3,Number> &t, + const SymmetricTensor<2,3,Number> &s) { const unsigned int dim = 3; @@ -2953,11 +2878,12 @@ double_contract (SymmetricTensor<2,3> &tmp, * @related SymmetricTensor * @author Wolfgang Bangerth, 2005 */ +template inline void -double_contract (SymmetricTensor<2,3> &tmp, - const SymmetricTensor<2,3> &s, - const SymmetricTensor<4,3> &t) +double_contract (SymmetricTensor<2,3,Number> &tmp, + const SymmetricTensor<2,3,Number> &s, + const SymmetricTensor<4,3,Number> &t) { const unsigned int dim = 3; @@ -2989,12 +2915,12 @@ double_contract (SymmetricTensor<2,3> &tmp, * @relates SymmetricTensor * @author Wolfgang Bangerth, 2005 */ -template -Tensor<1,dim> -operator * (const SymmetricTensor<2,dim> &src1, - const Tensor<1,dim> &src2) +template +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 &src1, * * @relates SymmetricTensor */ -template +template 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 +template 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 DEAL_II_NAMESPACE_OPEN -template class Tensor; -template class Tensor<1,dim>; +template class Tensor; +template class Tensor<1,dim,Number>; /** * Provide a general tensor class with an arbitrary rank, i.e. with @@ -33,10 +33,15 @@ template class Tensor<1,dim>; * produce far more efficient code than for matrices with * runtime-dependent dimension. * + * This class provides an optional template argument for the type of the + * underlying data. It defaults to @p double values. It can be used to base + * tensors on @p float or @p complex numbers or any other data type that + * implements basic arithmetic operations. + * * @ingroup geomprimitives * @author Wolfgang Bangerth, 1998-2005 */ -template +template class Tensor { public: @@ -69,7 +74,7 @@ class Tensor * Type of stored objects. This * is a tensor of lower rank. */ - typedef Tensor value_type; + typedef Tensor value_type; /** * Declare an array type which @@ -77,7 +82,7 @@ class Tensor * object of this type * statically. */ - typedef typename Tensor::array_type array_type[dim]; + typedef typename Tensor::array_type array_type[dim]; /** * Constructor. Initialize all entries @@ -94,17 +99,17 @@ class Tensor /** * Read-Write access operator. */ - Tensor &operator [] (const unsigned int i); + Tensor &operator [] (const unsigned int i); /** * Read-only access operator. */ - const Tensor &operator [] (const unsigned int i) const; + const Tensor &operator [] (const unsigned int i) const; /** * Assignment operator. */ - Tensor & operator = (const Tensor &); + Tensor & operator = (const Tensor &); /** * This operator assigns a scalar @@ -118,40 +123,40 @@ class Tensor * all elements of the tensor to * zero. */ - Tensor & operator = (const double d); + Tensor & operator = (const Number d); /** * Test for equality of two tensors. */ - bool operator == (const Tensor &) const; + bool operator == (const Tensor &) const; /** * Test for inequality of two tensors. */ - bool operator != (const Tensor &) const; + bool operator != (const Tensor &) const; /** * Add another tensor. */ - Tensor & operator += (const Tensor &); + Tensor & operator += (const Tensor &); /** * Subtract another tensor. */ - Tensor & operator -= (const Tensor &); + Tensor & operator -= (const Tensor &); /** * Scale the tensor by factor, * i.e. multiply all components by * factor. */ - Tensor & operator *= (const double factor); + Tensor & operator *= (const Number factor); /** * Scale the vector by * 1/factor. */ - Tensor & operator /= (const double factor); + Tensor & operator /= (const Number factor); /** * Add two tensors. If possible, you @@ -159,7 +164,7 @@ class Tensor * instead since this does not need the * creation of a temporary. */ - Tensor operator + (const Tensor &) const; + Tensor operator + (const Tensor &) const; /** * Subtract two tensors. If possible, @@ -167,20 +172,20 @@ class Tensor * instead since this does not need the * creation of a temporary. */ - Tensor operator - (const Tensor &) const; + Tensor operator - (const Tensor &) const; /** * Unary minus operator. Negate all * entries of a tensor. */ - Tensor operator - () const; + Tensor operator - () const; /** * Return the Frobenius-norm of a tensor, * i.e. the square root of the sum of * squares of all entries. */ - double norm () const; + Number norm () const; /** * Return the square of the @@ -194,7 +199,7 @@ class Tensor * may also be useful in other * contexts. */ - double norm_square () const; + Number norm_square () const; /** * Fill a vector with all tensor elements. @@ -205,7 +210,7 @@ class Tensor * usual in C++, the rightmost * index of the tensor marches fastest. */ - void unroll (Vector & result) const; + void unroll (Vector & result) const; /** @@ -254,12 +259,12 @@ class Tensor * Array of tensors holding the * subelements. */ - Tensor subtensor[dim]; + Tensor subtensor[dim]; /** * Help function for unroll. */ - void unroll_recursion(Vector &result, + void unroll_recursion(Vector &result, unsigned int &start_index) const; // make the following class a @@ -272,7 +277,7 @@ class Tensor // also, it would be sufficient to make // the function unroll_loops a friend, // but that seems to be impossible as well. - template friend class Tensor; + template friend class Tensor; }; @@ -280,9 +285,9 @@ class Tensor #ifndef DOXYGEN -template +template inline -Tensor::Tensor () +Tensor::Tensor () { // default constructor. not specifying an initializer list calls // the default constructor of the subobjects, which initialize them @@ -290,19 +295,19 @@ Tensor::Tensor () } -template +template inline -Tensor::Tensor (const array_type &initializer) +Tensor::Tensor (const array_type &initializer) { for (unsigned int i=0; i(initializer[i]); + subtensor[i] = Tensor(initializer[i]); } -template +template inline -typename Tensor::value_type& -Tensor::operator[] (const unsigned int i) +typename Tensor::value_type& +Tensor::operator[] (const unsigned int i) { Assert (i::operator[] (const unsigned int i) } -template +template inline -const typename Tensor::value_type& -Tensor::operator[] (const unsigned int i) const +const typename Tensor::value_type& +Tensor::operator[] (const unsigned int i) const { Assert (i::operator[] (const unsigned int i) const } -template +template inline -Tensor & -Tensor::operator = (const Tensor &t) +Tensor & +Tensor::operator = (const Tensor &t) { for (unsigned int i=0; i::operator = (const Tensor &t) } -template +template inline -Tensor & -Tensor::operator = (const double d) +Tensor & +Tensor::operator = (const Number d) { Assert (d==0, ExcMessage ("Only assignment with zero is allowed")); @@ -345,10 +350,10 @@ Tensor::operator = (const double d) } -template +template inline bool -Tensor::operator == (const Tensor &p) const +Tensor::operator == (const Tensor &p) const { for (unsigned int i=0; i::operator == (const Tensor &p) const } -template +template inline bool -Tensor::operator != (const Tensor &p) const +Tensor::operator != (const Tensor &p) const { return !((*this) == p); } -template +template inline -Tensor & -Tensor::operator += (const Tensor &p) +Tensor & +Tensor::operator += (const Tensor &p) { for (unsigned int i=0; i::operator += (const Tensor &p) } -template +template inline -Tensor & -Tensor::operator -= (const Tensor &p) +Tensor & +Tensor::operator -= (const Tensor &p) { for (unsigned int i=0; i::operator -= (const Tensor &p) } -template +template inline -Tensor & -Tensor::operator *= (const double s) +Tensor & +Tensor::operator *= (const Number s) { for (unsigned int i=0; i::operator *= (const double s) } -template +template inline -Tensor & -Tensor::operator /= (const double s) +Tensor & +Tensor::operator /= (const Number s) { for (unsigned int i=0; i::operator /= (const double s) } -template +template inline -Tensor -Tensor::operator + (const Tensor &t) const +Tensor +Tensor::operator + (const Tensor &t) const { - Tensor tmp(*this); + Tensor tmp(*this); for (unsigned int i=0; i::operator + (const Tensor &t) const } -template +template inline -Tensor -Tensor::operator - (const Tensor &t) const +Tensor +Tensor::operator - (const Tensor &t) const { - Tensor tmp(*this); + Tensor tmp(*this); for (unsigned int i=0; i::operator - (const Tensor &t) const } -template +template inline -Tensor -Tensor::operator - () const +Tensor +Tensor::operator - () const { - Tensor tmp; + Tensor tmp; for (unsigned int i=0; i::operator - () const } -template +template inline -double Tensor::norm () const +Number Tensor::norm () const { return std::sqrt (norm_square()); } -template +template inline -double Tensor::norm_square () const +Number Tensor::norm_square () const { - double s = 0; + Number s = 0; for (unsigned int i=0; i::norm_square () const } -template +template inline -void Tensor::clear () +void Tensor::clear () { for (unsigned int i=0; i +template inline std::size_t -Tensor::memory_consumption () +Tensor::memory_consumption () { - return sizeof(Tensor); + return sizeof(Tensor); } -template +template template inline void -Tensor::serialize(Archive & ar, const unsigned int) +Tensor::serialize(Archive & ar, const unsigned int) { ar & subtensor; } @@ -510,9 +515,9 @@ Tensor::serialize(Archive & ar, const unsigned int) * * @relates Tensor */ -template +template inline -std::ostream & operator << (std::ostream &out, const Tensor &p) +std::ostream & operator << (std::ostream &out, const Tensor &p) { for (unsigned int i=0; i &p) * @relates Tensor * @author Guido Kanschat, 2000 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template inline -double -operator * (const Tensor<1,dim> &src1, - const Tensor<1,dim> &src2) +Number +operator * (const Tensor<1,dim,Number> &src1, + const Tensor<1,dim,Number> &src2) { return contract(src1, src2); } @@ -591,11 +596,11 @@ operator * (const Tensor<1,dim> &src1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template 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 &dest, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template -Tensor<1,dim> -operator * (const Tensor<2,dim> &src1, - const Tensor<1,dim> &src2) +template +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 &src1, * @relates Tensor * @author Guido Kanschat, 2001 */ -template +template 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 &dest, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template 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 &dest, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template inline -void contract (Tensor<2,dim> &dest, - const Tensor<2,dim> &src1, const unsigned int index1, - const Tensor<2,dim> &src2, const unsigned int index2) +void contract (Tensor<2,dim,Number> &dest, + const Tensor<2,dim,Number> &src1, const unsigned int index1, + const Tensor<2,dim,Number> &src2, const unsigned int index2) { dest.clear (); @@ -777,7 +782,7 @@ void contract (Tensor<2,dim> &dest, default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2))); }; break; case 2: @@ -798,12 +803,12 @@ void contract (Tensor<2,dim> &dest, default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2))); }; break; default: - Assert (false, (typename Tensor<2,dim>::ExcInvalidTensorIndex (index1))); + Assert (false, (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1))); }; } @@ -819,11 +824,11 @@ void contract (Tensor<2,dim> &dest, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template inline -void contract (Tensor<2,dim> &dest, - const Tensor<3,dim> &src1, const unsigned int index1, - const Tensor<1,dim> &src2) +void contract (Tensor<2,dim,Number> &dest, + const Tensor<3,dim,Number> &src1, const unsigned int index1, + const Tensor<1,dim,Number> &src2) { dest.clear (); @@ -852,7 +857,7 @@ void contract (Tensor<2,dim> &dest, default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index1))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index1))); }; } @@ -864,11 +869,11 @@ void contract (Tensor<2,dim> &dest, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template 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 &dest, * * @relates Tensor */ -template +template inline -void contract (Tensor<3,dim> &dest, - const Tensor<3,dim> &src1, const unsigned int index1, - const Tensor<2,dim> &src2, const unsigned int index2) +void contract (Tensor<3,dim,Number> &dest, + const Tensor<3,dim,Number> &src1, const unsigned int index1, + const Tensor<2,dim,Number> &src2, const unsigned int index2) { dest.clear (); @@ -921,7 +926,7 @@ void contract (Tensor<3,dim> &dest, break; default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2))); } break; @@ -944,7 +949,7 @@ void contract (Tensor<3,dim> &dest, break; default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2))); } break; @@ -967,13 +972,13 @@ void contract (Tensor<3,dim> &dest, break; default: Assert (false, - (typename Tensor<2,dim>::ExcInvalidTensorIndex (index2))); + (typename Tensor<2,dim,Number>::ExcInvalidTensorIndex (index2))); } break; default: Assert (false, - (typename Tensor<3,dim>::ExcInvalidTensorIndex (index1))); + (typename Tensor<3,dim,Number>::ExcInvalidTensorIndex (index1))); } } @@ -994,13 +999,13 @@ void contract (Tensor<3,dim> &dest, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template 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 &dest, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template 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 &src1, * @relates Tensor * @author Wolfgang Bangerth, 2005 */ -template +template 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 &dest, * @relates Tensor * @author Guido Kanschat, 2004 */ -template +template 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& u, * @relates Tensor * @author Toby D. Young, 2011 */ -template +template 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 &t1, * @relates Tensor * @author Toby D. Young, 2011 */ -template +template 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 &t1, * @relates Tensor * @author Toby D. Young, 2011 */ -template +template 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 &t1, * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template -void outer_product (Tensor<2,dim> &dst, - const Tensor<1,dim> &src1, - const Tensor<1,dim> &src2) +template +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 &dst, * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template -void outer_product (Tensor<3,dim> &dst, - const Tensor<1,dim> &src1, - const Tensor<2,dim> &src2) +template +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 &dst, * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template -void outer_product (Tensor<3,dim> &dst, - const Tensor<2,dim> &src1, - const Tensor<1,dim> &src2) +template +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 &dst, * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template -void outer_product (Tensor<1,dim> &dst, - const double src1, - const Tensor<1,dim> &src2) +template +void outer_product (Tensor<1,dim,Number> &dst, + const Number src1, + const Tensor<1,dim,Number> &src2) { for (unsigned int i=0; i &dst, * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template -void outer_product (Tensor<1,dim> &dst, - const Tensor<1,dim> src1, - const double src2) +template +void outer_product (Tensor<1,dim,Number> &dst, + const Tensor<1,dim,Number> src1, + const Number src2) { for (unsigned int i=0; i &dst, * @relates Tensor * @author Guido Kanschat, 2001 */ -template +template inline void -cross_product (Tensor<1,dim> &dst, - const Tensor<1,dim> &src) +cross_product (Tensor<1,dim,Number> &dst, + const Tensor<1,dim,Number> &src) { - Assert (dim==2, ExcImpossibleInDim(dim)); + Assert (dim==2, ExcInternalError()); dst[0] = src[1]; dst[1] = -src[0]; @@ -1337,14 +1342,14 @@ cross_product (Tensor<1,dim> &dst, * @relates Tensor * @author Guido Kanschat, 2001 */ -template +template inline void -cross_product (Tensor<1,dim> &dst, - const Tensor<1,dim> &src1, - const Tensor<1,dim> &src2) +cross_product (Tensor<1,dim,Number> &dst, + const Tensor<1,dim,Number> &src1, + const Tensor<1,dim,Number> &src2) { - Assert (dim==3, ExcImpossibleInDim(dim)); + Assert (dim==3, ExcInternalError()); dst[0] = src1[1]*src2[2] - src1[2]*src2[1]; dst[1] = src1[2]*src2[0] - src1[0]*src2[2]; @@ -1362,13 +1367,13 @@ cross_product (Tensor<1,dim> &dst, * @relates Tensor * @author Wolfgang Bangerth, 2008 */ -template +template 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 &t1, * @relates Tensor * @author Wolfgang Bangerth, 1998 */ -template +template inline -double determinant (const Tensor &t) +Number determinant (const Tensor &t) { return determinant(t[0]); } @@ -1401,8 +1406,9 @@ double determinant (const Tensor &t) * @relates Tensor * @author Wolfgang Bangerth, 1998 */ +template inline -double determinant (const Tensor<1,1> &t) +Number determinant (const Tensor<1,1,Number> &t) { return t[0]; } @@ -1416,8 +1422,9 @@ double determinant (const Tensor<1,1> &t) * @relates Tensor * @author Wolfgang Bangerth, 1998 */ +template inline -double determinant (const Tensor<2,1> &t) +Number determinant (const Tensor<2,1,Number> &t) { return t[0][0]; } @@ -1430,8 +1437,9 @@ double determinant (const Tensor<2,1> &t) * @relates Tensor * @author Wolfgang Bangerth, 1998 */ +template inline -double determinant (const Tensor<2,2> &t) +Number determinant (const Tensor<2,2,Number> &t) { return ((t[0][0] * t[1][1]) - (t[1][0] * t[0][1])); @@ -1444,8 +1452,9 @@ double determinant (const Tensor<2,2> &t) * @relates Tensor * @author Wolfgang Bangerth, 1998 */ +template inline -double determinant (const Tensor<2,3> &t) +Number determinant (const Tensor<2,3,Number> &t) { // get this using Maple: // with(linalg); @@ -1469,9 +1478,9 @@ double determinant (const Tensor<2,3> &t) * @relates Tensor * @author Wolfgang Bangerth, 2009 */ -template +template inline -double determinant (const Tensor<2,dim> &t) +Number determinant (const Tensor<2,dim,Number> &t) { // compute the determinant using the // Laplace expansion of the @@ -1481,7 +1490,7 @@ double determinant (const Tensor<2,dim> &t) // // for some algorithmic simplicity, we use // the expansion along the last row - double det = 0; + Number det = 0; for (unsigned int k=0; k &t) for (unsigned int j=0; j(k+1)) * + const Number cofactor = std::pow (-1., static_cast(k+1)) * determinant (minor); det += t[dim-1][k] * cofactor; } - return std::pow (-1., static_cast(dim)) * det; + return std::pow (-1., static_cast(dim)) * det; } @@ -1508,10 +1517,10 @@ double determinant (const Tensor<2,dim> &t) * @relates Tensor * @author Wolfgang Bangerth, 2001 */ -template -double trace (const Tensor<2,dim> &d) +template +Number trace (const Tensor<2,dim,Number> &d) { - double t=0; + Number t=0; for (unsigned int i=0; i &d) * @relates Tensor * @author Wolfgang Bangerth, 2000 */ -template +template inline -Tensor<2,dim> -invert (const Tensor<2,dim> &t) +Tensor<2,dim,Number> +invert (const Tensor<2,dim,Number> &t) { - Tensor<2,dim> return_tensor; + Tensor<2,dim,Number> return_tensor; switch (dim) { case 1: @@ -1544,7 +1553,7 @@ invert (const Tensor<2,dim> &t) // this is Maple output, // thus a bit unstructured { - const double t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]); + const Number t4 = 1.0/(t[0][0]*t[1][1]-t[0][1]*t[1][0]); return_tensor[0][0] = t[1][1]*t4; return_tensor[0][1] = -t[0][1]*t4; return_tensor[1][0] = -t[1][0]*t4; @@ -1554,7 +1563,7 @@ invert (const Tensor<2,dim> &t) case 3: { - const double t4 = t[0][0]*t[1][1], + const Number t4 = t[0][0]*t[1][1], t6 = t[0][0]*t[1][2], t8 = t[0][1]*t[1][0], t00 = t[0][2]*t[1][0], @@ -1597,16 +1606,16 @@ invert (const Tensor<2,dim> &t) * @relates Tensor * @author Wolfgang Bangerth, 2002 */ -template +template 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 &t) * @relates Tensor * @author Wolfgang Bangerth, 2002 */ +template inline -Tensor<2,1> -transpose (const Tensor<2,1> &t) +Tensor<2,1,Number> +transpose (const Tensor<2,1,Number> &t) { return t; } @@ -1639,12 +1649,13 @@ transpose (const Tensor<2,1> &t) * @relates Tensor * @author Wolfgang Bangerth, 2002 */ +template inline -Tensor<2,2> -transpose (const Tensor<2,2> &t) +Tensor<2,2,Number> +transpose (const Tensor<2,2,Number> &t) { - const double x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}}; - return Tensor<2,2>(x); + const Number x[2][2] = {{t[0][0], t[1][0]}, {t[0][1], t[1][1]}}; + return Tensor<2,2,Number>(x); } @@ -1657,19 +1668,78 @@ transpose (const Tensor<2,2> &t) * @relates Tensor * @author Wolfgang Bangerth, 2002 */ +template 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 +inline +Tensor +operator * (const Tensor &t, + const Number factor) +{ + Tensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Multiplication of a tensor of general rank with a scalar Number + * from the left. + * + * @relates Tensor + */ +template +inline +Tensor +operator * (const Number factor, + const Tensor &t) +{ + Tensor tt = t; + tt *= factor; + return tt; +} + + + +/** + * Division of a tensor of general rank by a scalar Number. + * + * @relates Tensor + */ +template +inline +Tensor +operator / (const Tensor &t, + const Number factor) +{ + Tensor tt = t; + tt /= factor; + return tt; +} + + + + /** * Multiplication of a tensor of general rank with a scalar double * from the right. diff --git a/deal.II/include/deal.II/base/tensor_base.h b/deal.II/include/deal.II/base/tensor_base.h index 9a2e511bf8..bbfc19eb9a 100644 --- a/deal.II/include/deal.II/base/tensor_base.h +++ b/deal.II/include/deal.II/base/tensor_base.h @@ -13,10 +13,10 @@ #define __deal2__tensor_base_h -// single this file out from tensor.h, since we want to derive Point -// 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 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 @@ -35,33 +35,39 @@ DEAL_II_NAMESPACE_OPEN #endif template class Vector; -template 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 class Point; // general template; specialized for rank==1; the general template is in // tensor.h -template class Tensor; -template class Tensor<0,dim>; -template class Tensor<1,dim>; +template class Tensor; +template class Tensor<0,dim,Number>; +template class Tensor<1,dim,Number>; /** - * This class is a specialized version of the - * Tensor class. It handles tensors of rank zero, - * i.e. scalars. The second template argument is ignored. + * This class is a specialized version of the Tensor + * class. It handles tensors of rank zero, i.e. scalars. The second template + * argument is ignored. * * This class exists because in some cases we want to construct - * objects of type Tensor@, which should expand to + * objects of type Tensor@, which should expand to * scalars, vectors, matrices, etc, depending on the values of the * template arguments @p dim and @p spacedim. We therefore need a - * class that acts as a scalar (i.e. @p double) for all purposes but + * 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 -class Tensor<0,dim> +template +class Tensor<0,dim,Number> { public: /** @@ -91,11 +97,11 @@ class Tensor<0,dim> /** * Type of stored objects. This - * is a double for a rank 1 tensor. + * is a Number for a rank 1 tensor. */ - typedef double value_type; - + typedef Number value_type; + /** * Constructor. Set to zero. */ @@ -107,21 +113,21 @@ class Tensor<0,dim> * array. */ Tensor (const value_type &initializer); - + /** * Copy constructor. */ - Tensor (const Tensor<0,dim> &); + Tensor (const Tensor<0,dim,Number> &); /** - * Conversion to double. Since + * Conversion to Number. Since * rank-0 tensors are scalars, * this is a natural operation. */ - operator double () const; + operator Number () const; /** - * Conversion to double. Since + * Conversion to Number. Since * rank-0 tensors are scalars, * this is a natural operation. * @@ -129,59 +135,59 @@ class Tensor<0,dim> * conversion operator that * returns a writable reference. */ - operator double& (); + operator Number& (); /** * Assignment operator. */ - Tensor<0,dim> & operator = (const Tensor<0,dim> &); + Tensor<0,dim,Number> & operator = (const Tensor<0,dim,Number> &); /** * Assignment operator. */ - Tensor<0,dim> & operator = (const double d); + Tensor<0,dim,Number> & operator = (const Number d); /** * Test for equality of two * tensors. */ - bool operator == (const Tensor<0,dim> &) const; + bool operator == (const Tensor<0,dim,Number> &) const; /** * Test for inequality of two * tensors. */ - bool operator != (const Tensor<0,dim> &) const; + bool operator != (const Tensor<0,dim,Number> &) const; /** * Add another vector, i.e. move * this point by the given * offset. */ - Tensor<0,dim> & operator += (const Tensor<0,dim> &); - + Tensor<0,dim,Number> & operator += (const Tensor<0,dim,Number> &); + /** * Subtract another vector. */ - Tensor<0,dim> & operator -= (const Tensor<0,dim> &); + Tensor<0,dim,Number> & operator -= (const Tensor<0,dim,Number> &); /** * Scale the vector by * factor, i.e. multiply all * coordinates by factor. */ - Tensor<0,dim> & operator *= (const double factor); + Tensor<0,dim,Number> & operator *= (const Number factor); /** * Scale the vector by 1/factor. */ - Tensor<0,dim> & operator /= (const double factor); + Tensor<0,dim,Number> & operator /= (const Number factor); /** * Returns the scalar product of * two vectors. */ - double operator * (const Tensor<0,dim> &) const; + Number operator * (const Tensor<0,dim,Number> &) const; /** * Add two tensors. If possible, @@ -189,7 +195,7 @@ class Tensor<0,dim> * since this does not need to * copy a point at least once. */ - Tensor<0,dim> operator + (const Tensor<0,dim> &) const; + Tensor<0,dim,Number> operator + (const Tensor<0,dim,Number> &) const; /** * Subtract two tensors. If @@ -198,13 +204,13 @@ class Tensor<0,dim> * need to copy a point at least * once. */ - Tensor<0,dim> operator - (const Tensor<0,dim> &) const; + Tensor<0,dim,Number> operator - (const Tensor<0,dim,Number> &) const; /** * Tensor with inverted entries. */ - Tensor<0,dim> operator - () const; - + Tensor<0,dim,Number> operator - () const; + /** * Return the Frobenius-norm of a * tensor, i.e. the square root @@ -215,7 +221,7 @@ class Tensor<0,dim> * l2 norm of * the vector. */ - double norm () const; + Number norm () const; /** * Return the square of the @@ -229,7 +235,7 @@ class Tensor<0,dim> * may also be useful in other * contexts. */ - double norm_square () const; + Number norm_square () const; /** * Reset all values to zero. @@ -266,9 +272,9 @@ class Tensor<0,dim> << "dim must be positive, but was " << arg1); /** - * Read or write the data of this object to or + * Read or write the data of this object to or * from a stream for the purpose of serialization - */ + */ template void serialize(Archive & ar, const unsigned int version); @@ -276,33 +282,35 @@ class Tensor<0,dim> /** * The value of this scalar object. */ - double value; + Number value; }; + /** - * This class is a specialized version of the Tensor 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 Tensor + * 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 * Point is that we use the Point 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 - * Tensor<1,dim> class. In particular, there is a distance() function - * to compute the Euclidian distance between two points in space. + * Tensor<1,dim,Number> class. In particular, there is a distance() + * function to compute the Euclidian distance between two points in space. * * However, the Point 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 Tensor<1,dim> class. + * different physical units), we use the Tensor<1,dim,Number> class. * * @ingroup geomprimitives * @author Wolfgang Bangerth, 1998-2005 */ -template -class Tensor<1,dim> +template +class Tensor<1,dim,Number> { public: /** @@ -332,11 +340,11 @@ class Tensor<1,dim> /** * Type of stored objects. This - * is a double for a rank 1 tensor. + * is a Number for a rank 1 tensor. */ - typedef double value_type; - + typedef Number value_type; + /** * Declare an array type which can * be used to initialize statically @@ -348,8 +356,8 @@ class Tensor<1,dim> * likely to overflow memory * limits. */ - typedef double array_type[(dim!=0) ? dim : 100000000]; - + typedef Number array_type[(dim!=0) ? dim : 100000000]; + /** * Constructor. Initialize all entries * to zero if initialize==true; this @@ -363,11 +371,11 @@ class Tensor<1,dim> * array. */ Tensor (const array_type &initializer); - + /** * Copy constructor. */ - Tensor (const Tensor<1,dim> &); + Tensor (const Tensor<1,dim,Number> &); /** * Read access to the indexth @@ -379,7 +387,7 @@ class Tensor<1,dim> * operator for * backcompatibility. */ - double operator [] (const unsigned int index) const; + Number operator [] (const unsigned int index) const; /** * Read and write access to the @@ -391,12 +399,12 @@ class Tensor<1,dim> * operator for * backcompatibility. */ - double & operator [] (const unsigned int index); + Number & operator [] (const unsigned int index); /** * Assignment operator. */ - Tensor<1,dim> & operator = (const Tensor<1,dim> &); + Tensor<1,dim,Number> & operator = (const Tensor<1,dim,Number> &); /** * This operator assigns a scalar @@ -410,49 +418,49 @@ class Tensor<1,dim> * all elements of the tensor to * zero. */ - Tensor<1,dim> & operator = (const double d); + Tensor<1,dim,Number> & operator = (const Number d); /** * Test for equality of two * tensors. */ - bool operator == (const Tensor<1,dim> &) const; + bool operator == (const Tensor<1,dim,Number> &) const; /** * Test for inequality of two * tensors. */ - bool operator != (const Tensor<1,dim> &) const; + bool operator != (const Tensor<1,dim,Number> &) const; /** * Add another vector, i.e. move * this point by the given * offset. */ - Tensor<1,dim> & operator += (const Tensor<1,dim> &); - + Tensor<1,dim,Number> & operator += (const Tensor<1,dim,Number> &); + /** * Subtract another vector. */ - Tensor<1,dim> & operator -= (const Tensor<1,dim> &); + Tensor<1,dim,Number> & operator -= (const Tensor<1,dim,Number> &); /** * Scale the vector by * factor, i.e. multiply all * coordinates by factor. */ - Tensor<1,dim> & operator *= (const double factor); + Tensor<1,dim,Number> & operator *= (const Number factor); /** * Scale the vector by 1/factor. */ - Tensor<1,dim> & operator /= (const double factor); + Tensor<1,dim,Number> & operator /= (const Number factor); /** * Returns the scalar product of * two vectors. */ - double operator * (const Tensor<1,dim> &) const; + Number operator * (const Tensor<1,dim,Number> &) const; /** * Add two tensors. If possible, @@ -460,7 +468,7 @@ class Tensor<1,dim> * since this does not need to * copy a point at least once. */ - Tensor<1,dim> operator + (const Tensor<1,dim> &) const; + Tensor<1,dim,Number> operator + (const Tensor<1,dim,Number> &) const; /** * Subtract two tensors. If @@ -469,13 +477,13 @@ class Tensor<1,dim> * need to copy a point at least * once. */ - Tensor<1,dim> operator - (const Tensor<1,dim> &) const; + Tensor<1,dim,Number> operator - (const Tensor<1,dim,Number> &) const; /** * Tensor with inverted entries. */ - Tensor<1,dim> operator - () const; - + Tensor<1,dim,Number> operator - () const; + /** * Return the Frobenius-norm of a * tensor, i.e. the square root @@ -486,7 +494,7 @@ class Tensor<1,dim> * l2 norm of * the vector. */ - double norm () const; + Number norm () const; /** * Return the square of the @@ -500,7 +508,7 @@ class Tensor<1,dim> * may also be useful in other * contexts. */ - double norm_square () const; + Number norm_square () const; /** * Reset all values to zero. @@ -531,7 +539,7 @@ class Tensor<1,dim> * usual in C++, the rightmost * index marches fastest. */ - void unroll (Vector &result) const; + void unroll (Vector &result) const; /** * Determine an estimate for @@ -554,14 +562,14 @@ class Tensor<1,dim> DeclException1 (ExcDimTooSmall, int, << "dim must be positive, but was " << arg1); - + /** - * Read or write the data of this object to or + * Read or write the data of this object to or * from a stream for the purpose of serialization - */ + */ template void serialize(Archive & ar, const unsigned int version); - + private: /** * Store the values in a simple @@ -572,7 +580,7 @@ class Tensor<1,dim> * constructor to disallow the * creation of such an object. */ - double values[(dim!=0) ? (dim) : 1]; + Number values[(dim!=0) ? (dim) : 1]; #ifdef DEAL_II_TEMPLATE_SPEC_ACCESS_WORKAROUND public: @@ -588,9 +596,9 @@ class Tensor<1,dim> * case, even if it should be * public for your compiler. */ - void unroll_recursion (Vector &result, + void unroll_recursion (Vector &result, unsigned int &start_index) const; - + private: /** * Make the following classes @@ -605,7 +613,7 @@ class Tensor<1,dim> * that seems to be impossible as * well. */ - template friend class dealii::Tensor; + template friend class dealii::Tensor; /** * Point is allowed access to @@ -619,24 +627,24 @@ class Tensor<1,dim> /** * Prints the value of this scalar. */ -template -std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p); +template +std::ostream & operator << (std::ostream &out, const Tensor<0,dim,Number> &p); /** * Prints the values of this tensor in the * form x1 x2 x3 etc. */ -template -std::ostream & operator << (std::ostream &out, const Tensor<1,dim> &p); +template +std::ostream & operator << (std::ostream &out, const Tensor<1,dim,Number> &p); #ifndef DOXYGEN /*---------------------------- Inline functions: Tensor<0,dim> ------------------------*/ -template +template inline -Tensor<0,dim>::Tensor () +Tensor<0,dim,Number>::Tensor () { Assert (dim>0, ExcDimTooSmall(dim)); @@ -645,9 +653,9 @@ Tensor<0,dim>::Tensor () -template +template inline -Tensor<0,dim>::Tensor (const value_type &initializer) +Tensor<0,dim,Number>::Tensor (const value_type &initializer) { Assert (dim>0, ExcDimTooSmall(dim)); @@ -656,39 +664,39 @@ Tensor<0,dim>::Tensor (const value_type &initializer) -template +template 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 +template inline -Tensor<0,dim>::operator double () const +Tensor<0,dim,Number>::operator Number () const { return value; } -template +template inline -Tensor<0,dim>::operator double & () +Tensor<0,dim,Number>::operator Number & () { return value; } -template +template inline -Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Tensor<0,dim,Number> &p) { value = p.value; return *this; @@ -696,9 +704,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator = (const Tensor<0,dim> &p) -template +template inline -Tensor<0,dim> & Tensor<0,dim>::operator = (const double d) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator = (const Number d) { value = d; return *this; @@ -706,27 +714,27 @@ Tensor<0,dim> & Tensor<0,dim>::operator = (const double d) -template +template 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 +template 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 +template inline -Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator += (const Tensor<0,dim,Number> &p) { value += p.value; return *this; @@ -734,9 +742,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator += (const Tensor<0,dim> &p) -template +template inline -Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Number> &p) { value -= p.value; return *this; @@ -744,9 +752,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator -= (const Tensor<0,dim> &p) -template +template inline -Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator *= (const Number s) { value *= s; return *this; @@ -754,9 +762,9 @@ Tensor<0,dim> & Tensor<0,dim>::operator *= (const double s) -template +template inline -Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s) +Tensor<0,dim,Number> & Tensor<0,dim,Number>::operator /= (const Number s) { value /= s; return *this; @@ -764,81 +772,81 @@ Tensor<0,dim> & Tensor<0,dim>::operator /= (const double s) -template +template 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 +template 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 +template 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 +template inline -Tensor<0,dim> Tensor<0,dim>::operator - () const +Tensor<0,dim,Number> Tensor<0,dim,Number>::operator - () const { return -value; } -template +template inline -double Tensor<0,dim>::norm () const +Number Tensor<0,dim,Number>::norm () const { return std::abs (value); } -template +template inline -double Tensor<0,dim>::norm_square () const +Number Tensor<0,dim,Number>::norm_square () const { return value*value; } -template +template inline -void Tensor<0,dim>::clear () +void Tensor<0,dim,Number>::clear () { value = 0; } -template +template template 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 +template inline -Tensor<1,dim>::Tensor (const bool initialize) +Tensor<1,dim,Number>::Tensor (const bool initialize) { Assert (dim>0, ExcDimTooSmall(dim)); @@ -849,9 +857,9 @@ Tensor<1,dim>::Tensor (const bool initialize) -template +template inline -Tensor<1,dim>::Tensor (const array_type &initializer) +Tensor<1,dim,Number>::Tensor (const array_type &initializer) { Assert (dim>0, ExcDimTooSmall(dim)); @@ -861,12 +869,12 @@ Tensor<1,dim>::Tensor (const array_type &initializer) -template +template 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::Tensor (const Tensor<1,dim> &p) template <> inline -Tensor<1,0>::Tensor (const Tensor<1,0> &) +Tensor<1,0,double>::Tensor (const Tensor<1,0,double> &) { // at some places in the library, // we have Point<0> for formal @@ -891,9 +899,9 @@ Tensor<1,0>::Tensor (const Tensor<1,0> &) -template +template inline -double Tensor<1,dim>::operator [] (const unsigned int index) const +Number Tensor<1,dim,Number>::operator [] (const unsigned int index) const { Assert (index::operator [] (const unsigned int index) const -template +template inline -double & Tensor<1,dim>::operator [] (const unsigned int index) +Number & Tensor<1,dim,Number>::operator [] (const unsigned int index) { Assert (index::operator [] (const unsigned int index) template <> inline -Tensor<1,0> & Tensor<1,0>::operator = (const Tensor<1,0> &) +Tensor<1,0,double> & Tensor<1,0,double>::operator = (const Tensor<1,0,double> &) { // at some places in the library, // we have Point<0> for formal @@ -930,82 +938,56 @@ Tensor<1,0> & Tensor<1,0>::operator = (const Tensor<1,0> &) -template <> +template 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 +template 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 -inline -Tensor<1,dim> & Tensor<1,dim>::operator = (const Tensor<1,dim> &p) -{ - for (unsigned int i=0; i -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 +template 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::operator == (const Tensor<1,dim> &p) const template <> inline -bool Tensor<1,0>::operator == (const Tensor<1,0> &) const +bool Tensor<1,0,double>::operator == (const Tensor<1,0,double> &) const { return true; } -template +template 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 +template 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 & Tensor<1,dim>::operator += (const Tensor<1,dim> &p) -template +template 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 & Tensor<1,dim>::operator -= (const Tensor<1,dim> &p) -template +template 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 & Tensor<1,dim>::operator *= (const double s) -template +template 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 & Tensor<1,dim>::operator /= (const double s) -template <> +template 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 +template 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 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 +template 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 -inline -Tensor<1,dim> Tensor<1,dim>::operator + (const Tensor<1,dim> &p) const -{ - return (Tensor<1,dim>(*this) += p); -} - - - -template -inline -Tensor<1,dim> Tensor<1,dim>::operator - (const Tensor<1,dim> &p) const -{ - return (Tensor<1,dim>(*this) -= p); -} - - - -template -inline -Tensor<1,dim> Tensor<1,dim>::operator - () const -{ - Tensor<1,dim> result; + Tensor<1,dim,Number> result; for (unsigned int i=0; i Tensor<1,dim>::operator - () const -template +template inline -double Tensor<1,dim>::norm () const +Number Tensor<1,dim,Number>::norm () const { return std::sqrt (norm_square()); } -template +template 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::norm_square () const -template +template inline -void Tensor<1,dim>::clear () +void Tensor<1,dim,Number>::clear () { for (unsigned int i=0; i::clear () -template +template 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 +template template inline -void Tensor<1,dim>::serialize(Archive & ar, const unsigned int) +void Tensor<1,dim,Number>::serialize(Archive & ar, const unsigned int) { ar & values; } @@ -1219,13 +1178,13 @@ void Tensor<1,dim>::serialize(Archive & ar, const unsigned int) * Output operator for tensors of rank 0. Since such tensors are * scalars, we simply print this one value. * - * @relates Tensor<0,dim> + * @relates Tensor<0,dim,Number> */ -template +template 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(p); + out << static_cast(p); return out; } @@ -1235,11 +1194,11 @@ std::ostream & operator << (std::ostream &out, const Tensor<0,dim> &p) * Output operator for tensors of rank 1. Print the elements * consecutively, with a space in between. * - * @relates Tensor<1,dim> + * @relates Tensor<1,dim,Number> */ -template +template 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 &p) * implemented specialized from the general template in order to avoid * a compiler warning that the loop is empty. * - * @relates Tensor<1,dim> + * @relates Tensor<1,dim,Number> */ inline -std::ostream & operator << (std::ostream &out, const Tensor<1,1> &p) +std::ostream & operator << (std::ostream &out, const Tensor<1,1,double> &p) { out << p[0]; @@ -1267,10 +1226,67 @@ std::ostream & operator << (std::ostream &out, const Tensor<1,1> &p) +/** + * Multiplication of a tensor of rank 1 with a scalar Number from the right. + * + * @relates Tensor<1,dim,Number> + */ +template +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 + */ +template +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 + */ +template +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 + * @relates Tensor<1,dim,Number> */ template inline @@ -1289,7 +1305,7 @@ operator * (const Tensor<1,dim> &t, /** * Multiplication of a tensor of rank 1 with a scalar double from the left. * - * @relates Tensor<1,dim> + * @relates Tensor<1,dim,Number> */ template inline @@ -1308,7 +1324,7 @@ operator * (const double factor, /** * Division of a tensor of rank 1 by a scalar double. * - * @relates Tensor<1,dim> + * @relates Tensor<1,dim,Number> */ template inline @@ -1321,6 +1337,8 @@ operator / (const Tensor<1,dim> &t, tt[d] = t[d] / factor; return tt; } + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 2f9b90f2b3..61e71f093c 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -32,7 +33,6 @@ template class Table; class SparsityPattern; template class Vector; template class Function; -template class Point; template class FiniteElement; template class DoFHandler; namespace hp diff --git a/deal.II/include/deal.II/fe/fe_tools.h b/deal.II/include/deal.II/fe/fe_tools.h index 9fc27fb6ac..e93b1ec0da 100644 --- a/deal.II/include/deal.II/fe/fe_tools.h +++ b/deal.II/include/deal.II/fe/fe_tools.h @@ -17,6 +17,8 @@ #include #include #include +#include +#include #include @@ -24,8 +26,6 @@ DEAL_II_NAMESPACE_OPEN template class FullMatrix; template class Vector; -template class SymmetricTensor; -template class Tensor; template class Quadrature; template class FiniteElement; template class DoFHandler; diff --git a/deal.II/include/deal.II/grid/grid_in.h b/deal.II/include/deal.II/grid/grid_in.h index 9c64767dda..5f34bb8943 100644 --- a/deal.II/include/deal.II/grid/grid_in.h +++ b/deal.II/include/deal.II/grid/grid_in.h @@ -16,13 +16,13 @@ #include #include #include +#include #include #include #include DEAL_II_NAMESPACE_OPEN -template class Point; template class Triangulation; template class CellData; class SubCellData; diff --git a/deal.II/include/deal.II/grid/tria_accessor.h b/deal.II/include/deal.II/grid/tria_accessor.h index 5061bbadec..426ed34b51 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.h +++ b/deal.II/include/deal.II/grid/tria_accessor.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -28,8 +29,6 @@ namespace std DEAL_II_NAMESPACE_OPEN -template class Point; - template class Triangulation; template class TriaRawIterator; template class TriaIterator; diff --git a/deal.II/include/deal.II/lac/full_matrix.h b/deal.II/include/deal.II/lac/full_matrix.h index 94f3f70e19..1f61a511a8 100644 --- a/deal.II/include/deal.II/lac/full_matrix.h +++ b/deal.II/include/deal.II/lac/full_matrix.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -29,7 +30,6 @@ DEAL_II_NAMESPACE_OPEN // forward declarations template class Vector; template class LAPACKFullMatrix; -template class Tensor; /*! @addtogroup Matrix1 diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h index f1b4532938..3b74811600 100644 --- a/deal.II/include/deal.II/numerics/vectors.h +++ b/deal.II/include/deal.II/numerics/vectors.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -26,7 +27,6 @@ DEAL_II_NAMESPACE_OPEN -template class Point; template class Function; template class FunctionMap; template class Quadrature; @@ -1173,8 +1173,8 @@ class VectorTools const unsigned char boundary_component, ConstraintMatrix& constraints, const hp::MappingCollection& mapping_collection = hp::StaticMappingQ1::mapping_collection); - - + + /** * Compute the constraints that * correspond to boundary conditions of diff --git a/deal.II/source/base/config.cc b/deal.II/source/base/config.cc index e142d3be6b..2a8f0043b4 100644 --- a/deal.II/source/base/config.cc +++ b/deal.II/source/base/config.cc @@ -14,8 +14,8 @@ #include #include -#include #include +#include DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/source/base/function_parser.cc b/deal.II/source/base/function_parser.cc index c4a895f610..3b08e510be 100644 --- a/deal.II/source/base/function_parser.cc +++ b/deal.II/source/base/function_parser.cc @@ -13,7 +13,6 @@ #include #include -#include #include #ifndef DEAL_II_DISABLE_PARSER @@ -38,7 +37,7 @@ FunctionParser::FunctionParser(const unsigned int n_components, : Function(n_components, initial_time), fp (0) -{ +{ fp = new fparser::FunctionParser[n_components]; } @@ -64,7 +63,7 @@ void FunctionParser::initialize (const std::string &varia constants, std::map< std::string, double >(), time_dependent, - use_degrees); + use_degrees); } @@ -104,7 +103,7 @@ void FunctionParser::initialize (const std::string &variables, } - + // Add the various constants to // the parser. std::map< std::string, double >::const_iterator @@ -115,21 +114,21 @@ void FunctionParser::initialize (const std::string &variables, const bool success = fp[i].AddConstant(constant->first, constant->second); AssertThrow (success, ExcMessage("Invalid Constant Name")); } - - const int ret_value = fp[i].Parse(expressions[i], - variables, + + const int ret_value = fp[i].Parse(expressions[i], + variables, use_degrees); - AssertThrow (ret_value == -1, + AssertThrow (ret_value == -1, ExcParseError(ret_value, fp[i].ErrorMsg())); - + // The fact that the parser did // not throw an error does not // mean that everything went // ok... we can still have // problems with the number of // variables... - } - + } + // Now we define how many variables // we expect to read in. We // distinguish between two cases: @@ -143,16 +142,16 @@ void FunctionParser::initialize (const std::string &variables, // parsed the variables string, if // none of this is the case, then // an exception is thrown. - if (time_dependent) + if (time_dependent) n_vars = dim+1; - else + else n_vars = dim; - + // Let's check if the number of // variables is correct... - AssertThrow (n_vars == fp[0].NVars(), + AssertThrow (n_vars == fp[0].NVars(), ExcDimensionMismatch(n_vars,fp[0].NVars())); - + // Now set the initialization bit. initialized = true; } @@ -172,7 +171,7 @@ void FunctionParser::initialize (const std::string &variables, Utilities::split_string_list (expression, ';'), constants, time_dependent, - use_degrees); + use_degrees); } @@ -199,19 +198,19 @@ void FunctionParser::initialize (const std::string &variables, template double FunctionParser::value (const Point &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::value (const Point &p, // have already been thrown if (dim != n_vars) vars[dim] = this->get_time(); - + double my_value = fp[component].Eval((double*)vars); - + AssertThrow (fp[component].EvalError() == 0, ExcMessage(fp[component].ErrorMsg())); return my_value; @@ -232,19 +231,19 @@ double FunctionParser::value (const Point &p, template void FunctionParser::vector_value (const Point &p, - Vector &values) const + Vector &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::vector_value (const Point &p, // have already been thrown if(dim != n_vars) vars[dim] = this->get_time(); - + for(unsigned int component = 0; component < this->n_components; ++component) { values(component) = fp[component].Eval((double*)vars); AssertThrow(fp[component].EvalError() == 0, ExcMessage(fp[component].ErrorMsg())); - } + } } #else @@ -292,7 +291,7 @@ FunctionParser::initialize(const std::string &, template double FunctionParser::value ( - const Point &, unsigned int) const + const Point &, unsigned int) const { Assert(false, ExcDisabled("parser")); return 0.; @@ -301,7 +300,7 @@ double FunctionParser::value ( template void FunctionParser::vector_value ( - const Point&, Vector&) const + const Point&, Vector&) const { Assert(false, ExcDisabled("parser")); } diff --git a/deal.II/source/base/symmetric_tensor.cc b/deal.II/source/base/symmetric_tensor.cc index c9b009af88..30ef5f61bb 100644 --- a/deal.II/source/base/symmetric_tensor.cc +++ b/deal.II/source/base/symmetric_tensor.cc @@ -18,10 +18,10 @@ DEAL_II_NAMESPACE_OPEN template <> -SymmetricTensor<4,3> -invert (const SymmetricTensor<4,3> &t) +SymmetricTensor<4,3,double> +invert<3,double> (const SymmetricTensor<4,3,double> &t) { - SymmetricTensor<4,3> tmp = t; + SymmetricTensor<4,3,double> tmp = t; // this function follows the exact same // scheme as the 2d case, except that @@ -45,7 +45,7 @@ invert (const SymmetricTensor<4,3> &t) for (unsigned int i=0; i &t) // too small Assert(max > 1.e-16*typical_diagonal_element, ExcMessage("This tensor seems to be noninvertible")); - + // row interchange if (r>j) { @@ -120,7 +120,7 @@ invert (const SymmetricTensor<4,3> &t) for (unsigned int i=3; i<6; ++i) for (unsigned int j=3; j<6; ++j) tmp.data[i][j] /= 4; - + return tmp; } diff --git a/deal.II/source/base/tensor.cc b/deal.II/source/base/tensor.cc index 16cbecd6ef..5ade4a8554 100644 --- a/deal.II/source/base/tensor.cc +++ b/deal.II/source/base/tensor.cc @@ -12,24 +12,24 @@ //--------------------------------------------------------------------------- -#include +#include #include -#include +#include DEAL_II_NAMESPACE_OPEN // storage for static variables -template -const unsigned int Tensor<1,dim>::dimension; +template +const unsigned int Tensor<1,dim,Number>::dimension; -template -const unsigned int Tensor::dimension; +template +const unsigned int Tensor::dimension; -template +template void -Tensor<1,dim>::unroll (Vector &result) const +Tensor<1,dim,Number>::unroll (Vector &result) const { Assert (result.size()==dim, ExcDimensionMismatch(dim, result.size())); @@ -40,12 +40,12 @@ Tensor<1,dim>::unroll (Vector &result) const -template +template void -Tensor::unroll (Vector &result) const +Tensor::unroll (Vector &result) const { - Assert(result.size()==std::pow(static_cast(dim),rank_), - ExcDimensionMismatch(static_cast(std::pow(static_cast(dim),rank_)), + Assert(result.size()==std::pow(static_cast(dim),rank_), + ExcDimensionMismatch(static_cast(std::pow(static_cast(dim),rank_)), result.size())); unsigned index = 0; @@ -54,40 +54,40 @@ Tensor::unroll (Vector &result) const -template +template void -Tensor::unroll_recursion (Vector &result, - unsigned int &index) const +Tensor::unroll_recursion (Vector &result, + unsigned int &index) const { for (unsigned i=0; i +template void -Tensor<1,dim>::unroll_recursion (Vector &result, - unsigned int &index) const +Tensor<1,dim,Number>::unroll_recursion (Vector &result, + unsigned int &index) const { for (unsigned i=0; i; -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 diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc index fc2fa9901f..9b09cc5103 100644 --- a/deal.II/source/fe/fe_tools.cc +++ b/deal.II/source/fe/fe_tools.cc @@ -83,6 +83,27 @@ namespace FETools { return new FE_Q<3>(quad); } + + + // Specializations for FE_DGQArbitraryNodes. + template <> + FiniteElement<1, 1>* + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<1>(quad); + } + template <> + FiniteElement<2, 2>* + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<2>(quad); + } + template <> + FiniteElement<3, 3>* + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<3>(quad); + } } namespace diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index f612f60f31..ad814bd73b 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -1118,7 +1118,7 @@ namespace internal // FEFaceValues in 1d Assert (false, ExcNotImplemented()); cross_product (boundary_forms[point], - (face_no == 0 ? 1 : -1) * cell_normal); + (face_no == 0 ? 1. : -1.) * cell_normal); } else if ( (dim==2) && (spacedim==3) ) cross_product (boundary_forms[point], diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc index c3175e6d4c..3b9e8959f5 100644 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@ -1607,7 +1607,7 @@ namespace hp "this class does not currently support this.")); create_active_fe_table (); - + tria_listeners.push_back (tria.signals.pre_refinement .connect (std_cxx1x::bind (&DoFHandler::pre_refinement_action, @@ -1631,7 +1631,7 @@ namespace hp for (unsigned int i=0; i -#include -#include #include DEAL_II_NAMESPACE_OPEN