]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge unified_function_and_vectorfunction (mergepoint 1) into the main branch.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Oct 1999 11:19:00 +0000 (11:19 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Oct 1999 11:19:00 +0000 (11:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@1793 0785d39b-7218-0410-832d-ea1e28bc413d

23 files changed:
deal.II/base/Todo
deal.II/base/include/base/function.h
deal.II/base/include/base/tensor_function.h
deal.II/base/include/base/tensorindex.h
deal.II/base/source/function.cc
deal.II/base/source/tensor_function.cc
deal.II/deal.II/Attic/examples/convergence/convergence.cc
deal.II/deal.II/Attic/examples/error-estimation/error-estimation.cc
deal.II/deal.II/Attic/examples/multigrid/multigrid.cc
deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc
deal.II/deal.II/Attic/examples/poisson/problem.cc
deal.II/deal.II/include/grid/intergrid_map.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/source/numerics/error_estimator.cc
deal.II/deal.II/source/numerics/matrices.cc
deal.II/deal.II/source/numerics/vectors.cc
tests/big-tests/convergence/convergence.cc
tests/big-tests/error-estimation/error-estimation.cc
tests/big-tests/multigrid/multigrid.cc
tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc
tests/big-tests/poisson/problem.cc

index f6c223c79e44cd6d889aeccaa5772ae55aa36605..3f13b27f5af7d01702b2582892deea8a0a7a3502 100644 (file)
@@ -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.
index 48ddc05c260994d46a0dcfad99c3b42be27c3dc3..11bb0ccf3feebfe03c8c8ab78d2eff2e1a2efe77 100644 (file)
 /*----------------------------   function.h     ---------------------------*/
 
 
+#include <base/forward-declarations.h>
 #include <base/exceptions.h>
-#include <vector>
-#include <base/point.h>
 #include <base/functiontime.h>
-
+#include <base/subscriptor.h>
+#include <vector>
 
 
 
 /**
- *  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<dim>   &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<dim>   &p,
+ *                        Vector<double>     &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<dim> >  &point_list,
+ *                      vector<double>             &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<dim> >    &point_list,
+ *                             vector<Vector<double> >      &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 <int dim>
-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<dim> &p) const;
+    virtual double value (const Point<dim>   &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<dim>   &p,
+                                Vector<double>     &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<Point<dim> > &points,
-                            vector<double>            &values) const;
+                            vector<double>            &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<dim> &p) const;
+    virtual void vector_value_list (const vector<Point<dim> > &points,
+                                   vector<Vector<double> >   &values) const;
+    
+                                    /**
+                                     * Return the gradient of the
+                                     * specified component of the
+                                     * function at the given point.
+                                     */
+    virtual Tensor<1,dim> gradient (const Point<dim>   &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<dim>       &p,
+                                          vector<Tensor<1,dim> > &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<Point<dim> > &points,
+                               vector<Tensor<1,dim> >    &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<Point<dim> > &points,
-                               vector<Tensor<1,dim> >    &gradients) const;
-    
+    virtual void vector_gradient_list (const vector<Point<dim> >       &points,
+                                      vector<vector<Tensor<1,dim> > > &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 <int dim>
 class ZeroFunction : public Function<dim> {
   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<dim>   &p,
+                         const unsigned int  component) const;
+
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point for all
+                                     * components.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual void   vector_value (const Point<dim> &p,
+                                Vector<double>   &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<Point<dim> > &points,
-                            vector<double>            &values) const;
+                            vector<double>            &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<Point<dim> > &points,
+                                   vector<Vector<double> >   &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<dim> &p) const;
+    virtual Tensor<1,dim> gradient (const Point<dim> &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<dim>       &p,
+                                          vector<Tensor<1,dim> > &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<Point<dim> > &points,
-                               vector<Tensor<1,dim> >    &gradients) const;
+                               vector<Tensor<1,dim> >    &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<Point<dim> >       &points,
+                                      vector<vector<Tensor<1,dim> > > &gradients) const;
 };
 
 
@@ -159,46 +354,83 @@ class ZeroFunction : public Function<dim> {
 
 
 /**
- *  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 <int dim>
 class ConstantFunction : public ZeroFunction<dim> {
   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<dim> &p) const;
+    virtual double value (const Point<dim>   &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<dim> &p,
+                                Vector<double>   &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<Point<dim> > &points,
-                            vector<double>            &values) const;
+                            vector<double>            &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<Point<dim> > &points,
+                                   vector<Vector<double> >   &values) const;
 
   protected:
                                     /**
index e79c73d21bc28ba333c67178658e5428e47ba696..6d0cacbc78ed4a8bd9f1c9f23a9b4cd2ea79eef9 100644 (file)
 #include <base/function.h>
 #include <base/point.h>
 #include <base/functiontime.h>
-#include <base/tensorindex.h>
 #include <base/forward-declarations.h>
 
-template <typename number> class Vector;
-template <int dim> class VectorFunction;
-template <int rank, int dim> 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 <int dim>
-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<dim> &p,
-                       Vector<double>   &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<Point<dim> > &points,
-                            vector<Vector<double> >   &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<Point<dim> >       &points,
-                               vector<vector<Tensor<1,dim> > > &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<dim>
-    {
-      public:
-                                        /**
-                                         * Constructor.
-                                         * The arguments are the #VectorFunction# to be
-                                         * accessed and the component index.
-                                         */
-       Extractor(const VectorFunction<dim>& f, unsigned int index);
-
-                                        /**
-                                         * Compute function value.
-                                         */
-       virtual double operator() (const Point<dim>& p) const;
-
-                                        /**
-                                         * Compute several values.
-                                         */
-       virtual void value_list (const vector<Point<dim> > &points,
-                                vector<double> &values) const;
-       
-
-                                        /**
-                                         * Compute derivative.
-                                         */
-       virtual Tensor<1,dim> gradient (const Point<dim>& p) const;
-
-                                        /**
-                                         * Compute several derivatives.
-                                         */
-       virtual void gradient_list (const vector<Point<dim> > &points,
-                                   vector<Tensor<1,dim> > &gradients) const;
-       
-      private:
-                                        /**
-                                         * Pointer to the #VectorFunction#.
-                                         */
-       const SmartPointer<VectorFunction<dim> > 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 <int rank, int dim>
-class TensorFunction : public VectorFunction<dim>
+class TensorFunction : public FunctionTime,
+                      public Subscriptor
 {
   public:
                                     /**
@@ -204,7 +60,7 @@ class TensorFunction : public VectorFunction<dim>
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual Tensor<rank, dim> operator () (const Point<dim> &p) const;
+    virtual Tensor<rank, dim> value (const Point<dim> &p) const;
 
                                     /**
                                      * Set #values# to the point values
@@ -218,8 +74,8 @@ class TensorFunction : public VectorFunction<dim>
                             vector<Tensor<rank,dim> > &values) const;
 
                                     /**
-                                     * Return the gradient of the function
-                                     * at the given point.
+                                     * Return the gradient of the
+                                     * function at the given point.
                                      */
     virtual Tensor<rank+1,dim> gradient (const Point<dim> &p) const;
 
@@ -233,24 +89,6 @@ class TensorFunction : public VectorFunction<dim>
     virtual void gradient_list (const vector<Point<dim> > &points,
                                vector<Tensor<rank+1,dim> > &gradients) const;
 
-                                    /**
-                                     * See #VectorFunction#.
-                                     */  
-    virtual void value (const Point<dim> &points,
-                            Vector<double> &values) const;
-
-                                    /**
-                                     * See #VectorFunction#.
-                                     */  
-    virtual void value_list (const vector<Point<dim> > &points,
-                            vector<Vector<double> > &values) const;
-
-                                    /**
-                                     * See #VectorFunction#.
-                                     */  
-    virtual void gradient_list (const vector<Point<dim> > &points,
-                               vector<vector<Tensor<1,dim> > > &gradients) const;
-
                                     /**
                                      * Exception
                                      */
index bf4e0b598a0391846a2c64c9b96424075e2eb651..b9d0ddb81bf75eb1898ce1506832e0f8e99ea5a4 100644 (file)
@@ -6,8 +6,8 @@
 
 #include <base/exceptions.h>
 
+
 /**
- *
  * Rank-independent access to elements of #Tensor#.  A little class
  * template remembering #rank# integers.
  *
  * or contact the developer.
  *
  * @author Guido Kanschat, 1999
- *
  */
 template<int rank>
 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");  
 };
index d2905c2a5b5504dbef3af5d53569bb15be0075fe..9d0da2c5c2dbfc8dfabb178882f43797886776f5 100644 (file)
@@ -2,12 +2,16 @@
 
 
 #include <base/function.h>
+#include <base/point.h>
+#include <lac/vector.h>
 #include <vector>
 
 
 template <int dim>
-Function<dim>::Function (const double initial_time) :
-               FunctionTime(initial_time)
+Function<dim>::Function (const unsigned int n_components,
+                        const double       initial_time) :
+               FunctionTime(initial_time),
+               n_components(n_components)
 {};
 
 
@@ -19,46 +23,116 @@ Function<dim>::~Function ()
 
 
 template <int dim>
-double Function<dim>::operator () (const Point<dim> &) const {
+double Function<dim>::value (const Point<dim> &,
+                            const unsigned int) const
+{
   Assert (false, ExcPureFunctionCalled());
   return 0;
 };
 
 
 
+template <int dim>
+void Function<dim>::vector_value (const Point<dim> &,
+                                 Vector<double>   &) const
+{
+  Assert (false, ExcPureFunctionCalled());
+};
+
+
+
 template <int dim>
 void Function<dim>::value_list (const vector<Point<dim> > &points,
-                               vector<double>            &values) const {
+                               vector<double>            &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; i<points.size(); ++i)
-    values[i]  = this->operator() (points[i]);
+    values[i]  = this->value (points[i], component);
+};
+
+
+
+template <int dim>
+void Function<dim>::vector_value_list (const vector<Point<dim> > &points,
+                                      vector<Vector<double> >   &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; i<points.size(); ++i)
+    this->vector_value (points[i], values[i]);
 };
 
 
 
 
 template <int dim>
