From 8615dc74cc0c2cbe77b850e18d3b22041aa3dda9 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 26 Oct 1999 11:19:00 +0000 Subject: [PATCH] Merge unified_function_and_vectorfunction (mergepoint 1) into the main branch. git-svn-id: https://svn.dealii.org/trunk@1793 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/Todo | 3 + deal.II/base/include/base/function.h | 376 ++++++++--- deal.II/base/include/base/tensor_function.h | 172 +---- deal.II/base/include/base/tensorindex.h | 121 ++-- deal.II/base/source/function.cc | 211 +++++- deal.II/base/source/tensor_function.cc | 150 +---- .../Attic/examples/convergence/convergence.cc | 27 +- .../error-estimation/error-estimation.cc | 68 +- .../Attic/examples/multigrid/multigrid.cc | 33 +- .../fixed-point-iteration/nonlinear.cc | 2 +- .../deal.II/Attic/examples/poisson/problem.cc | 46 +- deal.II/deal.II/include/grid/intergrid_map.h | 2 +- deal.II/deal.II/include/grid/tria.h | 4 +- deal.II/deal.II/include/numerics/matrices.h | 20 +- deal.II/deal.II/include/numerics/vectors.h | 103 ++- .../source/numerics/error_estimator.cc | 11 +- deal.II/deal.II/source/numerics/matrices.cc | 236 +++++-- deal.II/deal.II/source/numerics/vectors.cc | 598 ++++++------------ tests/big-tests/convergence/convergence.cc | 27 +- .../error-estimation/error-estimation.cc | 68 +- tests/big-tests/multigrid/multigrid.cc | 33 +- .../fixed-point-iteration/nonlinear.cc | 2 +- tests/big-tests/poisson/problem.cc | 46 +- 23 files changed, 1272 insertions(+), 1087 deletions(-) diff --git a/deal.II/base/Todo b/deal.II/base/Todo index f6c223c79e..3f13b27f5a 100644 --- a/deal.II/base/Todo +++ b/deal.II/base/Todo @@ -26,3 +26,6 @@ Fill in docs for the timer class. Hopefully finally find a way to let it measure times larger than half an hour. +Review the TensorIndex class. Better documentation, remove general + template. Move constructors to the back of the file, rather than + inline in the classes. Find out whether it is really needed. diff --git a/deal.II/base/include/base/function.h b/deal.II/base/include/base/function.h index 48ddc05c26..11bb0ccf3f 100644 --- a/deal.II/base/include/base/function.h +++ b/deal.II/base/include/base/function.h @@ -6,49 +6,105 @@ /*---------------------------- function.h ---------------------------*/ +#include #include -#include -#include #include - +#include +#include /** - * This class is a model for a continuous function. It returns the value - * of the function at a given point through the #operator ()# member function, - * which is virtual. It also has a function to return a whole list of function - * values at different points to reduce the overhead of the virtual function - * calls; this function is preset to successively call the function returning - * one value at a time. + * This class is a model for a general function. It serves the purpose + * of representing scalar and vector valued functions. To this end, we + * consider scalar functions as a special case of vector valued + * functions, in the former case only having a single component return + * vector. Since handling with vectors is comparatively expensive, + * functions are provided which only ask for a single component of the + * function, which is what you will usually need in case you know that + * your function is scalar-valued. * - * There are other functions return the gradient of the function at one or - * several points. You only have to overload those functions you need; the - * functions returning several values at a time will call those returning - * only one value, while those ones will throw an exception when called but - * not overloaded. + * Access to function objects therefore is through the following + * methods: + * \begin{verbatim} + * // access to one component at one point + * double value (const Point &p, + * const unsigned int component = 0) const; * - * Unless only called a very small number of times, you should overload - * both those functions returning only one value as well as those returning - * a whole array, since the cost of evaluation of a point value is often - * less than the virtual function call itself. + * // return all components at one point + * void vector_value (const Point &p, + * Vector &value) const; + * \end{verbatim} * + * For more efficiency, there are other functions returning one or all + * components at a list of points at once: + * \begin{verbatim} + * // access to one component at several points + * void value_list (const vector > &point_list, + * vector &value_list, + * const unsigned int component = 0) const; * - * Support for time dependant functions can be found in the base - * class #FunctionTime#. - - * @author Wolfgang Bangerth, 1998, 1999 + * // return all components at several points + * void vector_value_list (const vector > &point_list, + * vector > &value_list) const; + * \end{verbatim} + * + * Furthermore, there are functions returning the gradient of the + * function at one or several points. + * + * You will usually only overload those functions you need; the + * functions returning several values at a time (#value_list#, + * #vector_value_list#, and gradient analoga) will call those + * returning only one value (#value#, #vector_value#, and gradient + * analoga), while those ones will throw an exception when called but + * not overloaded. + * + * Note however, that the functions returning all components of the + * function at one or several points (i.e. #vector_value#, + * #vector_value_list#), will not call the function returning one + * component at one point repeatedly, once for each point and + * component. The reason is efficiency: this would amount to too many + * virtual function calls. If you have vector-valued functions, you + * should therefore also provide overloads of the virtual functions + * for all components at a time. + * + * Also note, that unless only called a very small number of times, + * you should overload all sets of functions (returning only one + * value, as well as those returning a whole array), since the cost of + * evaluation of a point value is often less than the virtual function + * call itself. + * + * + * Support for time dependant functions can be found in the base + * class #FunctionTime#. + * + * {\bf Note}: if the functions you are dealing with have sizes which + * are a priori known (for example, #dim# elements), you might + * consider using the #TensorFunction# class instead. + * + * @author Wolfgang Bangerth, 1998, 1999 */ template -class Function : public FunctionTime +class Function : public FunctionTime, + public Subscriptor { public: + /** + * Number of vector components. + */ + const unsigned int n_components; + /** - * Constructor. May take an initial vakue - * for the time variable, which defaults - * to zero. + * Constructor. May take an + * initial value for the number + * of components (which defaults + * to one, i.e. a scalar + * function), and the time + * variable, which defaults to + * zero. */ - Function (const double initial_time = 0.0); + Function (const unsigned int n_components = 1, + const double initial_time = 0.0); /** * Virtual destructor; absolutely @@ -57,38 +113,111 @@ class Function : public FunctionTime virtual ~Function (); /** - * Return the value of the function - * at the given point. + * Return the value of the + * function at the given + * point. Unless there is only + * one component (i.e. the + * function is scalar), you + * should state the component you + * want to have evaluated; it + * defaults to zero, i.e. the + * first component. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; /** - * Set #values# to the point values - * of the function at the #points#. - * It is assumed that #values# - * already has the right size, i.e. - * the same size as the #points# - * array. + * Return all components of a + * vector-valued function at a + * given point. + * + * Be default, this function + * repeatedly calls the other + * #operator()# for each + * component separately, to fill + * the output array. + * + * #values# shall have the right + * size beforehand, + * i.e. #n_components#. + */ + virtual void vector_value (const Point &p, + Vector &values) const; + + /** + * Set #values# to the point + * values of the specified + * component of the function at + * the #points#. It is assumed + * that #values# already has the + * right size, i.e. the same + * size as the #points# array. */ virtual void value_list (const vector > &points, - vector &values) const; + vector &values, + const unsigned int component = 0) const; /** - * Return the gradient of the function - * at the given point. + * Set #values# to the point + * values of the function at the + * #points#. It is assumed that + * #values# already has the right + * size, i.e. the same size as + * the #points# array, and that + * all elements be vectors with + * the same number of components + * as this function has. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual void vector_value_list (const vector > &points, + vector > &values) const; + + /** + * Return the gradient of the + * specified component of the + * function at the given point. + */ + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + /** + * Return the gradient of the + * specified component of the + * function at the given point, + * for all components. + */ + virtual void vector_gradient (const Point &p, + vector > &gradients) const; + + /** + * Set #gradients# to the + * gradients of the specified + * component of the function at + * the #points#. It is assumed + * that #gradients# already has the + * right size, i.e. the same + * size as the #points# array. + */ + virtual void gradient_list (const vector > &points, + vector > &gradients, + const unsigned int component = 0) const; + /** * Set #gradients# to the gradients of - * the function at the #points#. - * It is assumed that #values# + * the function at the #points#, + * for all components. + * It is assumed that #gradients# * already has the right size, i.e. * the same size as the #points# array. + * + * The outer loop over + * #gradients# is over the points + * in the list, the inner loop + * over the different components + * of the function. */ - virtual void gradient_list (const vector > &points, - vector > &gradients) const; - + virtual void vector_gradient_list (const vector > &points, + vector > > &gradients) const; + /** * Exception */ @@ -106,52 +235,118 @@ class Function : public FunctionTime /** - * Provide a function which always returns zero. Obviously, also the derivates - * of this function are zero. + * Provide a function which always returns zero. Obviously, also the + * derivates of this function are zero. Also, it returns zero on all + * components in case the function is not a scalar one, which can be + * obtained by passing the constructor the appropriate number of + * components. + * + * This function is of use when you want to implement homogeneous boundary + * conditions, or zero initial conditions. * - * This function is of use when you want to implement homogeneous boundary - * conditions. + * @author Wolfgang Bangerth, 1998, 1999 */ template class ZeroFunction : public Function { public: + /** + * Constructor. The number of + * 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. + * at the given point for one + * component. + */ + virtual double value (const Point &p, + const unsigned int component) const; + + /** + * Return the value of the function + * at the given point for all + * components. */ - virtual double operator () (const Point &p) const; + virtual void vector_value (const Point &p, + Vector &return_value) const; /** * Set #values# to the point values - * of the function at the #points#. + * of the function at the #points#, + * for the given component. * It is assumed that #values# * already has the right size, i.e. * the same size as the #points# * array. */ virtual void value_list (const vector > &points, - vector &values) const; + vector &values, + const unsigned int component = 0) const; + /** + * Set #values# to the point values + * of the function at the #points#, + * for all components. + * It is assumed that #values# + * already has the right size, i.e. + * the same size as the #points# + * array. + */ + virtual void vector_value_list (const vector > &points, + vector > &values) const; + /** * Return the gradient of the function - * at the given point. + * at the given point, for the + * given component. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + /** + * Return the gradient of the + * specified component of the + * function at the given point, + * for all components. + */ + virtual void vector_gradient (const Point &p, + vector > &gradients) const; + /** * Set #gradients# to the gradients of - * the function at the #points#. + * the function at the #points#, + * for the given component. * It is assumed that #values# * already has the right size, i.e. * the same size as the #points# array. */ virtual void gradient_list (const vector > &points, - vector > &gradients) const; + 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# + * already has the right size, i.e. + * the same size as the #points# array. + * + * The outer loop over + * #gradients# is over the points + * in the list, the inner loop + * over the different components + * of the function. + */ + virtual void vector_gradient_list (const vector > &points, + vector > > &gradients) const; }; @@ -159,46 +354,83 @@ class ZeroFunction : public Function { /** - * Provide a function which always returns a constant value, which is delivered - * upon construction. Obviously, the derivates of this function are zero, which - * is why we derive this class from #ZeroFunction#: we then only have to - * overload th value functions, not all the derivatives. In some way, it would - * be more obvious to do the derivation in the opposite direction, i.e. let - * #ZeroFunction# be a more specialized version of #ConstantFunction#; however, - * this would be more inefficient, since we could not make use of the fact that - * the function value of the #ZeroFunction# is known at compile time and need - * not be looked up somewhere in memory. + * Provide a function which always returns a constant value, which is + * delivered upon construction. Obviously, the derivates of this + * function are zero, which is why we derive this class from + * #ZeroFunction#: we then only have to overload th value functions, + * not all the derivatives. In some way, it would be more obvious to + * do the derivation in the opposite direction, i.e. let + * #ZeroFunction# be a more specialized version of #ConstantFunction#; + * however, this would be more inefficient, since we could not make + * use of the fact that the function value of the #ZeroFunction# is + * known at compile time and need not be looked up somewhere in + * memory. + * + * You can pass to the constructor an integer denoting the number of + * components this function shall have. It defaults to one. If it is + * greater than one, then the function will return the constant value + * in all its components, which might not be overly useful a feature + * in most cases, however. + * + * @author Wolfgang Bangerth, 1998, 1999 */ template class ConstantFunction : public ZeroFunction { public: /** * Constructor; takes the constant function - * value as an argument. + * value as an argument. The number of + * components is preset to one. */ - ConstantFunction (const double value); + ConstantFunction (const double value, + const unsigned int n_components = 1); /** * Virtual destructor; absolutely * necessary in this case. */ virtual ~ConstantFunction (); + /** * Return the value of the function - * at the given point. + * at the given point for one + * component. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + + /** + * Return the value of the function + * at the given point for all + * components. + */ + virtual void vector_value (const Point &p, + Vector &return_value) const; /** * Set #values# to the point values - * of the function at the #points#. + * of the function at the #points#, + * for the given component. * It is assumed that #values# * already has the right size, i.e. * the same size as the #points# * array. */ virtual void value_list (const vector > &points, - vector &values) const; + vector &values, + const unsigned int component = 0) const; + + /** + * Set #values# to the point values + * of the function at the #points#, + * for all components. + * It is assumed that #values# + * already has the right size, i.e. + * the same size as the #points# + * array. + */ + virtual void vector_value_list (const vector > &points, + vector > &values) const; protected: /** diff --git a/deal.II/base/include/base/tensor_function.h b/deal.II/base/include/base/tensor_function.h index e79c73d21b..6d0cacbc78 100644 --- a/deal.II/base/include/base/tensor_function.h +++ b/deal.II/base/include/base/tensor_function.h @@ -13,154 +13,9 @@ #include #include #include -#include #include -template class Vector; -template class VectorFunction; -template class TensorFunction; -/** - * Base class for multi-valued functions. - * While #TensorFunction# provides a highly structured class for multi-valued - * functions, #VectorFunction# is on a lower level. The results are #Vectors# of - * values without further structure. The dimension of the result is determined at - * execution time. - * @author Guido Kanschat, 1999 - */ -template -class VectorFunction : public FunctionTime, - public Subscriptor -{ - public: - /** - * Number of vector components. - */ - const unsigned int n_components; - - /** - * Constructor. May take an initial vakue - * for the time variable, which defaults - * to zero. - */ - VectorFunction (const unsigned int n_components, - const double initial_time = 0.0); - - /** - * Virtual destructor; absolutely - * necessary in this case. - */ - virtual ~VectorFunction (); - - /** - * Set #values# to the point values - * of the function at points #p#. - * It is assumed that #values# - * already has the right size, i.e. - * the same size as the #n_components# - * array. - * - * Usually only #value_list# is called, - * e.g. by #FEValues#. So, to avoid - * multiple calling of this virtual function - * by #value_list#, implement the vectorfunction - * directly in #value_list# of the derived - * class. - */ - virtual void value (const Point &p, - Vector &values) const; - - /** - * Set #values# to the point values - * of the function at the #points#. - * It is assumed that #values# - * already has the right size, i.e. - * the same size as the #points# - * array. - * - * This function uses multiple calling - * of the virtual function #value# (see there). - * If possible, overload this function. - */ - virtual void value_list (const vector > &points, - vector > &values) const; - - /** - * Set #gradients# to the gradients of - * the function at the #points#. - * It is assumed that #values# - * already has the right size, i.e. - * the same size as the #points# array. - */ - virtual void gradient_list (const vector > &points, - vector > > &gradients) const; - - /** - * Access #VectorFunction# as a #Function#. - * This class allows to store a reference to a - * #VectorFunction# and an #index#. Later on, it - * can be used as a normal single valued #Function#. - */ - class Extractor : public Function - { - public: - /** - * Constructor. - * The arguments are the #VectorFunction# to be - * accessed and the component index. - */ - Extractor(const VectorFunction& f, unsigned int index); - - /** - * Compute function value. - */ - virtual double operator() (const Point& p) const; - - /** - * Compute several values. - */ - virtual void value_list (const vector > &points, - vector &values) const; - - - /** - * Compute derivative. - */ - virtual Tensor<1,dim> gradient (const Point& p) const; - - /** - * Compute several derivatives. - */ - virtual void gradient_list (const vector > &points, - vector > &gradients) const; - - private: - /** - * Pointer to the #VectorFunction#. - */ - const SmartPointer > vectorfunction; - - /** - * Index in #VectorFunction#. - */ - const unsigned int index; - }; - - - /** - * Exception - */ - DeclException0 (ExcPureFunctionCalled); - /** - * Exception - */ - DeclException2 (ExcVectorHasWrongSize, - int, int, - << "The vector has size " << arg1 << " but should have " - << arg2 << " elements."); - -}; - /** @@ -184,7 +39,8 @@ class VectorFunction : public FunctionTime, * @author Guido Kanschat, 1999 */ template -class TensorFunction : public VectorFunction +class TensorFunction : public FunctionTime, + public Subscriptor { public: /** @@ -204,7 +60,7 @@ class TensorFunction : public VectorFunction * Return the value of the function * at the given point. */ - virtual Tensor operator () (const Point &p) const; + virtual Tensor value (const Point &p) const; /** * Set #values# to the point values @@ -218,8 +74,8 @@ class TensorFunction : public VectorFunction vector > &values) const; /** - * Return the gradient of the function - * at the given point. + * Return the gradient of the + * function at the given point. */ virtual Tensor gradient (const Point &p) const; @@ -233,24 +89,6 @@ class TensorFunction : public VectorFunction virtual void gradient_list (const vector > &points, vector > &gradients) const; - /** - * See #VectorFunction#. - */ - virtual void value (const Point &points, - Vector &values) const; - - /** - * See #VectorFunction#. - */ - virtual void value_list (const vector > &points, - vector > &values) const; - - /** - * See #VectorFunction#. - */ - virtual void gradient_list (const vector > &points, - vector > > &gradients) const; - /** * Exception */ diff --git a/deal.II/base/include/base/tensorindex.h b/deal.II/base/include/base/tensorindex.h index bf4e0b598a..b9d0ddb81b 100644 --- a/deal.II/base/include/base/tensorindex.h +++ b/deal.II/base/include/base/tensorindex.h @@ -6,8 +6,8 @@ #include + /** - * * Rank-independent access to elements of #Tensor#. A little class * template remembering #rank# integers. * @@ -16,49 +16,52 @@ * or contact the developer. * * @author Guido Kanschat, 1999 - * */ template class TensorIndex { -private: - /** - * Field of indices. - */ - unsigned int index[rank]; -public: - /** - * Constructor taking #rank# indices. - */ - TensorIndex(...); - - /** - * Access operator returning index - * in #n#th component - */ - unsigned int operator () (unsigned int n) const; - + private: + /** + * Field of indices. + */ + unsigned int index[rank]; + public: + /** + * Constructor taking #rank# indices. + */ + TensorIndex(...); + + /** + * Access operator returning index + * in #n#th component + */ + unsigned int operator () (const unsigned int n) const; + /** * Exception. */ DeclException1(ExcRank, int, - << "Index " << arg1 << " higher than maximum " << rank-1); + << "Index " << arg1 << " higher than maximum " << rank-1); }; + template<> class TensorIndex<4> { -private: - /** - * Field of indices. - */ - unsigned int index[4]; -public: - /** - * Constructor taking #rank# indices. - */ - TensorIndex(unsigned int i0, unsigned int i1, unsigned int i2, unsigned int i3) + private: + /** + * Field of indices. + */ + unsigned int index[4]; + public: + /** + * Constructor taking #rank# indices. + */ + TensorIndex (const unsigned int i0, + const unsigned int i1, + const unsigned int i2, + const unsigned int i3) { index[0] = i0; index[1] = i1; @@ -67,33 +70,41 @@ public: } - /** - * Access operator returning index - * in #n#th component - */ - unsigned int operator () (unsigned int n) const + /** + * Access operator returning index + * in #n#th component + */ + unsigned int operator () (const unsigned int n) const { Assert(n<4, ExcRank(n)); return index[n]; } - DeclException1(ExcRank, unsigned int, - << "Index " << arg1 << " higher than maximum 3"); + + /** + * Exception + */ + DeclException1(ExcRank, unsigned int, + << "Index " << arg1 << " higher than maximum 3"); }; + + template<> class TensorIndex<3> { -private: - /** - * Field of indices. - */ + private: + /** + * Field of indices. + */ unsigned int index[3]; public: /** * Constructor taking #rank# indices. */ - TensorIndex(unsigned int i0, unsigned int i1, unsigned int i2) + TensorIndex(const unsigned int i0, + const unsigned int i1, + const unsigned int i2) { index[0] = i0; index[1] = i1; @@ -106,17 +117,22 @@ private: * in #n#th component * */ - unsigned int operator () (unsigned int n) const + unsigned int operator () (const unsigned int n) const { Assert(n<3, ExcRank(n)); return index[n]; } + /** + * Exception + */ DeclException1(ExcRank, unsigned int, << "Index " << arg1 << " higher than maximum 2"); }; + + template<> class TensorIndex<2> { @@ -129,7 +145,8 @@ class TensorIndex<2> /** * Constructor taking #rank# indices. */ - TensorIndex(unsigned int i0, unsigned int i1) + TensorIndex(const unsigned int i0, + const unsigned int i1) { index[0] = i0; index[1] = i1; @@ -140,16 +157,22 @@ class TensorIndex<2> * Access operator returning index * in #n#th component */ - unsigned int operator () (unsigned int n) const + unsigned int operator () (const unsigned int n) const { Assert(n<2, ExcRank(n)); return index[n]; } + + /** + * Exception + */ DeclException1(ExcRank, unsigned int, << "Index " << arg1 << " higher than maximum 1"); }; + + template<> class TensorIndex<1> { @@ -162,7 +185,7 @@ class TensorIndex<1> /** * Constructor taking #rank# indices. */ - TensorIndex(unsigned int i0) + TensorIndex(const unsigned int i0) { index[0] = i0; } @@ -172,12 +195,16 @@ class TensorIndex<1> * Access operator returning index * in #n#th component */ - unsigned int operator () (unsigned int n) const + unsigned int operator () (const unsigned int n) const { Assert(n<1, ExcRank(n)); return index[n]; } + + /** + * Exception + */ DeclException1(ExcRank, unsigned int, << "Index " << arg1 << " higher than maximum 0"); }; diff --git a/deal.II/base/source/function.cc b/deal.II/base/source/function.cc index d2905c2a5b..9d0da2c5c2 100644 --- a/deal.II/base/source/function.cc +++ b/deal.II/base/source/function.cc @@ -2,12 +2,16 @@ #include +#include +#include #include template -Function::Function (const double initial_time) : - FunctionTime(initial_time) +Function::Function (const unsigned int n_components, + const double initial_time) : + FunctionTime(initial_time), + n_components(n_components) {}; @@ -19,46 +23,116 @@ Function::~Function () template -double Function::operator () (const Point &) const { +double Function::value (const Point &, + const unsigned int) const +{ Assert (false, ExcPureFunctionCalled()); return 0; }; +template +void Function::vector_value (const Point &, + Vector &) const +{ + Assert (false, ExcPureFunctionCalled()); +}; + + + template void Function::value_list (const vector > &points, - vector &values) const { + vector &values, + const unsigned int component) const +{ + // check whether component is in + // the valid range is up to the + // derived class Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; ioperator() (points[i]); + values[i] = this->value (points[i], component); +}; + + + +template +void Function::vector_value_list (const vector > &points, + vector > &values) const +{ + // check whether component is in + // the valid range is up to the + // derived class + Assert (values.size() == points.size(), + ExcVectorHasWrongSize(values.size(), points.size())); + + for (unsigned int i=0; ivector_value (points[i], values[i]); }; template -Tensor<1,dim> Function::gradient (const Point &) const { +Tensor<1,dim> Function::gradient (const Point &, + const unsigned int) const +{ Assert (false, ExcPureFunctionCalled()); return Point(); }; +template +void Function::vector_gradient (const Point &, + vector > &) const +{ + Assert (false, ExcPureFunctionCalled()); +}; + + + template void Function::gradient_list (const vector > &points, - vector > &gradients) const { + vector > &gradients, + const unsigned int component) const +{ Assert (gradients.size() == points.size(), ExcVectorHasWrongSize(gradients.size(), points.size())); for (unsigned int i=0; i +void Function::vector_gradient_list (const vector > &points, + vector > > &gradients) const +{ + Assert (gradients.size() == points.size(), + ExcVectorHasWrongSize(gradients.size(), points.size())); + + for (unsigned int i=0; i +ZeroFunction::ZeroFunction (const unsigned int n_components) : + Function (n_components) +{}; + template @@ -67,15 +141,30 @@ ZeroFunction::~ZeroFunction () {}; template -double ZeroFunction::operator () (const Point &) const { +double ZeroFunction::value (const Point &, + const unsigned int) const +{ return 0.; }; +template +void ZeroFunction::vector_value (const Point &, + Vector &return_value) const +{ + Assert (return_value.size() == n_components, + ExcVectorHasWrongSize (return_value.size(), n_components)); + + fill_n (return_value.begin(), n_components, 0.0); +}; + + + template void ZeroFunction::value_list (const vector > &points, - vector &values) const { + vector &values, + const unsigned int /*component*/) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); @@ -85,27 +174,82 @@ void ZeroFunction::value_list (const vector > &points, template -Tensor<1,dim> ZeroFunction::gradient (const Point &) const { +void ZeroFunction::vector_value_list (const vector > &points, + vector > &values) const +{ + Assert (values.size() == points.size(), + ExcVectorHasWrongSize(values.size(), points.size())); + + for (unsigned int i=0; i +Tensor<1,dim> ZeroFunction::gradient (const Point &, + const unsigned int) const +{ return Tensor<1,dim>(); }; +template +void ZeroFunction::vector_gradient (const Point &, + vector > &gradients) const +{ + Assert (gradients.size() == n_components, + ExcVectorHasWrongSize(gradients.size(), n_components)); + + for (unsigned int c=0; c void ZeroFunction::gradient_list (const vector > &points, - vector > &gradients) const { + vector > &gradients, + const unsigned int /*component*/) const +{ Assert (gradients.size() == points.size(), ExcVectorHasWrongSize(gradients.size(), points.size())); - gradients.clear (); + for (unsigned int i=0; i +void ZeroFunction::vector_gradient_list (const vector > &points, + vector > > &gradients) const +{ + Assert (gradients.size() == points.size(), + ExcVectorHasWrongSize(gradients.size(), points.size())); + for (unsigned int i=0; i -ConstantFunction::ConstantFunction (const double value) : - function_value(value) {}; +ConstantFunction::ConstantFunction (const double value, + const unsigned int n_components) : + ZeroFunction (n_components), + function_value (value) {}; template @@ -114,15 +258,30 @@ ConstantFunction::~ConstantFunction () {}; template -double ConstantFunction::operator () (const Point &) const { +double ConstantFunction::value (const Point &, + const unsigned int) const +{ return function_value; }; +template +void ConstantFunction::vector_value (const Point &, + Vector &return_value) const +{ + Assert (return_value.size() == n_components, + ExcVectorHasWrongSize (return_value.size(), n_components)); + + fill_n (return_value.begin(), n_components, function_value); +}; + + + template void ConstantFunction::value_list (const vector > &points, - vector &values) const { + vector &values, + const unsigned int /*component*/) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); @@ -131,6 +290,24 @@ void ConstantFunction::value_list (const vector > &points, +template +void ConstantFunction::vector_value_list (const vector > &points, + vector > &values) const +{ + Assert (values.size() == points.size(), + ExcVectorHasWrongSize(values.size(), points.size())); + + for (unsigned int i=0; i; diff --git a/deal.II/base/source/tensor_function.cc b/deal.II/base/source/tensor_function.cc index 6b7a11ac25..34801c9c77 100644 --- a/deal.II/base/source/tensor_function.cc +++ b/deal.II/base/source/tensor_function.cc @@ -7,104 +7,7 @@ #include #include -template -VectorFunction::VectorFunction(unsigned n_components, const double initial_time) - : - FunctionTime(initial_time), - n_components(n_components) -{} - - -template -VectorFunction::~VectorFunction() -{} - -/* -template double -VectorFunction::operator () (const Point &, unsigned) const - -{ - Assert (false, ExcPureFunctionCalled()); - return 0.; -} -*/ - -template -void -VectorFunction::value (const Point &, Vector &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - -template -void -VectorFunction::value_list (const vector > &ps, - vector > &us) const -{ - for (unsigned int i=0 ; i -void -VectorFunction::gradient_list (const vector > &, - vector > > &) const -{ - Assert (false, ExcPureFunctionCalled()); -} -template -VectorFunction::Extractor::Extractor(const VectorFunction& f, - unsigned int index) - : - vectorfunction(f), - index(index) -{} - -template -double -VectorFunction::Extractor::operator() (const Point& p) const -{ - Vector v(vectorfunction->n_components); - vectorfunction->value(p,v); - return v(index); -} - - -template -Tensor<1,dim> -VectorFunction::Extractor::gradient (const Point&) const -{ - Assert(false, ExcNotImplemented()); - return Tensor<1,dim>(); -} - -template -void -VectorFunction::Extractor::value_list (const vector > &points, - vector &values) const -{ - vector > v(values.size(), - Vector(vectorfunction->n_components)); - vectorfunction->value_list(p,v); - for (unsigned int i=0 ; i -void -VectorFunction::Extractor::gradient_list (const vector > &points, - vector > &gradients) const -{ - vector > > v(values.size(), - vector >(vectorfunction->n_components)); - vectorfunction->value_list(p,v); - for (unsigned int i=0 ; i::Extractor::gradient_list (const vector > &points template TensorFunction::TensorFunction (const double initial_time) : - VectorFunction(pow(dim,rank), initial_time) + FunctionTime (initial_time) {}; @@ -124,41 +27,30 @@ TensorFunction::~TensorFunction () -// template -// double -// TensorFunction::operator () (TensorIndex i, -// const Point &) const -// { -// int k=i(0); -// k++; - -// Assert (false, ExcPureFunctionCalled()); -// return 0; -// }; - - template Tensor -TensorFunction::operator() (const Point &) const +TensorFunction::value (const Point &) const { Assert (false, ExcPureFunctionCalled()); return Tensor(); }; + template void TensorFunction::value_list (const vector > &points, - vector > &values) const + vector > &values) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; ioperator() (points[i]); + values[i] = this->value (points[i]); }; + template Tensor TensorFunction::gradient (const Point &) const @@ -182,36 +74,6 @@ TensorFunction::gradient_list (const vector > &points, }; -template void -TensorFunction::value (const Point &p, - Vector &erg) const -{ - Tensor h = operator()(p); - h.unroll(erg); -} - - -template void -TensorFunction::value_list (const vector > & points, - vector > & values) const -{ - Assert (values.size() == points.size(), - ExcVectorHasWrongSize(values.size(), points.size())); - - for (unsigned int i=0; i void -TensorFunction::gradient_list (const vector > &, - vector > > &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - template class TensorFunction<1,1>; template class TensorFunction<2,1>; diff --git a/deal.II/deal.II/Attic/examples/convergence/convergence.cc b/deal.II/deal.II/Attic/examples/convergence/convergence.cc index 1a5ad43c5c..afe9ac1f4b 100644 --- a/deal.II/deal.II/Attic/examples/convergence/convergence.cc +++ b/deal.II/deal.II/Attic/examples/convergence/convergence.cc @@ -94,7 +94,8 @@ class RHSPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; @@ -106,19 +107,22 @@ class Solution : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; /** * Return the gradient of the function * at the given point. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; template <> -double RHSPoly<2>::operator () (const Point<2> &p) const { +double RHSPoly<2>::value (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -128,7 +132,8 @@ double RHSPoly<2>::operator () (const Point<2> &p) const { template <> -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::value (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -137,7 +142,8 @@ double Solution<2>::operator () (const Point<2> &p) const { template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::gradient (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -162,7 +168,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, fe_values.shape_grad(j,point)) * fe_values.JxW(point); rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -376,8 +382,11 @@ int PoissonProblem::run (const unsigned int level) { if (dof->n_dofs()<=5000) { - Vector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; - Vector h1_seminorm_error_per_dof, h1_error_per_dof; + Vector l1_error_per_dof(dof->n_dofs()); + Vector l2_error_per_dof(dof->n_dofs()); + Vector linfty_error_per_dof(dof->n_dofs()); + Vector h1_seminorm_error_per_dof(dof->n_dofs()); + Vector h1_error_per_dof(dof->n_dofs()); dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc index 5158c8c901..aebd9cf999 100644 --- a/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc +++ b/deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc @@ -114,25 +114,32 @@ class Solution { class GaussShape : public Function { public: - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; class Singular : public Function { public: - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; class Kink : public Function { public: class Coefficient : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; }; @@ -150,7 +157,8 @@ class RHS { */ class GaussShape : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; /** @@ -160,7 +168,8 @@ class RHS { */ class Singular : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; /** @@ -171,7 +180,8 @@ class RHS { */ class Kink : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; }; @@ -179,13 +189,15 @@ class RHS { template <> -double Solution<2>::GaussShape::operator () (const Point<2> &p) const { +double Solution<2>::GaussShape::value (const Point<2> &p, + const unsigned int) const { return p(0)*p(1)*exp(-40*p.square()); }; template <> -Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p, + const unsigned int) const { return Point<2> ((1-80.*p(0)*p(0))*p(1)*exp(-40*p.square()), (1-80.*p(1)*p(1))*p(0)*exp(-40*p.square())); }; @@ -193,13 +205,15 @@ Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { template <> -double Solution<2>::Singular::operator () (const Point<2> &p) const { +double Solution<2>::Singular::value (const Point<2> &p, + const unsigned int) const { return pow(p.square(), 1./3.); }; template <> -Tensor<1,2> Solution<2>::Singular::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::Singular::gradient (const Point<2> &p, + const unsigned int) const { return 2./3.*pow(p.square(), -2./3.) * p; }; @@ -213,21 +227,24 @@ inline double theta(const double x) { template <> -double Solution<2>::Kink::operator () (const Point<2> &p) const { +double Solution<2>::Kink::value (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return (1+4*theta(s))*s; }; template <> -Tensor<1,2> Solution<2>::Kink::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::Kink::gradient (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return (1+4*theta(s))*Point<2>(-2*p(0),1); }; template <> -double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { +double Solution<2>::Kink::Coefficient::value (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return 1./(1.+4.*theta(s)); }; @@ -235,19 +252,22 @@ double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { template <> -double RHS<2>::GaussShape::operator () (const Point<2> &p) const { +double RHS<2>::GaussShape::value (const Point<2> &p, + const unsigned int) const { return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); }; template <> -double RHS<2>::Singular::operator () (const Point<2> &p) const { +double RHS<2>::Singular::value (const Point<2> &p, + const unsigned int) const { return -4./9. * pow(p.square(), -2./3.); }; template <> -double RHS<2>::Kink::operator () (const Point<2> &) const { +double RHS<2>::Kink::value (const Point<2> &, + const unsigned int) const { return 2; }; @@ -266,7 +286,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, for (unsigned int point=0; point::assemble (FullMatrix &cell_matrix, fe_values.JxW(point) * c; rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -535,8 +555,8 @@ void PoissonProblem::run (ParameterHandler &prm) { cout << estimated_error_per_cell.l2_norm() << endl; estimated_error.push_back (estimated_error_per_cell.l2_norm()); - Vector l2_error_per_dof, linfty_error_per_dof; - Vector h1_error_per_dof, estimated_error_per_dof; + Vector l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs()); + Vector h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs()); Vector error_ratio (dof->n_dofs()); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc index 572a209dea..317a7282c9 100644 --- a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc +++ b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc @@ -174,7 +174,8 @@ class RHSPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; @@ -186,19 +187,24 @@ class Solution : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; /** * Return the gradient of the function * at the given point. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; template <> -double RHSPoly<2>::operator () (const Point<2> &p) const { +double RHSPoly<2>::value (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -208,7 +214,10 @@ double RHSPoly<2>::operator () (const Point<2> &p) const { template <> -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::value (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -217,7 +226,10 @@ double Solution<2>::operator () (const Point<2> &p) const { template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::gradient (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -242,7 +254,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, fe_values.shape_grad(j,point)) * fe_values.JxW(point); rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -534,8 +546,11 @@ int PoissonProblem::run (const unsigned int level) { if (dof->DoFHandler::n_dofs()<=5000) { - Vector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; - Vector h1_seminorm_error_per_dof, h1_error_per_dof; + Vector l1_error_per_dof (dof->DoFHandler::n_dofs()); + Vector l2_error_per_dof (dof->DoFHandler::n_dofs()); + Vector linfty_error_per_dof (dof->DoFHandler::n_dofs()); + Vector h1_seminorm_error_per_dof (dof->DoFHandler::n_dofs()); + Vector h1_error_per_dof (dof->DoFHandler::n_dofs()); dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc b/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc index dd4053a9e2..0034eb8e9f 100644 --- a/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc +++ b/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc @@ -33,7 +33,7 @@ template class RightHandSide : public Function { public: - double operator () (const Point &p) const + double value (const Point &p) const { double x = 80; for (unsigned int d=0; d { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const { + virtual double value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + double x = 1; for (unsigned int i=0; i &p, + Vector &values) const { + Assert (values.size()==1, ExcVectorHasWrongSize (values.size(), 1)); + + double x = 1; + + for (unsigned int i=0; i { * empty. */ virtual void value_list (const vector > &points, - vector &values) const { + vector &values, + const unsigned int component) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; i::operator() (points[i]); + values[i] = BoundaryValuesSine::value (points[i], component); }; }; @@ -49,7 +69,9 @@ class BoundaryValuesJump : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const { + virtual double value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); switch (dim) { case 1: @@ -73,7 +95,8 @@ class RHSTrigPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int) const; }; @@ -89,7 +112,8 @@ class RHSPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int) const; }; @@ -200,7 +224,10 @@ CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_ template -double RHSTrigPoly::operator () (const Point &p) const { +double RHSTrigPoly::value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double pi = 3.1415926536; switch (dim) { @@ -219,7 +246,10 @@ double RHSTrigPoly::operator () (const Point &p) const { template -double RHSPoly::operator () (const Point &p) const { +double RHSPoly::value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + double ret_val = 0; for (unsigned int i=0; i# + * Note that this class could in principle be based on the C++ #map# * data type. Instead, it uses another data format which is more effective both * in terms of computing time for access as well as with regard to memory * consumpion. diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index 3d4484f876..17b9a83cc2 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -1186,8 +1186,10 @@ struct TriaNumberCache<3> : public TriaNumberCache<2> * void main () { * Triangulation<2> tria; * // set the boundary function + * // for all boundaries with + * // boundary indicator 0 * Ball ball; - * tria.set_boundary (&ball); + * tria.set_boundary (0, &ball); * * // read some coarse grid * diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 5e5b691b85..0b84b55ab6 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -67,21 +67,25 @@ * cross coupling of shape functions belonging to different subelements. * * If the finite element for which the mass matrix is to be built - * has more than one component, the resulting matrix will not - * couple the different components. It will furthermore accept - * a single coefficient through the #Function# parameter for all + * has more than one component, the resulting matrix will not couple + * the different components. It will furthermore accept a single + * coefficient through the #Function# parameter for all * components. If you want different coefficients for the different - * parameters, you need to call the respective function accepting - * a #VectorFunction# argument. + * parameters, you need to pass a function object representing the + * respective number of components. * * \item #create_laplace_matrix#: there are two versions of this; the * one which takes the #Function# object creates * $a_{ij} = \int_\Omega a(x) \nabla\phi_i(x) \nabla\phi_j(x) dx$, * $a$ being the given function, while the other one assumes that - * $a=1$ which enables some optimzations. In fact the two versions + * $a=1$ which enables some optimizations. In fact the two versions * are in one function, the coefficient being given as a defaulted * argument, which is a pointer to a function and defaults to zero. * This function uses the #LaplaceMatrix# class. + * + * If the finite element in use presently has more than only one + * component, this function may not be overly useful and presently + * throws an error. * \end{itemize} * * All created matrices are `raw': they are not condensed, i.e. hanging @@ -320,6 +324,10 @@ class MatrixCreator * Exception */ DeclException0 (ExcInvalidFE); + /** + * Exception + */ + DeclException0 (ExcComponentMismatch); }; diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 0e798566d0..6cbb45f2ec 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -269,36 +269,20 @@ class VectorTools */ typedef map*> FunctionMap; - /** - * Data type for vector valued boundary function map. - */ - typedef map*> VectorFunctionMap; - /** * Compute the interpolation of - * #function# at the support points to - * the finite element space. + * #function# at the support + * points to the finite element + * space. It is assumed that the + * number of components of + * #function# matches that of the + * finite element used by #dof#. * * See the general documentation of this * class for further information. */ static void interpolate (const DoFHandler &dof, const Function &function, - Vector &vec); - - /** - * Compute the interpolation of - * #vectorfunction# at the support points to - * the finite element space. This is the - * analogue for vectorfunctions - * to the #interpolate# function for scalar - * functions above. - * - * See the general documentation of this - * class for further information. - */ - static void interpolate (const DoFHandler &dof, - const VectorFunction&vectorfunction, Vector &vec); /** @@ -383,6 +367,12 @@ class VectorTools * of the boundary part to be projected * on already was in the variable. * + * It is assumed that the number + * of components of the functions + * in #dirichlet_bc# matches that + * of the finite element used by + * #dof#. + * * See the general doc for more * information. */ @@ -390,15 +380,6 @@ class VectorTools const FunctionMap &dirichlet_bc, map &boundary_values); - /** - * Create boundary value information for vector - * valued functions. - * See the other #interpolate_boundary_values#. - */ - static void interpolate_boundary_values (const DoFHandler &dof, - const VectorFunctionMap &dirichlet_bc, - map &boundary_values); - /** * Project #function# to the boundary * of the domain, using the given quadrature @@ -406,9 +387,15 @@ class VectorTools * #boundary_values# contained values * before, the new ones are added, or * the old one overwritten if a node - * of the boundary part to be prjected + * of the boundary part to be projected * on already was in the variable. * + * It is assumed that the number + * of components of the functions + * in #boundary_functions# + * matches that of the finite + * element used by #dof#. + * * See the general documentation of this * class for further information. */ @@ -434,18 +421,28 @@ class VectorTools * accuracy of the #double# data type is * used. * - * The additional argument #weight# allows - * to evaluate weighted norms. This is useful - * for weighting the error of different parts - * differently. A special use is - * to have #weight=0# in some parts of the - * domain, e.g. at - * the location of a shock and #weight=1# - * elsewhere. This allows convergence tests - * in smooth parts of in general discontinuous - * solutions. - * By default, no weighting function is given, - * i.e. weight=1 in the whole domain. + * The additional argument + * #weight# allows to evaluate + * weighted norms. This is useful + * for weighting the error of + * different parts differently. A + * special use is to have + * #weight=0# in some parts of + * the domain, e.g. at the + * location of a shock and + * #weight=1# elsewhere. This + * allows convergence tests in + * smooth parts of in general + * discontinuous solutions. By + * default, no weighting function + * is given, i.e. weight=1 in the + * whole domain. + * + * It is assumed that the number + * of components of the function + * #exact_solution# matches that + * of the finite element used by + * #dof#. * * See the general documentation of this * class for more information. @@ -458,18 +455,6 @@ class VectorTools const NormType &norm, const Function *weight=0); - /** - * Compute the error for the solution of a system. - * See the other #integrate_difference#. - */ - static void integrate_difference (const DoFHandler &dof, - const Vector &fe_function, - const VectorFunction&exact_solution, - Vector &difference, - const Quadrature &q, - const NormType &norm, - const Function *weight=0); - /** * Mean-value filter for Stokes. * The pressure in Stokes' @@ -497,16 +482,18 @@ class VectorTools * Exception */ DeclException0 (ExcNotUseful); - /** * Exception */ DeclException0 (ExcInvalidFE); - /** * Exception */ DeclException0 (ExcInvalidBoundaryIndicator); + /** + * Exception + */ + DeclException0 (ExcComponentMismatch); }; diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index a1ecfd31b6..835c193af7 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -39,7 +39,9 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, { Assert (selected_component < dof.get_fe().n_components, ExcInvalidComponent (selected_component, dof.get_fe().n_components)); - + Assert (coefficient->n_components == 1, + ExcInternalError()); + const unsigned int dim=1; // reserve one slot for each cell and set @@ -68,7 +70,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, // now get the gradients on the // both sides of the point - vector > > gradients (2, vector >(dof.get_fe().n_components)); + vector > > + gradients (2, vector >(dof.get_fe().n_components)); fe_values.reinit (cell); fe_values.get_function_grads (solution, gradients); @@ -83,13 +86,13 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof, } else if (neumann_bc.find(n) != neumann_bc.end()) - grad_neighbor = neumann_bc.find(n)->second->operator()(cell->vertex(0)); + grad_neighbor = neumann_bc.find(n)->second->value(cell->vertex(0)); else grad_neighbor = 0; const double jump = (grad_here - grad_neighbor) * (coefficient != 0 ? - (*coefficient)(cell->vertex(n)) : + coefficient->value(cell->vertex(n)) : 1); error(cell_index) += jump*jump * cell->diameter(); }; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index edbd3f1f2a..c409fc3dfc 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -128,7 +128,9 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & vector &dof_to_boundary_mapping, const Function *a) { const FiniteElement &fe = dof.get_fe(); - + const unsigned int n_components = fe.n_components; + const bool fe_is_system = (n_components != 1); + Assert (matrix.n() == dof.n_boundary_dofs(rhs), ExcInternalError()); Assert (matrix.n() == matrix.m(), ExcInternalError()); Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); @@ -138,20 +140,34 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & Assert (*max_element(dof_to_boundary_mapping.begin(),dof_to_boundary_mapping.end()) == (signed int)matrix.n()-1, ExcInternalError()); + Assert (n_components == rhs.begin()->second->n_components, + ExcComponentMismatch()); const unsigned int dofs_per_cell = fe.total_dofs, dofs_per_face = fe.dofs_per_face; - const unsigned int n_components = fe.n_components; - Assert (n_components == 1, ExcNotImplemented()); - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); Vector cell_vector(dofs_per_cell); UpdateFlags update_flags = UpdateFlags (update_JxW_values | update_q_points); FEFaceValues fe_values (fe, q, update_flags); - + + // two variables for the coefficient, + // one for the two cases indicated in + // the name + vector coefficient_values_scalar (fe_values.n_quadrature_points); + vector > coefficient_values_system (fe_values.n_quadrature_points, + Vector(n_components)); + + vector rhs_values_scalar (fe_values.n_quadrature_points); + vector > rhs_values_system (fe_values.n_quadrature_points, + Vector(n_components)); + + vector dofs (dofs_per_cell); + vector dofs_on_face_vector (dofs_per_face); + set dofs_on_face; + DoFHandler::active_cell_iterator cell = dof.begin_active (), endc = dof.end (); for (; cell!=endc; ++cell) @@ -167,39 +183,95 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & const FullMatrix &values = fe_values.get_shape_values (); const vector &weights = fe_values.get_JxW_values (); - vector rhs_values (fe_values.n_quadrature_points); - rhs.find(cell->face(face)->boundary_indicator()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values); - - if (a != 0) + + if (fe_is_system) + // FE has several components { - vector coefficient_values (fe_values.n_quadrature_points); - a->value_list (fe_values.get_quadrature_points(), coefficient_values); - for (unsigned int point=0; pointface(face)->boundary_indicator()) + ->second->vector_value_list (fe_values.get_quadrature_points(), + rhs_values_system); + + if (a != 0) + { + a->vector_value_list (fe_values.get_quadrature_points(), + coefficient_values_system); + for (unsigned int point=0; pointface(face)->boundary_indicator()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + + if (a != 0) { - for (unsigned int j=0; jvalue_list (fe_values.get_quadrature_points(), + coefficient_values_scalar); + for (unsigned int point=0; point::create_boundary_mass_matrix (const DoFHandler & // inefficient, so we copy the dofs // into a set, which enables binary // searches. - vector dofs (dofs_per_cell); cell->get_dof_indices (dofs); - - vector dofs_on_face_vector (dofs_per_face); cell->face(face)->get_dof_indices (dofs_on_face_vector); - set dofs_on_face (dofs_on_face_vector.begin(), - dofs_on_face_vector.end()); + + dofs_on_face.clear (); + dofs_on_face.insert (dofs_on_face_vector.begin(), + dofs_on_face_vector.end()); + #ifdef DEBUG // in debug mode: compute an element // in the matrix which is @@ -315,8 +387,6 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & - - template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, @@ -593,24 +663,53 @@ void MassMatrix::assemble (FullMatrix &cell_matrix, if (coefficient != 0) { - vector coefficient_values (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int i=0; in_components == 1) + // scalar coefficient given + { + vector coefficient_values (fe_values.n_quadrature_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int i=0; i > coefficient_values (fe_values.n_quadrature_points, + Vector(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int i=0; i::assemble (FullMatrix &cell_matrix, const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components; - // for system elements: need - // VectorFunction for rhs + // for system elements: not + // implemented at present Assert (n_components==1, ExcNotImplemented()); Assert (cell_matrix.n() == total_dofs, @@ -704,8 +803,8 @@ void MassMatrix::assemble (Vector &rhs, const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components; - // for system elements: need - // VectorFunction for rhs + // for system elements: not + // implemented at present Assert (n_components==1, ExcNotImplemented()); Assert (rhs.size() == total_dofs, @@ -749,8 +848,9 @@ void LaplaceMatrix::assemble (FullMatrix &cell_matrix, const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components; - // for system elements: need - // VectorFunction for rhs + // for system elements: might be + // not so useful, not implemented + // at present Assert (n_components==1, ExcNotImplemented()); Assert (cell_matrix.n() == total_dofs, @@ -815,8 +915,9 @@ void LaplaceMatrix::assemble (FullMatrix &cell_matrix, const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components; - // for system elements: need - // VectorFunction for coefficient + // for system elements: might be + // not so useful, not implemented + // at present Assert ((n_components==1) || (coefficient==0), ExcNotImplemented()); Assert (cell_matrix.n() == total_dofs, @@ -870,8 +971,9 @@ void LaplaceMatrix::assemble (Vector &rhs, const FiniteElement &fe = fe_values.get_fe(); const unsigned int n_components = fe.n_components; - // for system elements: need - // VectorFunction for rhs + // for system elements: might be + // not so useful, not implemented + // at present Assert (n_components==1, ExcNotImplemented()); Assert (rhs.size() == total_dofs, diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index ab82810c21..6597ff783b 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -4,7 +4,6 @@ #include -#include #include #include #include @@ -42,56 +41,18 @@ inline double sqr_point (const Tensor<1,dim> &p) { - - template -void VectorTools::interpolate (const DoFHandler &dof, - const Function &function, - Vector &vec) +void VectorTools::interpolate (const DoFHandler &dof, + const Function &function, + Vector &vec) { - const FiniteElement &fe = dof.get_fe(); - - // use #interpolate# function with - // #VectorFunction# param for system - // elements - Assert (fe.n_components == 1, ExcNotUseful()); + Assert (dof.get_fe().n_components == function.n_components, + ExcComponentMismatch()); - DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - vector dofs_on_cell (fe.total_dofs); - vector dof_values_on_cell (fe.total_dofs); - vector > support_points (fe.total_dofs); - for (; cell!=endc; ++cell) - { - // for each cell: - // get location of finite element - // off-points - fe.get_support_points (cell, support_points); - // get function values at these points - function.value_list (support_points, dof_values_on_cell); - // get indices of the dofs on this cell - cell->get_dof_indices (dofs_on_cell); - // distribute function values to the - // whole vector - for (unsigned int i=0; i -void VectorTools::interpolate (const DoFHandler &dof, - const VectorFunction&vectorfunction, - Vector &vec) -{ - const FiniteElement &fe = dof.get_fe(); + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components; + const bool fe_is_system = (n_components != 1); - // use #interpolate# function with - // #Function# param for non-system - // elements - Assert (fe.n_components == vectorfunction.n_components, ExcNotUseful()); - DoFHandler::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); @@ -185,7 +146,7 @@ void VectorTools::interpolate (const DoFHandler &dof, // } // The following is more general. - // It also works if #dofs_per_cell>1#, + // It also works if #dofs_per_x>1#, // i.e. it is usable also for systems // including // FEQ3, FEQ4, FEDG_Qx. @@ -248,8 +209,16 @@ void VectorTools::interpolate (const DoFHandler &dof, vector > support_points (fe.total_dofs); vector > rep_points (n_rep_points); - vector > function_values_at_rep_points ( - n_rep_points, Vector(fe.n_components)); + + // get space for the values of the + // function at the rep support points. + // + // have two versions, one for system fe + // and one for scalar ones, to take the + // more efficient one respectively + vector function_values_scalar (n_rep_points); + vector > function_values_system (n_rep_points, + Vector(fe.n_components)); for (; cell!=endc; ++cell) { @@ -263,26 +232,48 @@ void VectorTools::interpolate (const DoFHandler &dof, for (unsigned int j=0; jget_dof_indices (dofs_on_cell); - // distribute the function values to - // the global vector - for (unsigned int i=0; i void VectorTools::interpolate(const DoFHandler &high_dof, const DoFHandler &low_dof, @@ -342,7 +333,8 @@ void VectorTools::project (const DoFHandler &dof, const Quadrature &q_boundary, const bool project_to_boundary_first) { - Assert (dof.get_fe().n_components == 1, ExcNotUseful()); + Assert (dof.get_fe().n_components == function.n_components, + ExcInvalidFE()); const FiniteElement &fe = dof.get_fe(); @@ -364,6 +356,10 @@ void VectorTools::project (const DoFHandler &dof, for (unsigned int i=0; i::project (const DoFHandler &dof, constraints.condense (mass_matrix); constraints.condense (tmp); - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp); + if (boundary_values.size() != 0) + MatrixTools::apply_boundary_values (boundary_values, + mass_matrix, vec, tmp); SolverControl control(1000,1e-16); PrimitiveVectorMemory > memory; @@ -438,7 +435,8 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof, const Function &rhs, Vector &rhs_vector) { - Assert (dof.get_fe().n_components == 1, ExcNotUseful()); + Assert (dof.get_fe().n_components == rhs.n_components, + ExcComponentMismatch()); UpdateFlags update_flags = UpdateFlags(update_q_points | update_JxW_values); @@ -468,15 +466,17 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof, template <> void VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &dof, - const FunctionMap &dirichlet_bc, - map &boundary_values) + const FunctionMap &dirichlet_bc, + map &boundary_values) { Assert (dirichlet_bc.find(255) == dirichlet_bc.end(), ExcInvalidBoundaryIndicator()); const FiniteElement<1> &fe = dof.get_fe(); - Assert (fe.dofs_per_vertex == 1, ExcInvalidFE()); - Assert (fe.n_components == 1, ExcInvalidFE()); + Assert (fe.n_components == dirichlet_bc.begin()->second->n_components, + ExcComponentMismatch()); + Assert (fe.dofs_per_vertex == fe.n_components, + ExcComponentMismatch()); // check whether boundary values at the // left boundary of the line are requested @@ -495,8 +495,9 @@ VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &dof, // now set the value of the leftmost // degree of freedom - boundary_values[leftmost_cell->vertex_dof_index(0,0)] - = dirichlet_bc.find(0)->second->operator()(leftmost_cell->vertex(0)); + for (unsigned int i=0; ivertex_dof_index(0,i)] + = dirichlet_bc.find(0)->second->value(leftmost_cell->vertex(0), i); }; // same for the right boundary of @@ -516,22 +517,14 @@ VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &dof, // now set the value of the rightmost // degree of freedom - boundary_values[rightmost_cell->vertex_dof_index(1,0)] - = dirichlet_bc.find(1)->second->operator()(rightmost_cell->vertex(1)); + for (unsigned int i=0; ivertex_dof_index(1,i)] + = dirichlet_bc.find(1)->second->value(rightmost_cell->vertex(1), i); }; }; - -template <> -void VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &, - const VectorFunctionMap&, - map&) -{ - Assert (false, ExcNotImplemented()); -}; - #endif @@ -540,64 +533,32 @@ template void VectorTools::interpolate_boundary_values (const DoFHandler &dof, const FunctionMap &dirichlet_bc, - map &boundary_values) { - Assert (dirichlet_bc.find(255) == dirichlet_bc.end(), - ExcInvalidBoundaryIndicator()); - - const FiniteElement &fe = dof.get_fe(); - Assert (fe.dofs_per_vertex == 1, ExcInvalidFE()); - Assert (fe.n_components == 1, ExcInvalidFE()); - - typename FunctionMap::const_iterator function_ptr; - - // field to store the indices of dofs - vector face_dofs (fe.dofs_per_face); - vector > dof_locations (face_dofs.size(), Point()); - vector dof_values (fe.dofs_per_face); - - DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - for (; face!=endf; ++face) - if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) != - dirichlet_bc.end()) - // face is subject to one of the - // bc listed in #dirichlet_bc# - { - // get indices, physical location and - // boundary values of dofs on this - // face - face->get_dof_indices (face_dofs); - fe.get_face_support_points (face, dof_locations); - function_ptr->second->value_list (dof_locations, dof_values); - - // enter into list - for (unsigned int i=0; i -void -VectorTools::interpolate_boundary_values (const DoFHandler &dof, - const VectorFunctionMap &dirichlet_bc, - map &boundary_values) + map &boundary_values) { Assert (dirichlet_bc.find(255) == dirichlet_bc.end(), ExcInvalidBoundaryIndicator()); - const FiniteElement &fe = dof.get_fe(); - Assert (fe.n_components == dirichlet_bc.begin()->second->n_components, + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components; + const bool fe_is_system = (n_components != 1); + + Assert (n_components == dirichlet_bc.begin()->second->n_components, ExcInvalidFE()); - typename VectorFunctionMap::const_iterator function_ptr; + typename FunctionMap::const_iterator function_ptr; // field to store the indices of dofs vector face_dofs (fe.dofs_per_face, -1); vector > dof_locations (face_dofs.size(), Point()); - vector< Vector > dof_values (fe.dofs_per_face, - Vector(fe.n_components)); + // array to store the values of + // the boundary function at the + // boundary points. have to arrays + // for scalar and vector functions + // to use the more efficient one + // respectively + vector dof_values_scalar (fe.dofs_per_face); + vector > dof_values_system (fe.dofs_per_face, + Vector(fe.n_components)); DoFHandler::active_face_iterator face = dof.begin_active_face(), endf = dof.end_face(); @@ -612,13 +573,32 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // face face->get_dof_indices (face_dofs); fe.get_face_support_points (face, dof_locations); - function_ptr->second->value_list (dof_locations, dof_values); - // enter into list - - for (unsigned int i=0; isecond->vector_value_list (dof_locations, dof_values_system); + + // enter into list + + for (unsigned int i=0; isecond->value_list (dof_locations, + dof_values_scalar, + 0); + + // enter into list + + for (unsigned int i=0; i::project_boundary_values (const DoFHandler &dof, const FunctionMap &boundary_functions, const Quadrature &q, map &boundary_values) { - Assert (dof.get_fe().n_components == 1, ExcInvalidFE()); + Assert (dof.get_fe().n_components == boundary_functions.begin()->second->n_components, + ExcComponentMismatch()); vector dof_to_boundary_mapping; dof.map_dof_to_boundary_indices (boundary_functions, dof_to_boundary_mapping); @@ -703,233 +684,24 @@ VectorTools::project_boundary_values (const DoFHandler &dof, - -template -void VectorTools::integrate_difference (const DoFHandler &dof, - const Vector &fe_function, - const Function &exact_solution, - Vector &difference, - const Quadrature &q, - const NormType &norm, - const Function *weight=0) -{ - const FiniteElement &fe = dof.get_fe(); - - difference.reinit (dof.get_tria().n_active_cells()); - - UpdateFlags update_flags = UpdateFlags (update_q_points | - update_JxW_values); - if ((norm==H1_seminorm) || (norm==H1_norm)) - update_flags = UpdateFlags (update_flags | update_gradients); - FEValues fe_values(fe, q, update_flags); - - // loop over all cells - DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (unsigned int index=0; cell != endc; ++cell, ++index) - { - double diff=0; - // initialize for this cell - fe_values.reinit (cell); - - switch (norm) - { - case mean: - case L1_norm: - case L2_norm: - case Linfty_norm: - case H1_norm: - { - // we need the finite element - // function \psi at the different - // integration points. Compute - // it like this: - // \psi(x_j)=\sum_i v_i \phi_i(x_j) - // with v_i the nodal values of the - // fe_function and \phi_i(x_j) the - // matrix of the trial function - // values at the integration point - // x_j. Then the vector - // of the \psi(x_j) is v*Phi with - // v being the vector of nodal - // values on this cell and Phi - // the matrix. - // - // we then need the difference: - // reference_function(x_j)-\psi_j - // and assign that to the vector - // \psi. - const unsigned int n_q_points = q.n_quadrature_points; - vector psi (n_q_points); - - // in praxi: first compute - // exact fe_function vector - exact_solution.value_list (fe_values.get_quadrature_points(), - psi); - // then subtract finite element - // fe_function - if (true) - { - vector function_values (n_q_points, 0); - fe_values.get_function_values (fe_function, function_values); - - transform (psi.begin(), psi.end(), - function_values.begin(), - psi.begin(), - minus()); - }; - - // for L1_norm and Linfty_norm: - // take absolute - // value, for the L2_norm take - // square of psi - switch (norm) - { - case mean: - break; - case L1_norm: - case Linfty_norm: - transform (psi.begin(), psi.end(), - psi.begin(), ptr_fun(fabs)); - break; - case L2_norm: - case H1_norm: - transform (psi.begin(), psi.end(), - psi.begin(), ptr_fun(sqr)); - break; - default: - Assert (false, ExcNotImplemented()); - }; - - // now weight the values - // at the quadrature points due - // to the weighting function - if (weight) - { - vector w(n_q_points); - weight->value_list(fe_values.get_quadrature_points(),w); - for (unsigned int q=0; q > psi (n_q_points); - - // in praxi: first compute - // exact fe_function vector - exact_solution.gradient_list (fe_values.get_quadrature_points(), - psi); - - // then subtract finite element - // fe_function - if (true) - { - vector > function_grads (n_q_points, Tensor<1,dim>()); - fe_values.get_function_grads (fe_function, function_grads); - - transform (psi.begin(), psi.end(), - function_grads.begin(), - psi.begin(), - minus >()); - }; - // take square of integrand - vector psi_square (psi.size(), 0.0); - for (unsigned int i=0; i w(n_q_points); - weight->value_list(fe_values.get_quadrature_points(),w); - for (unsigned int q=0; q void VectorTools::integrate_difference (const DoFHandler &dof, const Vector &fe_function, - const VectorFunction&exact_solution, + const Function &exact_solution, Vector &difference, const Quadrature &q, const NormType &norm, const Function *weight) { - Assert(norm != mean , ExcNotUseful()); + const unsigned int n_q_points = q.n_quadrature_points; + const FiniteElement &fe = dof.get_fe(); + const unsigned int n_components = fe.n_components; + const bool fe_is_system = (n_components != 1); + + Assert( !((n_components == 1) && (norm == mean)), + ExcNotUseful()); - const FiniteElement &fe = dof.get_fe(); - difference.reinit (dof.get_tria().n_active_cells()); UpdateFlags update_flags = UpdateFlags (update_q_points | @@ -937,6 +709,25 @@ VectorTools::integrate_difference (const DoFHandler &dof, if ((norm==H1_seminorm) || (norm==H1_norm)) update_flags = UpdateFlags (update_flags | update_gradients); FEValues fe_values(fe, q, update_flags); + + vector< Vector > function_values (n_q_points, + Vector(n_components)); + vector > > function_grads (n_q_points, + vector >(n_components)); + vector weight_values (n_q_points); + + vector > psi_values (n_q_points, + Vector(n_components)); + vector > > psi_grads (n_q_points, + vector >(n_components)); + vector psi_scalar (n_q_points); + vector psi_square (n_q_points); + + // tmp vector when we use the + // Function functions for + // scalar functions + vector tmp_values (fe_values.n_quadrature_points); + vector > tmp_gradients (fe_values.n_quadrature_points); // loop over all cells DoFHandler::active_cell_iterator cell = dof.begin_active(), @@ -950,30 +741,34 @@ VectorTools::integrate_difference (const DoFHandler &dof, switch (norm) { case mean: - break; case L1_norm: case L2_norm: case Linfty_norm: case H1_norm: { - const unsigned int n_q_points = q.n_quadrature_points; - vector > psi (n_q_points, Vector(fe.n_components)); - // first compute the exact solution // (vectors) at the quadrature points - exact_solution.value_list (fe_values.get_quadrature_points(), psi); + // try to do this as efficient as + // possible by avoiding a second + // virtual function call in case + // the function really has only + // one component + if (fe_is_system) + exact_solution.vector_value_list (fe_values.get_quadrature_points(), + psi_values); + else + { + exact_solution.value_list (fe_values.get_quadrature_points(), + tmp_values); + for (unsigned int i=0; i > function_values ( - n_q_points, Vector(fe.n_components)); - - fe_values.get_function_values (fe_function, function_values); - - for (unsigned int q=0; q::integrate_difference (const DoFHandler &dof, // Use psi_scalar to store the squares // of the vectors or the vector norms // respectively. - vector psi_scalar (n_q_points); switch (norm) { case mean: @@ -994,7 +788,7 @@ VectorTools::integrate_difference (const DoFHandler &dof, case L2_norm: case H1_norm: for (unsigned int q=0; q::integrate_difference (const DoFHandler &dof, // to the weighting function if (weight) { - vector w(n_q_points); - weight->value_list(fe_values.get_quadrature_points(),w); + weight->value_list(fe_values.get_quadrature_points(), + weight_values); for (unsigned int q=0; q::integrate_difference (const DoFHandler &dof, switch (norm) { case mean: - break; case L1_norm: case L2_norm: case H1_norm: @@ -1057,45 +850,48 @@ VectorTools::integrate_difference (const DoFHandler &dof, // Until now, #diff# includes the // square of the L2_norm. - // same procedure as above, but now - // psi is a vector of Jacobians - // i.e. psi is a vector of vectors of - // gradients. - const unsigned int n_q_points = q.n_quadrature_points; - vector > > psi ( - n_q_points, vector >(fe.n_components, Tensor<1,dim>())); - // in praxi: first compute // exact fe_function vector - exact_solution.gradient_list (fe_values.get_quadrature_points(), psi); + // + // try to be a little clever + // to avoid recursive virtual + // function calls when calling + // #gradient_list# for functions + // that are really scalar + // functions + if (fe_is_system) + exact_solution.vector_gradient_list (fe_values.get_quadrature_points(), + psi_grads); + else + { + exact_solution.gradient_list (fe_values.get_quadrature_points(), + tmp_gradients); + for (unsigned int i=0; i > > function_grads ( - n_q_points, vector >(fe.n_components, Tensor<1,dim>())); - fe_values.get_function_grads (fe_function, function_grads); + fe_values.get_function_grads (fe_function, function_grads); + for (unsigned int k=0; k psi_square (psi.size(), 0.0); - for (unsigned int q=0; q w(n_q_points); - weight->value_list(fe_values.get_quadrature_points(),w); + weight->value_list(fe_values.get_quadrature_points(), + weight_values); for (unsigned int q=0; q { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; @@ -106,19 +107,22 @@ class Solution : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; /** * Return the gradient of the function * at the given point. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; template <> -double RHSPoly<2>::operator () (const Point<2> &p) const { +double RHSPoly<2>::value (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -128,7 +132,8 @@ double RHSPoly<2>::operator () (const Point<2> &p) const { template <> -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::value (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -137,7 +142,8 @@ double Solution<2>::operator () (const Point<2> &p) const { template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::gradient (const Point<2> &p, + const unsigned int) const { const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -162,7 +168,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, fe_values.shape_grad(j,point)) * fe_values.JxW(point); rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -376,8 +382,11 @@ int PoissonProblem::run (const unsigned int level) { if (dof->n_dofs()<=5000) { - Vector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; - Vector h1_seminorm_error_per_dof, h1_error_per_dof; + Vector l1_error_per_dof(dof->n_dofs()); + Vector l2_error_per_dof(dof->n_dofs()); + Vector linfty_error_per_dof(dof->n_dofs()); + Vector h1_seminorm_error_per_dof(dof->n_dofs()); + Vector h1_error_per_dof(dof->n_dofs()); dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/tests/big-tests/error-estimation/error-estimation.cc b/tests/big-tests/error-estimation/error-estimation.cc index 5158c8c901..aebd9cf999 100644 --- a/tests/big-tests/error-estimation/error-estimation.cc +++ b/tests/big-tests/error-estimation/error-estimation.cc @@ -114,25 +114,32 @@ class Solution { class GaussShape : public Function { public: - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; class Singular : public Function { public: - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; class Kink : public Function { public: class Coefficient : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; - virtual double operator () (const Point &p) const; - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; }; @@ -150,7 +157,8 @@ class RHS { */ class GaussShape : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; /** @@ -160,7 +168,8 @@ class RHS { */ class Singular : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; /** @@ -171,7 +180,8 @@ class RHS { */ class Kink : public Function { public: - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; }; @@ -179,13 +189,15 @@ class RHS { template <> -double Solution<2>::GaussShape::operator () (const Point<2> &p) const { +double Solution<2>::GaussShape::value (const Point<2> &p, + const unsigned int) const { return p(0)*p(1)*exp(-40*p.square()); }; template <> -Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p, + const unsigned int) const { return Point<2> ((1-80.*p(0)*p(0))*p(1)*exp(-40*p.square()), (1-80.*p(1)*p(1))*p(0)*exp(-40*p.square())); }; @@ -193,13 +205,15 @@ Tensor<1,2> Solution<2>::GaussShape::gradient (const Point<2> &p) const { template <> -double Solution<2>::Singular::operator () (const Point<2> &p) const { +double Solution<2>::Singular::value (const Point<2> &p, + const unsigned int) const { return pow(p.square(), 1./3.); }; template <> -Tensor<1,2> Solution<2>::Singular::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::Singular::gradient (const Point<2> &p, + const unsigned int) const { return 2./3.*pow(p.square(), -2./3.) * p; }; @@ -213,21 +227,24 @@ inline double theta(const double x) { template <> -double Solution<2>::Kink::operator () (const Point<2> &p) const { +double Solution<2>::Kink::value (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return (1+4*theta(s))*s; }; template <> -Tensor<1,2> Solution<2>::Kink::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::Kink::gradient (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return (1+4*theta(s))*Point<2>(-2*p(0),1); }; template <> -double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { +double Solution<2>::Kink::Coefficient::value (const Point<2> &p, + const unsigned int) const { const double s = p(1)-p(0)*p(0); return 1./(1.+4.*theta(s)); }; @@ -235,19 +252,22 @@ double Solution<2>::Kink::Coefficient::operator () (const Point<2> &p) const { template <> -double RHS<2>::GaussShape::operator () (const Point<2> &p) const { +double RHS<2>::GaussShape::value (const Point<2> &p, + const unsigned int) const { return (480.-6400.*p.square())*p(0)*p(1)*exp(-40.*p.square()); }; template <> -double RHS<2>::Singular::operator () (const Point<2> &p) const { +double RHS<2>::Singular::value (const Point<2> &p, + const unsigned int) const { return -4./9. * pow(p.square(), -2./3.); }; template <> -double RHS<2>::Kink::operator () (const Point<2> &) const { +double RHS<2>::Kink::value (const Point<2> &, + const unsigned int) const { return 2; }; @@ -266,7 +286,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, for (unsigned int point=0; point::assemble (FullMatrix &cell_matrix, fe_values.JxW(point) * c; rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -535,8 +555,8 @@ void PoissonProblem::run (ParameterHandler &prm) { cout << estimated_error_per_cell.l2_norm() << endl; estimated_error.push_back (estimated_error_per_cell.l2_norm()); - Vector l2_error_per_dof, linfty_error_per_dof; - Vector h1_error_per_dof, estimated_error_per_dof; + Vector l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs()); + Vector h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs()); Vector error_ratio (dof->n_dofs()); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc index 572a209dea..317a7282c9 100644 --- a/tests/big-tests/multigrid/multigrid.cc +++ b/tests/big-tests/multigrid/multigrid.cc @@ -174,7 +174,8 @@ class RHSPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; }; @@ -186,19 +187,24 @@ class Solution : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int component) const; /** * Return the gradient of the function * at the given point. */ - virtual Tensor<1,dim> gradient (const Point &p) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component) const; }; template <> -double RHSPoly<2>::operator () (const Point<2> &p) const { +double RHSPoly<2>::value (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -208,7 +214,10 @@ double RHSPoly<2>::operator () (const Point<2> &p) const { template <> -double Solution<2>::operator () (const Point<2> &p) const { +double Solution<2>::value (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -217,7 +226,10 @@ double Solution<2>::operator () (const Point<2> &p) const { template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p) const { +Tensor<1,2> Solution<2>::gradient (const Point<2> &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double x = p(0), y = p(1); const double pi= 3.1415926536; @@ -242,7 +254,7 @@ void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, fe_values.shape_grad(j,point)) * fe_values.JxW(point); rhs(i) += fe_values.shape_value(i,point) * - right_hand_side(fe_values.quadrature_point(point)) * + right_hand_side.value(fe_values.quadrature_point(point)) * fe_values.JxW(point); }; }; @@ -534,8 +546,11 @@ int PoissonProblem::run (const unsigned int level) { if (dof->DoFHandler::n_dofs()<=5000) { - Vector l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof; - Vector h1_seminorm_error_per_dof, h1_error_per_dof; + Vector l1_error_per_dof (dof->DoFHandler::n_dofs()); + Vector l2_error_per_dof (dof->DoFHandler::n_dofs()); + Vector linfty_error_per_dof (dof->DoFHandler::n_dofs()); + Vector h1_seminorm_error_per_dof (dof->DoFHandler::n_dofs()); + Vector h1_error_per_dof (dof->DoFHandler::n_dofs()); dof->distribute_cell_to_dof_vector (l1_error_per_cell, l1_error_per_dof); dof->distribute_cell_to_dof_vector (l2_error_per_cell, l2_error_per_dof); dof->distribute_cell_to_dof_vector (linfty_error_per_cell, diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc index dd4053a9e2..0034eb8e9f 100644 --- a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc +++ b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc @@ -33,7 +33,7 @@ template class RightHandSide : public Function { public: - double operator () (const Point &p) const + double value (const Point &p) const { double x = 80; for (unsigned int d=0; d { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const { + virtual double value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + double x = 1; for (unsigned int i=0; i &p, + Vector &values) const { + Assert (values.size()==1, ExcVectorHasWrongSize (values.size(), 1)); + + double x = 1; + + for (unsigned int i=0; i { * empty. */ virtual void value_list (const vector > &points, - vector &values) const { + vector &values, + const unsigned int component) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); for (unsigned int i=0; i::operator() (points[i]); + values[i] = BoundaryValuesSine::value (points[i], component); }; }; @@ -49,7 +69,9 @@ class BoundaryValuesJump : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const { + virtual double value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); switch (dim) { case 1: @@ -73,7 +95,8 @@ class RHSTrigPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int) const; }; @@ -89,7 +112,8 @@ class RHSPoly : public Function { * Return the value of the function * at the given point. */ - virtual double operator () (const Point &p) const; + virtual double value (const Point &p, + const unsigned int) const; }; @@ -200,7 +224,10 @@ CurvedLine::get_new_point_on_quad (const typename Triangulation::quad_ template -double RHSTrigPoly::operator () (const Point &p) const { +double RHSTrigPoly::value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + const double pi = 3.1415926536; switch (dim) { @@ -219,7 +246,10 @@ double RHSTrigPoly::operator () (const Point &p) const { template -double RHSPoly::operator () (const Point &p) const { +double RHSPoly::value (const Point &p, + const unsigned int component) const { + Assert (component==0, ExcIndexRange (component, 0, 1)); + double ret_val = 0; for (unsigned int i=0; i