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.
/*---------------------------- 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
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
*/
/**
- * 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;
};
/**
- * 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:
/**
#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.");
-
-};
-
/**
* @author Guido Kanschat, 1999
*/
template <int rank, int dim>
-class TensorFunction : public VectorFunction<dim>
+class TensorFunction : public FunctionTime,
+ public Subscriptor
{
public:
/**
* 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
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;
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
*/
#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;
}
- /**
- * 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;
* 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>
{
/**
* 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;
* 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>
{
/**
* Constructor taking #rank# indices.
*/
- TensorIndex(unsigned int i0)
+ TensorIndex(const unsigned int i0)
{
index[0] = i0;
}
* 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");
};
#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)
{};
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>
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()));
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>
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()));
+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>;
#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
template <int rank, int dim>
TensorFunction<rank, dim>::TensorFunction (const double initial_time)
:
- VectorFunction<dim>(pow(dim,rank), initial_time)
+ FunctionTime (initial_time)
{};
-// 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
};
-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>;
* 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 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;
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;
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;
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);
};
};
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,
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;
};
};
*/
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;
};
/**
*/
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;
};
/**
*/
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;
};
};
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()));
};
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;
};
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));
};
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;
};
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)
{
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);
};
};
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,
* 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 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;
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;
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;
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);
};
};
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,
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)
* 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;
+ };
/**
* 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);
};
};
* 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:
* 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;
};
* 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;
};
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)
{
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));
* 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.
* 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
*
* 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
* Exception
*/
DeclException0 (ExcInvalidFE);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcComponentMismatch);
};
*/
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);
/**
* 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.
*/
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
* #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.
*/
* 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.
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'
* Exception
*/
DeclException0 (ExcNotUseful);
-
/**
* Exception
*/
DeclException0 (ExcInvalidFE);
-
/**
* Exception
*/
DeclException0 (ExcInvalidBoundaryIndicator);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcComponentMismatch);
};
{
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
// 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);
}
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();
};
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());
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)
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
// 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
-
-
template <int dim>
void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
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)
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,
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,
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,
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,
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,
#include <base/function.h>
-#include <base/tensorfunction.h>
#include <grid/dof.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-
-
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();
// }
// 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.
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)
{
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,
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();
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.;
};
}
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;
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);
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
// 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
// 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
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();
// 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];
+ };
};
}
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);
-
-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 |
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(),
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
// 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:
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(),
// 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,
switch (norm)
{
case mean:
- break;
case L1_norm:
case L2_norm:
case H1_norm:
// 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
* 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 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;
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;
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;
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);
};
};
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,
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;
};
};
*/
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;
};
/**
*/
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;
};
/**
*/
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;
};
};
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()));
};
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;
};
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));
};
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;
};
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)
{
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);
};
};
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,
* 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 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;
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;
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;
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);
};
};
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,
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)
* 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;
+ };
/**
* 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);
};
};
* 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:
* 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;
};
* 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;
};
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)
{
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));