-Tensor<1,dim> Function<dim>::gradient (const Point<dim> &) const {
+Tensor<1,dim> Function<dim>::gradient (const Point<dim> &,
+                                      const unsigned int) const
+{
   Assert (false, ExcPureFunctionCalled());
   return Point<dim>();
 };
 
 
 
+template <int dim>
+void Function<dim>::vector_gradient (const Point<dim>       &,
+                                    vector<Tensor<1,dim> > &) const
+{
+  Assert (false, ExcPureFunctionCalled());
+};
+
+
+
 template <int dim>
 void Function<dim>::gradient_list (const vector<Point<dim> > &points,
-                                  vector<Tensor<1,dim> >    &gradients) const {
+                                  vector<Tensor<1,dim> >    &gradients,
+                                  const unsigned int         component) const
+{
   Assert (gradients.size() == points.size(),
          ExcVectorHasWrongSize(gradients.size(), points.size()));
 
   for (unsigned int i=0; i<points.size(); ++i)
-    gradients[i] = gradient(points[i]);
+    gradients[i] = gradient(points[i], component);
 };
 
 
 
+template <int dim>
+void Function<dim>::vector_gradient_list (const vector<Point<dim> >       &points,
+                                         vector<vector<Tensor<1,dim> > > &gradients) const
+{
+  Assert (gradients.size() == points.size(),
+         ExcVectorHasWrongSize(gradients.size(), points.size()));
+
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      Assert (gradients[i].size() == n_components,
+             ExcVectorHasWrongSize(gradients[i].size(), n_components));
+      for (unsigned int component=0; component<n_components; ++component)
+       gradients[i][component] = gradient(points[i], component);
+    };
+};
+
+
+
+
+
+template <int dim>
+ZeroFunction<dim>::ZeroFunction (const unsigned int n_components) :
+               Function<dim> (n_components)
+{};
+
 
 
 template <int dim>
@@ -67,15 +141,30 @@ ZeroFunction<dim>::~ZeroFunction () {};
 
 
 template <int dim>
-double ZeroFunction<dim>::operator () (const Point<dim> &) const {
+double ZeroFunction<dim>::value (const Point<dim> &,
+                                      const unsigned int) const
+{
   return 0.;
 };
 
 
 
+template <int dim>
+void ZeroFunction<dim>::vector_value (const Point<dim> &,
+                                     Vector<double>   &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 <int dim>
 void ZeroFunction<dim>::value_list (const vector<Point<dim> > &points,
-                                   vector<double>            &values) const {
+                                   vector<double>            &values,
+                                   const unsigned int         /*component*/) const {
   Assert (values.size() == points.size(),
          ExcVectorHasWrongSize(values.size(), points.size()));
 
@@ -85,27 +174,82 @@ void ZeroFunction<dim>::value_list (const vector<Point<dim> > &points,
 
 
 template <int dim>
-Tensor<1,dim> ZeroFunction<dim>::gradient (const Point<dim> &) const {
+void ZeroFunction<dim>::vector_value_list (const vector<Point<dim> > &points,
+                                          vector<Vector<double> >   &values) const
+{
+  Assert (values.size() == points.size(),
+         ExcVectorHasWrongSize(values.size(), points.size()));
+
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      Assert (values[i].size() == n_components,
+             ExcVectorHasWrongSize(values[i].size(), n_components));
+      fill_n (values[i].begin(), n_components, 0.);
+    };
+};
+
+
+
+template <int dim>
+Tensor<1,dim> ZeroFunction<dim>::gradient (const Point<dim> &,
+                                          const unsigned int) const
+{
   return Tensor<1,dim>();
 };
 
 
 
+template <int dim>
+void ZeroFunction<dim>::vector_gradient (const Point<dim>       &,
+                                        vector<Tensor<1,dim> > &gradients) const
+{
+  Assert (gradients.size() == n_components,
+         ExcVectorHasWrongSize(gradients.size(), n_components));
+
+  for (unsigned int c=0; c<n_components; ++c)
+    gradients[c].clear ();
+};
+
+
+
 template <int dim>
 void ZeroFunction<dim>::gradient_list (const vector<Point<dim> > &points,
-                                      vector<Tensor<1,dim> >    &gradients) const {
+                                      vector<Tensor<1,dim> >    &gradients,
+                                      const unsigned int         /*component*/) const
+{
   Assert (gradients.size() == points.size(),
          ExcVectorHasWrongSize(gradients.size(), points.size()));
 
-  gradients.clear ();
+  for (unsigned int i=0; i<points.size(); ++i)
+    gradients[i].clear ();
+};
+
+
+
+template <int dim>
+void ZeroFunction<dim>::vector_gradient_list (const vector<Point<dim> >       &points,
+                                             vector<vector<Tensor<1,dim> > > &gradients) const
+{
+  Assert (gradients.size() == points.size(),
+         ExcVectorHasWrongSize(gradients.size(), points.size()));
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      Assert (gradients[i].size() == n_components,
+             ExcVectorHasWrongSize(gradients[i].size(), n_components));
+      for (unsigned int c=0; c<n_components; ++c)
+       gradients[i][c].clear ();
+    };
 };
 
 
 
 
+
 template <int dim>
-ConstantFunction<dim>::ConstantFunction (const double value) :
-               function_value(value) {};
+ConstantFunction<dim>::ConstantFunction (const double value,
+                                        const unsigned int n_components) :
+               ZeroFunction<dim> (n_components),
+               function_value    (value) {};
 
 
 template <int dim>
@@ -114,15 +258,30 @@ ConstantFunction<dim>::~ConstantFunction () {};
 
 
 template <int dim>
-double ConstantFunction<dim>::operator () (const Point<dim> &) const {
+double ConstantFunction<dim>::value (const Point<dim> &,
+                                          const unsigned int) const
+{
   return function_value;
 };
 
 
 
+template <int dim>
+void ConstantFunction<dim>::vector_value (const Point<dim> &,
+                                         Vector<double>   &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 <int dim>
 void ConstantFunction<dim>::value_list (const vector<Point<dim> > &points,
-                                   vector<double>            &values) const {
+                                       vector<double>            &values,
+                                       const unsigned int         /*component*/) const {
   Assert (values.size() == points.size(),
          ExcVectorHasWrongSize(values.size(), points.size()));
 
@@ -131,6 +290,24 @@ void ConstantFunction<dim>::value_list (const vector<Point<dim> > &points,
 
 
 
+template <int dim>
+void ConstantFunction<dim>::vector_value_list (const vector<Point<dim> > &points,
+                                              vector<Vector<double> >   &values) const
+{
+  Assert (values.size() == points.size(),
+         ExcVectorHasWrongSize(values.size(), points.size()));
+
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      Assert (values[i].size() == n_components,
+             ExcVectorHasWrongSize(values[i].size(), n_components));
+      fill_n (values[i].begin(), n_components, function_value);
+    };
+};
+
+
+
+
 // explicit instantiations
 
 template class Function<1>;
index 6b7a11ac25693e75b13495a77d4412dca578a055..34801c9c77734bad6281dc8bda0df437ecae741d 100644 (file)
@@ -7,104 +7,7 @@
 #include <cmath>
 #include <lac/vector.h>
 
-template <int dim>
-VectorFunction<dim>::VectorFunction(unsigned n_components, const double initial_time)
-               :
-               FunctionTime(initial_time),
-               n_components(n_components)
-{}
-
-
-template <int dim>
-VectorFunction<dim>::~VectorFunction()
-{}
-
-/*
-template <int dim> double
-VectorFunction<dim>::operator () (const Point<dim> &, unsigned) const
-
-{
-  Assert (false, ExcPureFunctionCalled());
-  return 0.;
-}
-*/
-
-template <int dim>
-void
-VectorFunction<dim>::value (const Point<dim>  &, Vector<double> &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
-
-
-template <int dim>
-void
-VectorFunction<dim>::value_list (const vector<Point<dim> > &ps,
-                                vector<Vector<double> > &us) const
-{
-  for (unsigned int i=0 ; i<ps.size() ; ++i)
-    value(ps[i], us[i]);
-}
-
-
-template <int dim>
-void
-VectorFunction<dim>::gradient_list (const vector<Point<dim> > &,
-                                   vector<vector<Tensor<1,dim> > > &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
 
-template <int dim>
-VectorFunction<dim>::Extractor::Extractor(const VectorFunction<dim>& f,
-                                         unsigned int index)
-               :
-               vectorfunction(f),
-               index(index)
-{}
-
-template <int dim>
-double
-VectorFunction<dim>::Extractor::operator() (const Point<dim>& p) const
-{
-  Vector<double> v(vectorfunction->n_components);
-  vectorfunction->value(p,v);
-  return v(index);
-}
-
-
-template <int dim>
-Tensor<1,dim>
-VectorFunction<dim>::Extractor::gradient (const Point<dim>&) const
-{
-  Assert(false, ExcNotImplemented());
-  return Tensor<1,dim>();
-}
-
-template <int dim>
-void
-VectorFunction<dim>::Extractor::value_list (const vector<Point<dim> > &points,
-                                           vector<double> &values) const
-{
-  vector<Vector<double> > v(values.size(),
-                           Vector<double>(vectorfunction->n_components));
-  vectorfunction->value_list(p,v);
-  for (unsigned int i=0 ; i<values.size() ; ++i)
-    values[i] = v[i](index);
-}
-
-
-template <int dim>
-void
-VectorFunction<dim>::Extractor::gradient_list (const vector<Point<dim> > &points,
-                                              vector<Tensor<1,dim> > &gradients) const
-{
-  vector<vector<Tensor<1,dim> > > v(values.size(),
-                           vector<Tensor<1,dim> >(vectorfunction->n_components));
-  vectorfunction->value_list(p,v);
-  for (unsigned int i=0 ; i<values.size() ; ++i)
-    values[i] = v[i][index];
-}
 
 //////////////////////////////////////////////////////////////////////
 // TensorFunction
@@ -113,7 +16,7 @@ VectorFunction<dim>::Extractor::gradient_list (const vector<Point<dim> > &points
 template <int rank, int dim>
 TensorFunction<rank, dim>::TensorFunction (const double initial_time)
                :
-               VectorFunction<dim>(pow(dim,rank), initial_time)
+               FunctionTime (initial_time)
 {};
 
 
@@ -124,41 +27,30 @@ TensorFunction<rank, dim>::~TensorFunction ()
 
 
 
-// template <int rank, int dim>
-// double
-// TensorFunction<rank, dim>::operator () (TensorIndex<rank> i,
-//                                      const Point<dim> &) const
-// {
-//   int k=i(0);
-//   k++;
-  
-//   Assert (false, ExcPureFunctionCalled());
-//   return 0;
-// };
-
-
 template <int rank, int dim>
 Tensor<rank,dim>
-TensorFunction<rank, dim>::operator() (const Point<dim> &) const
+TensorFunction<rank, dim>::value (const Point<dim> &) const
 {
   Assert (false, ExcPureFunctionCalled());
   return Tensor<rank,dim>();
 };
 
 
+
 template <int rank, int dim>
 void
 TensorFunction<rank, dim>::value_list (const vector<Point<dim> > &points,
-                                vector<Tensor<rank,dim> > &values) const
+                                      vector<Tensor<rank,dim> > &values) const
 {
   Assert (values.size() == points.size(),
          ExcVectorHasWrongSize(values.size(), points.size()));
 
   for (unsigned int i=0; i<points.size(); ++i)
-    values[i]  = this->operator() (points[i]);
+    values[i]  = this->value (points[i]);
 };
 
 
+
 template <int rank, int dim>
 Tensor<rank+1,dim>
 TensorFunction<rank, dim>::gradient (const Point<dim> &) const
@@ -182,36 +74,6 @@ TensorFunction<rank, dim>::gradient_list (const vector<Point<dim> > &points,
 };
 
 
-template <int rank, int dim> void
-TensorFunction<rank, dim>::value (const Point<dim>  &p,
-                                  Vector<double> &erg) const
-{
-  Tensor<rank,dim> h = operator()(p);
-  h.unroll(erg);
-}
-
-
-template <int rank, int dim> void
-TensorFunction<rank, dim>::value_list (const vector<Point<dim> > & points,
-                                vector<Vector<double> > & values) const
-{
-  Assert (values.size() == points.size(),
-         ExcVectorHasWrongSize(values.size(), points.size()));
-
-  for (unsigned int i=0; i<points.size(); ++i)
-    operator() (points[i]).unroll(values[i]);
-  
-}
-
-
-template <int rank, int dim> void
-TensorFunction<rank, dim>::gradient_list (const vector<Point<dim> > &,
-                                   vector<vector<Tensor<1,dim> > > &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
-
-
 
 template class TensorFunction<1,1>;
 template class TensorFunction<2,1>;
index 1a5ad43c5c4f5166c420d38dd692b7a93b528fb8..afe9ac1f4b4fae4deb4d485e251169ceffb96905 100644 (file)
@@ -94,7 +94,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
 };
 
 
@@ -106,19 +107,22 @@ class Solution : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
                                     /**
                                      * Return the gradient of the function
                                      * at the given point.
                                      */
-    virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+    virtual Tensor<1,dim> gradient (const Point<dim> &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<double>  &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<dim>::run (const unsigned int level) {
 
   if (dof->n_dofs()<=5000) 
     {
-      Vector<double> l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_seminorm_error_per_dof, h1_error_per_dof;
+      Vector<double> l1_error_per_dof(dof->n_dofs());
+      Vector<double> l2_error_per_dof(dof->n_dofs());
+      Vector<double> linfty_error_per_dof(dof->n_dofs());
+      Vector<double> h1_seminorm_error_per_dof(dof->n_dofs());
+      Vector<double> 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,
index 5158c8c901bf0015c0ab5bc20c5ff309c4e1ac71..aebd9cf99929204c72355474039bad510cff107c 100644 (file)
@@ -114,25 +114,32 @@ class Solution {
 
     class GaussShape : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 
     class Singular : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 
     class Kink : public Function<dim> {
       public:
        class Coefficient : public Function<dim> {
          public:
-           virtual double operator () (const Point<dim> &p) const;
+           virtual double value (const Point<dim> &p,
+                                 const unsigned int component) const;
        };
        
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 };
 
@@ -150,7 +157,8 @@ class RHS {
                                      */
     class GaussShape : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
     };
 
                                     /**
@@ -160,7 +168,8 @@ class RHS {
                                      */
     class Singular : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
     };
 
                                     /**
@@ -171,7 +180,8 @@ class RHS {
                                      */
     class Kink : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &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<double>  &cell_matrix,
   for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) 
     {
       const double c = (use_coefficient ?
-                       coefficient(fe_values.quadrature_point(point)) :
+                       coefficient.value(fe_values.quadrature_point(point)) :
                        1);
       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
        {
@@ -276,7 +296,7 @@ void PoissonEquation<2>::assemble (FullMatrix<double>  &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<dim>::run (ParameterHandler &prm) {
       cout << estimated_error_per_cell.l2_norm() << endl;
       estimated_error.push_back (estimated_error_per_cell.l2_norm());
 
-      Vector<double> l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_error_per_dof, estimated_error_per_dof;
+      Vector<double> l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs());
+      Vector<double> h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs());
       Vector<double> 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,
index 572a209dea2a0c9636d6060e19cbe1250e6e9da7..317a7282c971032ab984a5aafb380f07155d223e 100644 (file)
@@ -174,7 +174,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
 };
 
 
@@ -186,19 +187,24 @@ class Solution : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
                                     /**
                                      * Return the gradient of the function
                                      * at the given point.
                                      */
-    virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+    virtual Tensor<1,dim> gradient (const Point<dim> &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<double>  &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<dim>::run (const unsigned int level) {
 
   if (dof->DoFHandler<dim>::n_dofs()<=5000) 
     {
-      Vector<double> l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_seminorm_error_per_dof, h1_error_per_dof;
+      Vector<double> l1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> l2_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> linfty_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> h1_seminorm_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> h1_error_per_dof (dof->DoFHandler<dim>::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,
index dd4053a9e219ae0fe64638783514772f3c11f275..0034eb8e9fce44ab0f5106e88ef92a89bda9a49c 100644 (file)
@@ -33,7 +33,7 @@ template <int dim>
 class RightHandSide :  public Function<dim> 
 {
   public:
-    double operator () (const Point<dim> &p) const 
+    double value (const Point<dim> &p) const 
       {
        double x = 80;
        for (unsigned int d=0; d<dim; ++d)
index 09ea3ce09a759f4f84c0f94f3f259bed7339a8fd..4c48dabf0b323b577fa5a3de8f879f1c7287c590 100644 (file)
@@ -16,13 +16,32 @@ class BoundaryValuesSine : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const {
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const {
+      Assert (component==0, ExcIndexRange (component, 0, 1));
+
       double x = 1;
       
       for (unsigned int i=0; i<dim; ++i)
        x *= cos(2*3.1415926536*p(i));
       return x;
     };
+
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual void value (const Point<dim> &p,
+                       Vector<double>   &values) const {
+      Assert (values.size()==1, ExcVectorHasWrongSize (values.size(), 1));
+
+      double x = 1;
+      
+      for (unsigned int i=0; i<dim; ++i)
+       x *= cos(2*3.1415926536*p(i));
+
+      values(0) = x;
+    };
     
 
                                     /**
@@ -32,11 +51,12 @@ class BoundaryValuesSine : public Function<dim> {
                                      * empty.
                                      */
     virtual void value_list (const vector<Point<dim> > &points,
-                            vector<double>            &values) const {
+                            vector<double>            &values,
+                            const unsigned int         component) const {
       Assert (values.size() == points.size(),
              ExcVectorHasWrongSize(values.size(), points.size()));
       for (unsigned int i=0; i<points.size(); ++i) 
-       values[i] = BoundaryValuesSine<dim>::operator() (points[i]);
+       values[i] = BoundaryValuesSine<dim>::value (points[i], component);
     };
 };
 
@@ -49,7 +69,9 @@ class BoundaryValuesJump : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const {
+    virtual double value (const Point<dim>   &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<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int) const;
 };
 
 
@@ -89,7 +112,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int) const;
 };
 
 
@@ -200,7 +224,10 @@ CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_
 
 
 template <int dim>
-double RHSTrigPoly<dim>::operator () (const Point<dim> &p) const {
+double RHSTrigPoly<dim>::value (const Point<dim>   &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<dim>::operator () (const Point<dim> &p) const {
 
 
 template <int dim>
-double RHSPoly<dim>::operator () (const Point<dim> &p) const {
+double RHSPoly<dim>::value (const Point<dim>   &p,
+                           const unsigned int  component) const {
+  Assert (component==0, ExcIndexRange (component, 0, 1));
+
   double ret_val = 0;
   for (unsigned int i=0; i<dim; ++i)
     ret_val += 2*p(i)*(1.-p(i));
index 95bd7cd45ba18f7d9e4dbbc2ea60e6276b372bfe..e07c6768f9ac30ed6cb1bb905b973d5fb2b329e6 100644 (file)
@@ -55,7 +55,7 @@
  * offering iterator functions and some minor additional requirements is
  * simple.
  *
- * Note that this class could in principle be base on the C++ #map<Key,Value>#
+ * Note that this class could in principle be based on the C++ #map<Key,Value>#
  * 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.
index 3d4484f8766d06ce37207860dc924071d3538806..17b9a83cc26f87978535e98129afc268d7597a93 100644 (file)
@@ -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
  *
index 5e5b691b85979280939863207949f5efc6e96351..0b84b55ab613e48a2bc32578f368953e372bbd1f 100644 (file)
  *   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<dim># 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);
 };
 
 
index 0e798566d0443dc14526d34078f828227399d1fb..6cbb45f2ec16c8fa6393f753526591986f587ced 100644 (file)
@@ -269,36 +269,20 @@ class VectorTools
                                      */
     typedef map<unsigned char,const Function<dim>*> FunctionMap;
 
-                                    /**
-                                     * Data type for vector valued boundary function map.
-                                     */
-    typedef map<unsigned char,const VectorFunction<dim>*> 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<dim>    &dof,
                             const Function<dim>      &function,
-                            Vector<double>           &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<dim>    &dof,
-                            const VectorFunction<dim>&vectorfunction,
                             Vector<double>           &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<int,double>       &boundary_values);
 
-                                    /**
-                                     * Create boundary value information for vector
-                                     * valued functions.
-                                     * See the other #interpolate_boundary_values#.
-                                     */
-    static void interpolate_boundary_values (const DoFHandler<dim> &dof,
-                                            const VectorFunctionMap     &dirichlet_bc,
-                                            map<int,double>       &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<dim>      *weight=0);
 
-                                    /**
-                                     * Compute the error for the solution of a system.
-                                     * See the other #integrate_difference#.
-                                     */
-    static void integrate_difference (const DoFHandler<dim>   &dof,
-                                     const Vector<double>     &fe_function,
-                                     const VectorFunction<dim>&exact_solution,
-                                     Vector<float>            &difference,
-                                     const Quadrature<dim>    &q,
-                                     const NormType           &norm,
-                                     const Function<dim> *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);
 };
 
 
index a1ecfd31b6f39c610744f803d5c59a95363402e3..835c193af702fccc45c0b90af71d14fc74b8fe0e 100644 (file)
@@ -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<vector<Tensor<1,dim> > > gradients (2, vector<Tensor<1,1> >(dof.get_fe().n_components));
+         vector<vector<Tensor<1,dim> > >
+           gradients (2, vector<Tensor<1,1> >(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();
        };
index edbd3f1f2a48e71e2fb33596dd96976f5a6ca230..c409fc3dfcf902fd8f485b0bdb751b0a1f471aaf 100644 (file)
@@ -128,7 +128,9 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
                                                      vector<int>              &dof_to_boundary_mapping,
                                                      const Function<dim>      *a) {
   const FiniteElement<dim> &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<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
   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<double> cell_matrix(dofs_per_cell, dofs_per_cell);
   Vector<double>     cell_vector(dofs_per_cell);
   
   
   UpdateFlags update_flags = UpdateFlags (update_JxW_values | update_q_points);
   FEFaceValues<dim> fe_values (fe, q, update_flags);
-  
+
+                                  // two variables for the coefficient,
+                                  // one for the two cases indicated in
+                                  // the name
+  vector<double>          coefficient_values_scalar (fe_values.n_quadrature_points);
+  vector<Vector<double> > coefficient_values_system (fe_values.n_quadrature_points,
+                                                    Vector<double>(n_components));
+
+  vector<double>          rhs_values_scalar (fe_values.n_quadrature_points);
+  vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
+                                            Vector<double>(n_components));
+
+  vector<int> dofs (dofs_per_cell);
+  vector<int> dofs_on_face_vector (dofs_per_face);
+  set<int> dofs_on_face;
+
   DoFHandler<dim>::active_cell_iterator cell = dof.begin_active (),
                                        endc = dof.end ();
   for (; cell!=endc; ++cell)
@@ -167,39 +183,95 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
 
          const FullMatrix<double> &values    = fe_values.get_shape_values ();
          const vector<double>     &weights   = fe_values.get_JxW_values ();
-         vector<double>           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<double> coefficient_values (fe_values.n_quadrature_points);
-             a->value_list (fe_values.get_quadrature_points(), coefficient_values);
-             for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-               for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
-                 {
-                   for (unsigned int j=0; j<fe_values.total_dofs; ++j)
-                     cell_matrix(i,j) += (values(i,point) *
-                                          values(j,point) *
-                                          weights[point] *
-                                          coefficient_values[point]);
-                   cell_vector(i) += values(i,point) *
-                                     rhs_values[point] *
-                                     weights[point];
-                 };
+             rhs.find(cell->face(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; point<fe_values.n_quadrature_points; ++point)
+                   for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                     {
+                       for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                         if (fe.system_to_component_index(i).first ==
+                             fe.system_to_component_index(j).first)
+                           {
+                             cell_matrix(i,j)
+                               += (values(i,point) *
+                                   values(j,point) *
+                                   weights[point] *
+                                   coefficient_values_system[point](
+                                     fe.system_to_component_index(i).first));
+                           };
+                       
+                       cell_vector(i) += values(i,point) *
+                                         rhs_values_system[point](
+                                           fe.system_to_component_index(i).first) *
+                                         weights[point];
+                     };
+               }
+             else
+               for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                 for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                   {
+                     for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                       if (fe.system_to_component_index(i).first ==
+                           fe.system_to_component_index(j).first)
+                         {
+                           cell_matrix(i,j) += (values(i,point) *
+                                                values(j,point) *
+                                                weights[point]);
+                         };
+                     
+                     cell_vector(i) += values(i,point) *
+                                       rhs_values_system[point](
+                                         fe.system_to_component_index(i).first) *
+                                       weights[point];
+                   };
            }
          else
-           for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-             for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                                            // FE is a scalar one
+           {
+             rhs.find(cell->face(face)->boundary_indicator())
+               ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+             if (a != 0)
                {
-                 for (unsigned int j=0; j<fe_values.total_dofs; ++j)
-                   cell_matrix(i,j) += (values(i,point) *
-                                        values(j,point) *
-                                        weights[point]);
-                 cell_vector(i) += values(i,point) *
-                                   rhs_values[point] *
-                                   weights[point];
-               };
+                 a->value_list (fe_values.get_quadrature_points(),
+                                coefficient_values_scalar);
+                 for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                   for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                     {
+                       for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                         cell_matrix(i,j) += (values(i,point) *
+                                              values(j,point) *
+                                              weights[point] *
+                                              coefficient_values_scalar[point]);
+                       cell_vector(i) += values(i,point) *
+                                         rhs_values_scalar[point] *
+                                         weights[point];
+                     };
+               }
+             else
+               for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                 for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
+                   {
+                     for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+                       cell_matrix(i,j) += (values(i,point) *
+                                            values(j,point) *
+                                            weights[point]);
+                     cell_vector(i) += values(i,point) *
+                                       rhs_values_scalar[point] *
+                                       weights[point];
+                   };
+           };
+         
 
 
                                           // now transfer cell matrix and vector
@@ -249,13 +321,13 @@ void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
                                           // inefficient, so we copy the dofs
                                           // into a set, which enables binary
                                           // searches.
-         vector<int> dofs (dofs_per_cell);
          cell->get_dof_indices (dofs);
-
-         vector<int> dofs_on_face_vector (dofs_per_face);
          cell->face(face)->get_dof_indices (dofs_on_face_vector);
-         set<int> 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<dim>::create_boundary_mass_matrix (const DoFHandler<dim>    &
 
 
 
-
-
 template <int dim>
 void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim>    &dof,
                                                const Quadrature<dim>    &q,
@@ -593,24 +663,53 @@ void MassMatrix<dim>::assemble (FullMatrix<double>      &cell_matrix,
   
   if (coefficient != 0)
     {
-      vector<double> coefficient_values (fe_values.n_quadrature_points);
-      coefficient->value_list (fe_values.get_quadrature_points(),
-                              coefficient_values);
-      for (unsigned int i=0; i<total_dofs; ++i) 
-       for (unsigned int j=0; j<total_dofs; ++j)
-         if ((n_components == 1)
-             ||
-             (fe.system_to_component_index(i).first ==
-              fe.system_to_component_index(j).first))
-           {
-             for (unsigned int point=0; point<n_q_points; ++point)
-               cell_matrix(i,j) += (values(i,point) *
-                                    values(j,point) *
-                                    weights[point] *
-                                    coefficient_values[point]);
-           };
+      if (coefficient->n_components == 1)
+                                        // scalar coefficient given
+       {
+         vector<double> coefficient_values (fe_values.n_quadrature_points);
+         coefficient->value_list (fe_values.get_quadrature_points(),
+                                  coefficient_values);
+         for (unsigned int i=0; i<total_dofs; ++i) 
+           for (unsigned int j=0; j<total_dofs; ++j)
+             if ((n_components == 1)
+                 ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               {
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point] *
+                                        coefficient_values[point]);
+               };
+       }
+      else
+                                        // vectorial coefficient
+                                        // given
+       {
+         vector<Vector<double> > coefficient_values (fe_values.n_quadrature_points,
+                                                     Vector<double>(n_components));
+         coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                         coefficient_values);
+         for (unsigned int i=0; i<total_dofs; ++i) 
+           for (unsigned int j=0; j<total_dofs; ++j)
+             if ((n_components == 1)
+                 ||
+                 (fe.system_to_component_index(i).first ==
+                  fe.system_to_component_index(j).first))
+               {
+                 for (unsigned int point=0; point<n_q_points; ++point)
+                   cell_matrix(i,j) += (values(i,point) *
+                                        values(j,point) *
+                                        weights[point] *
+                                        coefficient_values[point](
+                                          fe.system_to_component_index(i).first));
+               };
+       };
+      
     }
   else
+                                    // no coefficient given
     for (unsigned int i=0; i<total_dofs; ++i) 
       for (unsigned int j=0; j<total_dofs; ++j)
        if ((n_components == 1)
@@ -639,8 +738,8 @@ void MassMatrix<dim>::assemble (FullMatrix<double>  &cell_matrix,
   const FiniteElement<dim>    &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<dim>::assemble (Vector<double>      &rhs,
   const FiniteElement<dim>    &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<dim>::assemble (FullMatrix<double>         &cell_matrix,
   const FiniteElement<dim>    &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<dim>::assemble (FullMatrix<double>  &cell_matrix,
   const FiniteElement<dim>    &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<dim>::assemble (Vector<double>      &rhs,
   const FiniteElement<dim>    &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,
index ab82810c21668c13ca6beeeff6f9faafab45fa81..6597ff783b7448431ae67f6aea8f9e71eadac86b 100644 (file)
@@ -4,7 +4,6 @@
 
 
 #include <base/function.h>
-#include <base/tensorfunction.h>
 #include <grid/dof.h>
 #include <grid/dof_accessor.h>
 #include <grid/tria_iterator.h>
@@ -42,56 +41,18 @@ inline double sqr_point (const Tensor<1,dim> &p) {
 
 
 
-
-
 template <int dim>
-void VectorTools<dim>::interpolate (const DoFHandler<dim>    &dof,
-                                   const Function<dim>      &function,
-                                   Vector<double>           &vec)
+void VectorTools<dim>::interpolate (const DoFHandler<dim> &dof,
+                                   const Function<dim>   &function,
+                                   Vector<double>        &vec)
 {
-  const FiniteElement<dim> &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<dim>::active_cell_iterator cell = dof.begin_active(),
-                                       endc = dof.end();
-  vector<int>         dofs_on_cell (fe.total_dofs);
-  vector<double>      dof_values_on_cell (fe.total_dofs);
-  vector<Point<dim> > 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<fe.total_dofs; ++i)
-       vec(dofs_on_cell[i]) = dof_values_on_cell[i];
-    };
-};
-
-
-
-template <int dim>
-void VectorTools<dim>::interpolate (const DoFHandler<dim>    &dof,
-                                   const VectorFunction<dim>&vectorfunction,
-                                   Vector<double>           &vec)
-{
-  const FiniteElement<dim> &fe = dof.get_fe();
+  const FiniteElement<dim> &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<dim>::active_cell_iterator cell = dof.begin_active(),
                                        endc = dof.end();
 
@@ -185,7 +146,7 @@ void VectorTools<dim>::interpolate (const DoFHandler<dim>    &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<dim>::interpolate (const DoFHandler<dim>    &dof,
   vector<Point<dim> > support_points (fe.total_dofs);
 
   vector<Point<dim> > rep_points (n_rep_points);
-  vector<Vector<double> > function_values_at_rep_points (
-    n_rep_points, Vector<double>(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<double>          function_values_scalar (n_rep_points);
+  vector<Vector<double> > function_values_system (n_rep_points,
+                                                 Vector<double>(fe.n_components));
 
   for (; cell!=endc; ++cell)
     {
@@ -263,26 +232,48 @@ void VectorTools<dim>::interpolate (const DoFHandler<dim>    &dof,
       for (unsigned int j=0; j<dofs_of_rep_points.size(); ++j)
        rep_points[j]=support_points[dofs_of_rep_points[j]];
 
-                                      // get function values at these points
-      vectorfunction.value_list (rep_points, function_values_at_rep_points);
-
                                       // get indices of the dofs on this cell
       cell->get_dof_indices (dofs_on_cell);
 
-                                      // distribute the function values to
-                                      // the global vector
-      for (unsigned int i=0; i<fe.total_dofs; ++i)
+
+                                      
+      if (fe_is_system)
        {
-         const unsigned int component
-           = fe.system_to_component_index(i).first;
-         const unsigned int rep_dof=dof_to_rep_index_table[i];
-         vec(dofs_on_cell[i])
-           = function_values_at_rep_points[rep_dof](component);
+                                          // get function values at
+                                          // these points. Here: get
+                                          // all components
+         function.vector_value_list (rep_points, function_values_system);
+                                          // distribute the function
+                                          // values to the global
+                                          // vector
+         for (unsigned int i=0; i<fe.total_dofs; ++i)
+           {
+             const unsigned int component
+               = fe.system_to_component_index(i).first;
+             const unsigned int rep_dof=dof_to_rep_index_table[i];
+             vec(dofs_on_cell[i])
+               = function_values_system[rep_dof](component);
+           };
        }
+      
+      else
+       {
+                                          // get first component only,
+                                          // which is the only component
+                                          // in the function anyway
+         function.value_list (rep_points, function_values_scalar, 0);
+                                          // distribute the function
+                                          // values to the global
+                                          // vector
+         for (unsigned int i=0; i<fe.total_dofs; ++i)
+           vec(dofs_on_cell[i]) 
+             = function_values_scalar[dof_to_rep_index_table[i]];
+       };
     }
 }
 
 
+
 template <int dim> void
 VectorTools<dim>::interpolate(const DoFHandler<dim>           &high_dof,
                              const DoFHandler<dim>           &low_dof,
@@ -342,7 +333,8 @@ void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
                                const Quadrature<dim-1>  &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<dim> &fe = dof.get_fe();
 
@@ -364,6 +356,10 @@ void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
            for (unsigned int i=0; i<fe.dofs_per_face; ++i)
                                               // enter zero boundary values
                                               // for all boundary nodes
+                                              //
+                                              // we need not care about
+                                              // vector valued elements here,
+                                              // since we set all components
              boundary_values[face_dof_indices[i]] = 0.;
          };
     }
@@ -412,8 +408,9 @@ void VectorTools<dim>::project (const DoFHandler<dim>    &dof,
 
   constraints.condense (mass_matrix);
   constraints.condense (tmp);
-  MatrixTools<dim>::apply_boundary_values (boundary_values,
-                                          mass_matrix, vec, tmp);
+  if (boundary_values.size() != 0)
+    MatrixTools<dim>::apply_boundary_values (boundary_values,
+                                            mass_matrix, vec, tmp);
 
   SolverControl                    control(1000,1e-16);
   PrimitiveVectorMemory<Vector<double> >   memory;
@@ -438,7 +435,8 @@ void VectorTools<dim>::create_right_hand_side (const DoFHandler<dim>    &dof,
                                               const Function<dim>      &rhs,
                                               Vector<double>           &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<dim>::create_right_hand_side (const DoFHandler<dim>    &dof,
 template <>
 void
 VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &dof,
-                                            const FunctionMap &dirichlet_bc,
-                                            map<int,double>   &boundary_values)
+                                            const FunctionMap   &dirichlet_bc,
+                                            map<int,double>     &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; i<fe.dofs_per_vertex; ++i)
+       boundary_values[leftmost_cell->vertex_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; i<fe.dofs_per_vertex; ++i)
+       boundary_values[rightmost_cell->vertex_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<int,double>&)
-{
-  Assert (false, ExcNotImplemented());
-};
-
 #endif
 
 
@@ -540,64 +533,32 @@ template <int dim>
 void
 VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
                                               const FunctionMap     &dirichlet_bc,
-                                              map<int,double>   &boundary_values) {
-  Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
-         ExcInvalidBoundaryIndicator());
-
-  const FiniteElement<dim> &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<int>         face_dofs (fe.dofs_per_face);
-  vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
-  vector<double>      dof_values (fe.dofs_per_face);
-       
-  DoFHandler<dim>::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<face_dofs.size(); ++i)
-         boundary_values[face_dofs[i]] = dof_values[i];
-      };
-};
-
-
-
-template <int dim>
-void
-VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
-                                              const VectorFunctionMap     &dirichlet_bc,
-                                              map<int,double>   &boundary_values)
+                                              map<int,double>       &boundary_values)
 {
   Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
          ExcInvalidBoundaryIndicator());
 
-  const FiniteElement<dim> &fe = dof.get_fe();
-  Assert (fe.n_components == dirichlet_bc.begin()->second->n_components,
+  const FiniteElement<dim> &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<int>         face_dofs (fe.dofs_per_face, -1);
   vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
-  vector< Vector<double> > dof_values (fe.dofs_per_face,
-                                      Vector<double>(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<double>          dof_values_scalar (fe.dofs_per_face);
+  vector<Vector<double> > dof_values_system (fe.dofs_per_face,
+                                            Vector<double>(fe.n_components));
        
   DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
                                        endf = dof.end_face();
@@ -612,13 +573,32 @@ VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &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; i<face_dofs.size(); ++i)
-         boundary_values[face_dofs[i]]
-           = dof_values[i](fe.face_system_to_component_index(i).first);
+       if (fe_is_system)
+         {
+           function_ptr->second->vector_value_list (dof_locations, dof_values_system);
+           
+                                            // enter into list
+           
+           for (unsigned int i=0; i<face_dofs.size(); ++i)
+             boundary_values[face_dofs[i]]
+               = dof_values_system[i](fe.face_system_to_component_index(i).first);
+         }
+       else
+                                          // fe has only one component,
+                                          // so save some computations
+         {
+                                            // get only the one component that
+                                            // this function has
+           function_ptr->second->value_list (dof_locations,
+                                             dof_values_scalar,
+                                             0);
+           
+                                            // enter into list
+           
+           for (unsigned int i=0; i<face_dofs.size(); ++i)
+             boundary_values[face_dofs[i]] = dof_values_scalar[i];
+         };
       };
 }
 
@@ -630,7 +610,8 @@ VectorTools<dim>::project_boundary_values (const DoFHandler<dim>    &dof,
                                           const FunctionMap        &boundary_functions,
                                           const Quadrature<dim-1>  &q,
                                           map<int,double>          &boundary_values) {
-  Assert (dof.get_fe().n_components == 1, ExcInvalidFE());
+  Assert (dof.get_fe().n_components == boundary_functions.begin()->second->n_components,
+         ExcComponentMismatch());
   
   vector<int>    dof_to_boundary_mapping;
   dof.map_dof_to_boundary_indices (boundary_functions, dof_to_boundary_mapping);
@@ -703,233 +684,24 @@ VectorTools<dim>::project_boundary_values (const DoFHandler<dim>    &dof,
 
 
 
-
-template <int dim>
-void VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
-                                            const Vector<double>     &fe_function,
-                                            const Function<dim>      &exact_solution,
-                                            Vector<float>            &difference,
-                                            const Quadrature<dim>    &q,
-                                            const NormType           &norm,
-                                            const Function<dim>      *weight=0)
-{
-  const FiniteElement<dim> &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<dim> fe_values(fe, q, update_flags);
-  
-                                  // loop over all cells
-  DoFHandler<dim>::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<double>   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<double> 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<double>());
-             };            
-
-                                            // 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<double> w(n_q_points);
-               weight->value_list(fe_values.get_quadrature_points(),w);
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 psi[q]*=w[q];
-             }
-
-                                            // ok, now we have the integrand,
-                                            // let's compute the integral,
-                                            // which is
-                                            // sum_j psi_j JxW_j
-                                            // (or |psi_j| or |psi_j|^2
-           switch (norm) 
-             {
-               case mean:
-               case L1_norm:
-                     diff = inner_product (psi.begin(), psi.end(),
-                                           fe_values.get_JxW_values().begin(),
-                                           0.0);
-                     break;
-               case L2_norm:
-               case H1_norm:
-                     diff = sqrt(inner_product (psi.begin(), psi.end(),
-                                                fe_values.get_JxW_values().begin(),
-                                                0.0));
-                     break;
-               case Linfty_norm:
-                     diff = *max_element (psi.begin(), psi.end());
-                     break;
-               default:
-                     Assert (false, ExcNotImplemented());
-             };
-
-                                            // note: the H1_norm uses the result
-                                            // of the L2_norm and control goes
-                                            // over to the next case statement!
-           if (norm != H1_norm)
-             break;
-         };
-
-         case H1_seminorm:
-         {
-                                            // note: the computation of the
-                                            // H1_norm starts at the previous
-                                            // case statement, but continues
-                                            // here!
-
-                                            // for H1_norm: re-square L2_norm.
-           diff = sqr(diff);
-
-                                            // same procedure as above, but now
-                                            // psi is a vector of gradients
-           const unsigned int n_q_points = q.n_quadrature_points;
-           vector<Tensor<1,dim> >   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<Tensor<1,dim> > 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<Tensor<1,dim> >());
-             };
-                                            // take square of integrand
-           vector<double> psi_square (psi.size(), 0.0);
-           for (unsigned int i=0; i<n_q_points; ++i)
-             psi_square[i] = sqr_point(psi[i]);
-
-                                            // now weight the values
-                                            // at the quadrature points due
-                                            // to the weighting function
-           if (weight)
-             {
-               vector<double> w(n_q_points);
-               weight->value_list(fe_values.get_quadrature_points(),w);
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 psi_square[q]*=w[q];
-             }
-
-                                            // add seminorm to L_2 norm or
-                                            // to zero
-           diff += inner_product (psi_square.begin(), psi_square.end(),
-                                  fe_values.get_JxW_values().begin(),
-                                  0.0);
-           diff = sqrt(diff);
-
-           break;
-         };
-                                            
-         default:
-               Assert (false, ExcNotImplemented());
-       };
-
-      
-                                      // append result of this cell
-                                      // to the end of the vector
-      difference(index) = diff;
-    };
-};
-
-
-
 template <int dim>
 void
 VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                                        const Vector<double>     &fe_function,
-                                       const VectorFunction<dim>&exact_solution,
+                                       const Function<dim>      &exact_solution,
                                        Vector<float>            &difference,
                                        const Quadrature<dim>    &q,
                                        const NormType           &norm,
                                        const Function<dim>      *weight)
 {
-  Assert(norm != mean , ExcNotUseful());
+  const unsigned int        n_q_points   = q.n_quadrature_points;
+  const FiniteElement<dim> &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<dim> &fe = dof.get_fe();
-  
   difference.reinit (dof.get_tria().n_active_cells());
   
   UpdateFlags update_flags = UpdateFlags (update_q_points  |
@@ -937,6 +709,25 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
   if ((norm==H1_seminorm) || (norm==H1_norm))
     update_flags = UpdateFlags (update_flags | update_gradients);
   FEValues<dim> fe_values(fe, q, update_flags);
+
+  vector< Vector<double> >        function_values (n_q_points,
+                                                  Vector<double>(n_components));
+  vector<vector<Tensor<1,dim> > > function_grads (n_q_points,
+                                                 vector<Tensor<1,dim> >(n_components));
+  vector<double> weight_values (n_q_points);
+  
+  vector<Vector<double> >         psi_values (n_q_points,
+                                             Vector<double>(n_components));
+  vector<vector<Tensor<1,dim> > > psi_grads (n_q_points,
+                                            vector<Tensor<1,dim> >(n_components));
+  vector<double> psi_scalar (n_q_points);
+  vector<double> psi_square (n_q_points);
+           
+                                  // tmp vector when we use the
+                                  // Function<dim> functions for
+                                  // scalar functions
+  vector<double>         tmp_values (fe_values.n_quadrature_points);
+  vector<Tensor<1,dim> > tmp_gradients (fe_values.n_quadrature_points);
   
                                   // loop over all cells
   DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
@@ -950,30 +741,34 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &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<Vector<double> >  psi (n_q_points, Vector<double>(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<n_q_points; ++i)
+                 psi_values[i](0) = tmp_values[i];
+             };
+           
                                             // then subtract finite element
                                             // fe_function
-           if (true) 
-             {
-               vector< Vector<double> > function_values (
-                 n_q_points, Vector<double>(fe.n_components));
-
-               fe_values.get_function_values (fe_function, function_values);
-
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 psi[q] -= function_values[q];
-             };            
+           fe_values.get_function_values (fe_function, function_values);
+           for (unsigned int q=0; q<n_q_points; ++q)
+             psi_values[q] -= function_values[q];
 
                                             // for L1_norm, Linfty_norm, L2_norm
                                             // and H1_norm take square of the
@@ -984,7 +779,6 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                                             // Use psi_scalar to store the squares
                                             // of the vectors or the vector norms
                                             // respectively.
-           vector<double>  psi_scalar (n_q_points);
            switch (norm)
              {
                case mean:
@@ -994,7 +788,7 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                case L2_norm:
                case H1_norm:
                      for (unsigned int q=0; q<n_q_points; ++q)
-                       psi_scalar[q]=psi[q].norm_sqr();
+                       psi_scalar[q] = psi_values[q].norm_sqr();
                      
                      if (norm == L1_norm || norm == Linfty_norm)
                        transform (psi_scalar.begin(), psi_scalar.end(),
@@ -1009,10 +803,10 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
                                             // to the weighting function
            if (weight)
              {
-               vector<double> 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<n_q_points; ++q)
-                 psi_scalar[q]*=w[q];
+                 psi_scalar[q] *= weight_values[q];
              }
 
                                             // ok, now we have the integrand,
@@ -1023,7 +817,6 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &dof,
            switch (norm)
              {
                case mean:
-                     break;      
                case L1_norm:
                case L2_norm:
                case H1_norm:
@@ -1057,45 +850,48 @@ VectorTools<dim>::integrate_difference (const DoFHandler<dim>    &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<vector<Tensor<1,dim> > >   psi (
-             n_q_points, vector<Tensor<1,dim> >(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<n_q_points; ++i)
+                 psi_grads[i][0] = tmp_gradients[i];
+             };
 
                                             // then subtract finite element
                                             // function_grads
-           if (true) 
-             {
-               vector<vector<Tensor<1,dim> > > function_grads (
-                 n_q_points, vector<Tensor<1,dim> >(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<n_components; ++k)
+             for (unsigned int q=0; q<n_q_points; ++q)
+               psi_grads[q][k] -= function_grads[q][k];
 
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 for (unsigned int k=0; k<fe.n_components; ++k)
-                   psi[q][k] -= function_grads[q][k];
-             };
-                                            // take square of integrand
-           vector<double> psi_square (psi.size(), 0.0);
-           for (unsigned int q=0; q<n_q_points; ++q)
-             for (unsigned int k=0; k<fe.n_components; ++k)
-               psi_square[q] += sqr_point(psi[q][k]);
+                                            // take square of integrand
+           fill_n (psi_square.begin(), n_q_points, 0.0);
+           for (unsigned int k=0; k<n_components; ++k)
+             for (unsigned int q=0; q<n_q_points; ++q)
+               psi_square[q] += sqr_point(psi_grads[q][k]);
 
                                             // now weight the values
                                             // at the quadrature points due
                                             // to the weighting function
            if (weight)
              {
-               vector<double> 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<n_q_points; ++q)
-                 psi_square[q]*=w[q];
+                 psi_square[q] *= weight_values[q];
              }
 
                                             // add seminorm to L_2 norm or
index 1a5ad43c5c4f5166c420d38dd692b7a93b528fb8..afe9ac1f4b4fae4deb4d485e251169ceffb96905 100644 (file)
@@ -94,7 +94,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
 };
 
 
@@ -106,19 +107,22 @@ class Solution : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
                                     /**
                                      * Return the gradient of the function
                                      * at the given point.
                                      */
-    virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+    virtual Tensor<1,dim> gradient (const Point<dim> &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<double>  &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<dim>::run (const unsigned int level) {
 
   if (dof->n_dofs()<=5000) 
     {
-      Vector<double> l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_seminorm_error_per_dof, h1_error_per_dof;
+      Vector<double> l1_error_per_dof(dof->n_dofs());
+      Vector<double> l2_error_per_dof(dof->n_dofs());
+      Vector<double> linfty_error_per_dof(dof->n_dofs());
+      Vector<double> h1_seminorm_error_per_dof(dof->n_dofs());
+      Vector<double> 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,
index 5158c8c901bf0015c0ab5bc20c5ff309c4e1ac71..aebd9cf99929204c72355474039bad510cff107c 100644 (file)
@@ -114,25 +114,32 @@ class Solution {
 
     class GaussShape : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 
     class Singular : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 
     class Kink : public Function<dim> {
       public:
        class Coefficient : public Function<dim> {
          public:
-           virtual double operator () (const Point<dim> &p) const;
+           virtual double value (const Point<dim> &p,
+                                 const unsigned int component) const;
        };
        
-       virtual double operator () (const Point<dim> &p) const;
-       virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
+       virtual Tensor<1,dim> gradient (const Point<dim> &p,
+                                       const unsigned int component) const;
     };
 };
 
@@ -150,7 +157,8 @@ class RHS {
                                      */
     class GaussShape : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
     };
 
                                     /**
@@ -160,7 +168,8 @@ class RHS {
                                      */
     class Singular : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &p,
+                             const unsigned int component) const;
     };
 
                                     /**
@@ -171,7 +180,8 @@ class RHS {
                                      */
     class Kink : public Function<dim> {
       public:
-       virtual double operator () (const Point<dim> &p) const;
+       virtual double value (const Point<dim> &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<double>  &cell_matrix,
   for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point) 
     {
       const double c = (use_coefficient ?
-                       coefficient(fe_values.quadrature_point(point)) :
+                       coefficient.value(fe_values.quadrature_point(point)) :
                        1);
       for (unsigned int i=0; i<fe_values.total_dofs; ++i) 
        {
@@ -276,7 +296,7 @@ void PoissonEquation<2>::assemble (FullMatrix<double>  &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<dim>::run (ParameterHandler &prm) {
       cout << estimated_error_per_cell.l2_norm() << endl;
       estimated_error.push_back (estimated_error_per_cell.l2_norm());
 
-      Vector<double> l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_error_per_dof, estimated_error_per_dof;
+      Vector<double> l2_error_per_dof(dof->n_dofs()), linfty_error_per_dof(dof->n_dofs());
+      Vector<double> h1_error_per_dof(dof->n_dofs()), estimated_error_per_dof(dof->n_dofs());
       Vector<double> 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,
index 572a209dea2a0c9636d6060e19cbe1250e6e9da7..317a7282c971032ab984a5aafb380f07155d223e 100644 (file)
@@ -174,7 +174,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
 };
 
 
@@ -186,19 +187,24 @@ class Solution : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int component) const;
                                     /**
                                      * Return the gradient of the function
                                      * at the given point.
                                      */
-    virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+    virtual Tensor<1,dim> gradient (const Point<dim> &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<double>  &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<dim>::run (const unsigned int level) {
 
   if (dof->DoFHandler<dim>::n_dofs()<=5000) 
     {
-      Vector<double> l1_error_per_dof, l2_error_per_dof, linfty_error_per_dof;
-      Vector<double> h1_seminorm_error_per_dof, h1_error_per_dof;
+      Vector<double> l1_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> l2_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> linfty_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> h1_seminorm_error_per_dof (dof->DoFHandler<dim>::n_dofs());
+      Vector<double> h1_error_per_dof (dof->DoFHandler<dim>::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,
index dd4053a9e219ae0fe64638783514772f3c11f275..0034eb8e9fce44ab0f5106e88ef92a89bda9a49c 100644 (file)
@@ -33,7 +33,7 @@ template <int dim>
 class RightHandSide :  public Function<dim> 
 {
   public:
-    double operator () (const Point<dim> &p) const 
+    double value (const Point<dim> &p) const 
       {
        double x = 80;
        for (unsigned int d=0; d<dim; ++d)
index 09ea3ce09a759f4f84c0f94f3f259bed7339a8fd..4c48dabf0b323b577fa5a3de8f879f1c7287c590 100644 (file)
@@ -16,13 +16,32 @@ class BoundaryValuesSine : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const {
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const {
+      Assert (component==0, ExcIndexRange (component, 0, 1));
+
       double x = 1;
       
       for (unsigned int i=0; i<dim; ++i)
        x *= cos(2*3.1415926536*p(i));
       return x;
     };
+
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual void value (const Point<dim> &p,
+                       Vector<double>   &values) const {
+      Assert (values.size()==1, ExcVectorHasWrongSize (values.size(), 1));
+
+      double x = 1;
+      
+      for (unsigned int i=0; i<dim; ++i)
+       x *= cos(2*3.1415926536*p(i));
+
+      values(0) = x;
+    };
     
 
                                     /**
@@ -32,11 +51,12 @@ class BoundaryValuesSine : public Function<dim> {
                                      * empty.
                                      */
     virtual void value_list (const vector<Point<dim> > &points,
-                            vector<double>            &values) const {
+                            vector<double>            &values,
+                            const unsigned int         component) const {
       Assert (values.size() == points.size(),
              ExcVectorHasWrongSize(values.size(), points.size()));
       for (unsigned int i=0; i<points.size(); ++i) 
-       values[i] = BoundaryValuesSine<dim>::operator() (points[i]);
+       values[i] = BoundaryValuesSine<dim>::value (points[i], component);
     };
 };
 
@@ -49,7 +69,9 @@ class BoundaryValuesJump : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const {
+    virtual double value (const Point<dim>   &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<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int) const;
 };
 
 
@@ -89,7 +112,8 @@ class RHSPoly : public Function<dim> {
                                      * Return the value of the function
                                      * at the given point.
                                      */
-    virtual double operator () (const Point<dim> &p) const;
+    virtual double value (const Point<dim> &p,
+                         const unsigned int) const;
 };
 
 
@@ -200,7 +224,10 @@ CurvedLine<dim>::get_new_point_on_quad (const typename Triangulation<dim>::quad_
 
 
 template <int dim>
-double RHSTrigPoly<dim>::operator () (const Point<dim> &p) const {
+double RHSTrigPoly<dim>::value (const Point<dim>   &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<dim>::operator () (const Point<dim> &p) const {
 
 
 template <int dim>
-double RHSPoly<dim>::operator () (const Point<dim> &p) const {
+double RHSPoly<dim>::value (const Point<dim>   &p,
+                           const unsigned int  component) const {
+  Assert (component==0, ExcIndexRange (component, 0, 1));
+
   double ret_val = 0;
   for (unsigned int i=0; i<dim; ++i)
     ret_val += 2*p(i)*(1.-p(i));

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.