+/**
+ * Namespace implementing some concrete classes derived from the
+ * @ref{Function} class that describe actual functions. This is rather
+ * a collection of classes that we have needed for our own programs
+ * once and thought they might be useful to others as well at some
+ * point.
+ */
+namespace Functions
+{
+
+
/**
* The distance to the origin squared.
*
*
* @author: Guido Kanschat, 1999
*/
-template<int dim>
-class SquareFunction : public Function<dim>
-{
- public:
- /**
- * Function value at one point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Function values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at one point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at one point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
-
+ template<int dim>
+ class SquareFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Function value at one point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Function values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at one point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at one point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
+
/**
* The function @p{xy}. This function serves as an example for
* a vanishing Laplacian.
*
* @author: Guido Kanschat, 2000
*/
-template<int dim>
-class Q1WedgeFunction : public Function<dim>
-{
- public:
- /**
- * Function value at one point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Function values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at one point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at one point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
-
+ template<int dim>
+ class Q1WedgeFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Function value at one point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Function values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at one point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at one point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
+
/**
* d-quadratic pillow on the unit hypercube.
*
*
* @author: Guido Kanschat, 1999
*/
-template<int dim>
-class PillowFunction : public Function<dim>
-{
- public:
- /**
- * Constructor. Provide a
- * constant that will be added to
- * each function value.
- */
- PillowFunction (const double offset=0.);
-
- /**
- * The value at a single point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at a single point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
- private:
- const double offset;
-};
-
-
-
+ template<int dim>
+ class PillowFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Constructor. Provide a
+ * constant that will be added to
+ * each function value.
+ */
+ PillowFunction (const double offset=0.);
+
+ /**
+ * The value at a single point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ private:
+ const double offset;
+ };
+
+
+
/**
* Cosine-shaped pillow function.
* This is another function with zero boundary values on $[-1,1]^d$. In the interior
* it is the product of $\cos(\pi/2 x_i)$.
* @author Guido Kanschat, 1999
*/
-template<int dim>
-class CosineFunction : public Function<dim>
-{
- public:
- /**
- * The value at a single point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at a single point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<2,dim> hessian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void hessian_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<2,dim> > &hessians,
+ template<int dim>
+ class CosineFunction : public Function<dim>
+ {
+ public:
+ /**
+ * The value at a single point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
const unsigned int component = 0) const;
-};
-
-
-
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<2,dim> hessian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void hessian_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<2,dim> > &hessians,
+ const unsigned int component = 0) const;
+ };
+
+
+
/**
* Product of exponential functions in each coordinate direction.
* @author Guido Kanschat, 1999
*/
-template<int dim>
-class ExpFunction : public Function<dim>
-{
- public:
- /**
- * The value at a single point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at a single point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
-
+ template<int dim>
+ class ExpFunction : public Function<dim>
+ {
+ public:
+ /**
+ * The value at a single point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
+
/**
* Singularity on the L-shaped domain in 2D.
*
*
* @author Guido Kanschat, 1999
*/
-class LSingularityFunction : public Function<2>
-{
- public:
- /**
- * The value at a single point.
- */
- virtual double value (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Values at multiple points.
- */
- virtual void value_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<1,2> gradient (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at a single point.
- */
- virtual double laplacian (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at multiple points.
- */
- virtual void laplacian_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
-
+ class LSingularityFunction : public Function<2>
+ {
+ public:
+ /**
+ * The value at a single point.
+ */
+ virtual double value (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Values at multiple points.
+ */
+ virtual void value_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<1,2> gradient (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void gradient_list (const std::vector<Point<2> > &points,
+ std::vector<Tensor<1,2> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ virtual double laplacian (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ virtual void laplacian_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
+
/**
* Singularity on the slit domain in 2D.
*
*
* @author Guido Kanschat, 1999
*/
-class SlitSingularityFunction : public Function<2>
-{
- public:
- /**
- * The value at a single point.
- */
- virtual double value (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Values at multiple points.
- */
- virtual void value_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at a single point.
- */
- virtual Tensor<1,2> gradient (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Gradients at multiple points.
- */
- virtual void gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at a single point.
- */
- virtual double laplacian (const Point<2> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian at multiple points.
- */
- virtual void laplacian_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-};
-
-
+ class SlitSingularityFunction : public Function<2>
+ {
+ public:
+ /**
+ * The value at a single point.
+ */
+ virtual double value (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Values at multiple points.
+ */
+ virtual void value_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at a single point.
+ */
+ virtual Tensor<1,2> gradient (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradients at multiple points.
+ */
+ virtual void gradient_list (const std::vector<Point<2> > &points,
+ std::vector<Tensor<1,2> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ virtual double laplacian (const Point<2> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ virtual void laplacian_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+ };
+
+
/**
* A jump in x-direction transported into some direction.
*
*
* @author: Guido Kanschat, 2000
*/
-template<int dim>
-class JumpFunction : public Function<dim>
-{
- public:
- /**
- * Constructor. Provide the
- * advection direction here and
- * the steepness of the slope.
- */
- JumpFunction (const Point<dim> &direction,
- const double steepness);
-
- /**
- * Function value at one point.
- */
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Function values at multiple points.
- */
- virtual void value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Gradient at one point.
- */
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- Gradients at multiple points.
- */
- virtual void gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at one point.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- /**
- * Laplacian of the function at multiple points.
- */
- virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int component = 0) const;
-
- /**
- * Determine an estimate for
- * the memory consumption (in
- * bytes) of this
- * object. Since sometimes
- * the size of objects can
- * not be determined exactly
- * (for example: what is the
- * memory consumption of an
- * STL @p{std::map} type with a
- * certain number of
- * elements?), this is only
- * an estimate. however often
- * quite close to the true
- * value.
- */
- unsigned int memory_consumption () const;
-
- protected:
- /**
- * Advection vector.
- */
- const Point<dim> direction;
-
- /**
- * Steepness (maximal derivative)
- * of the slope.
- */
- const double steepness;
-
- /**
- * Advection angle.
- */
- double angle;
-
- /**
- * Sine of @p{angle}.
- */
- double sine;
-
- /**
- * Cosine of @p{angle}.
- */
- double cosine;
-};
-
-
-
+ template<int dim>
+ class JumpFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Constructor. Provide the
+ * advection direction here and
+ * the steepness of the slope.
+ */
+ JumpFunction (const Point<dim> &direction,
+ const double steepness);
+
+ /**
+ * Function value at one point.
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Function values at multiple points.
+ */
+ virtual void value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Gradient at one point.
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ Gradients at multiple points.
+ */
+ virtual void gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at one point.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian of the function at multiple points.
+ */
+ virtual void laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ /**
+ * Determine an estimate for
+ * the memory consumption (in
+ * bytes) of this
+ * object. Since sometimes
+ * the size of objects can
+ * not be determined exactly
+ * (for example: what is the
+ * memory consumption of an
+ * STL @p{std::map} type with a
+ * certain number of
+ * elements?), this is only
+ * an estimate. however often
+ * quite close to the true
+ * value.
+ */
+ unsigned int memory_consumption () const;
+
+ protected:
+ /**
+ * Advection vector.
+ */
+ const Point<dim> direction;
+
+ /**
+ * Steepness (maximal derivative)
+ * of the slope.
+ */
+ const double steepness;
+
+ /**
+ * Advection angle.
+ */
+ double angle;
+
+ /**
+ * Sine of @p{angle}.
+ */
+ double sine;
+
+ /**
+ * Cosine of @p{angle}.
+ */
+ double cosine;
+ };
+
+
+
/**
* Given a wavenumber vector generate a cosine function. The
* wavenumber coefficient is given as a @p{d}-dimensional point @p{k}
*
* @author Wolfgang Bangerth, 2001
*/
-template <int dim>
-class FourierCosineFunction : public Function<dim>
-{
- public:
- /**
- * Constructor. Take the Fourier
- * coefficients in each space
- * direction as argument.
- */
- FourierCosineFunction (const Point<dim> &fourier_coefficients);
-
- /**
- * 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 value (const Point<dim> &p,
- const unsigned int component = 0) 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;
-
- /**
- * Compute the Laplacian of a
- * given component at point @p{p}.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
- private:
- /**
- * Stored Fourier coefficients.
- */
- const Point<dim> fourier_coefficients;
-};
-
-
-
+ template <int dim>
+ class FourierCosineFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Constructor. Take the Fourier
+ * coefficients in each space
+ * direction as argument.
+ */
+ FourierCosineFunction (const Point<dim> &fourier_coefficients);
+
+ /**
+ * 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 value (const Point<dim> &p,
+ const unsigned int component = 0) 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;
+
+ /**
+ * Compute the Laplacian of a
+ * given component at point @p{p}.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ private:
+ /**
+ * Stored Fourier coefficients.
+ */
+ const Point<dim> fourier_coefficients;
+ };
+
+
+
/**
* Given a wavenumber vector generate a sine function. The
* wavenumber coefficient is given as a @p{d}-dimensional point @p{k}
*
* @author Wolfgang Bangerth, 2001
*/
-template <int dim>
-class FourierSineFunction : public Function<dim>
-{
- public:
- /**
- * Constructor. Take the Fourier
- * coefficients in each space
- * direction as argument.
- */
- FourierSineFunction (const Point<dim> &fourier_coefficients);
-
- /**
- * 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 value (const Point<dim> &p,
- const unsigned int component = 0) 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;
-
- /**
- * Compute the Laplacian of a
- * given component at point @p{p}.
- */
- virtual double laplacian (const Point<dim> &p,
- const unsigned int component = 0) const;
- private:
- /**
- * Stored Fourier coefficients.
- */
- const Point<dim> fourier_coefficients;
+ template <int dim>
+ class FourierSineFunction : public Function<dim>
+ {
+ public:
+ /**
+ * Constructor. Take the Fourier
+ * coefficients in each space
+ * direction as argument.
+ */
+ FourierSineFunction (const Point<dim> &fourier_coefficients);
+
+ /**
+ * 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 value (const Point<dim> &p,
+ const unsigned int component = 0) 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;
+
+ /**
+ * Compute the Laplacian of a
+ * given component at point @p{p}.
+ */
+ virtual double laplacian (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ private:
+ /**
+ * Stored Fourier coefficients.
+ */
+ const Point<dim> fourier_coefficients;
+ };
+
+
};
-
-
#endif
-template<int dim>
-double
-SquareFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
+namespace Functions
{
- return p.square();
-}
-
-
+
+
+ template<int dim>
+ double
+ SquareFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ return p.square();
+ }
+
+
// if necessary try to work around a bug in the IBM xlC compiler
#ifdef XLC_WORK_AROUND_STD_BUG
-using namespace std;
+ using namespace std;
#endif
-
-template<int dim>
-void
-SquareFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+
+ template<int dim>
+ void
+ SquareFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ values[i] = p.square();
+ }
+ }
+
+
+ template<int dim>
+ double
+ SquareFunction<dim>::laplacian (const Point<dim> &,
+ const unsigned int) const
+ {
+ return 2*dim;
+ }
+
+
+ template<int dim>
+ void
+ SquareFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ values[i] = 2*dim;
+ }
+
+
+
+ template<int dim>
+ Tensor<1,dim>
+ SquareFunction<dim>::gradient (const Point<dim> &p,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- values[i] = p.square();
- }
-}
-
-
-template<int dim>
-double
-SquareFunction<dim>::laplacian (const Point<dim> &,
- const unsigned int) const
-{
- return 2*dim;
-}
-
-
-template<int dim>
-void
-SquareFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- values[i] = 2*dim;
-}
-
-
-
-template<int dim>
-Tensor<1,dim>
-SquareFunction<dim>::gradient (const Point<dim> &p,
+ {
+ return p*2;
+ }
+
+
+
+ template<int dim>
+ void
+ SquareFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ gradients[i] = points[i]*2;
+ }
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+ template<int dim>
+ double
+ Q1WedgeFunction<dim>::value (const Point<dim> &p,
const unsigned int) const
-{
- return p*2;
-}
-
-
-
-template<int dim>
-void
-SquareFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
+ {
+ return p(0)*p(1);
+ }
+
+
+
+ template<int dim>
+ void
+ Q1WedgeFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0; i<points.size(); ++i)
- gradients[i] = points[i]*2;
-}
-
-
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ values[i] = p(0)*p(1);
+ }
+ }
+
+
+ template<int dim>
+ double
+ Q1WedgeFunction<dim>::laplacian (const Point<dim> &,
+ const unsigned int) const
+ {
+ return 0.;
+ }
+
+
+ template<int dim>
+ void
+ Q1WedgeFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ values[i] = 0.;
+ }
+
+
+
+ template<int dim>
+ Tensor<1,dim>
+ Q1WedgeFunction<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> erg;
+ erg[0] = p(1);
+ erg[1] = p(0);
+ return erg;
+ }
+
+
+
+ template<int dim>
+ void
+ Q1WedgeFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ {
+ gradients[i][0] = points[i](1);
+ gradients[i][1] = points[i](0);
+ }
+ }
+
+
//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-double
-Q1WedgeFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- return p(0)*p(1);
-}
-
-
-
-template<int dim>
-void
-Q1WedgeFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+
+
+ template<int dim>
+ PillowFunction<dim>::PillowFunction (const double offset)
+ :
+ offset(offset)
+ {}
+
+
+ template<int dim>
+ double
+ PillowFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ switch(dim)
+ {
+ case 1:
+ return 1.-p(0)*p(0)+offset;
+ case 2:
+ return (1.-p(0)*p(0))*(1.-p(1)*p(1))+offset;
+ case 3:
+ return (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2))+offset;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ PillowFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ values[i] = 1.-p(0)*p(0)+offset;
+ break;
+ case 2:
+ values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1))+offset;
+ break;
+ case 3:
+ values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2))+offset;
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+
+
+ template<int dim>
+ double
+ PillowFunction<dim>::laplacian (const Point<dim> &p,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- values[i] = p(0)*p(1);
- }
-}
-
-
-template<int dim>
-double
-Q1WedgeFunction<dim>::laplacian (const Point<dim> &,
+ {
+ switch(dim)
+ {
+ case 1:
+ return -2.;
+ case 2:
+ return -2.*((1.-p(0)*p(0))+(1.-p(1)*p(1)));
+ case 3:
+ return -2.*((1.-p(0)*p(0))*(1.-p(1)*p(1))
+ +(1.-p(1)*p(1))*(1.-p(2)*p(2))
+ +(1.-p(2)*p(2))*(1.-p(0)*p(0)));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ PillowFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ values[i] = -2.;
+ break;
+ case 2:
+ values[i] = -2.*((1.-p(0)*p(0))+(1.-p(1)*p(1)));
+ break;
+ case 3:
+ values[i] = -2.*((1.-p(0)*p(0))*(1.-p(1)*p(1))
+ +(1.-p(1)*p(1))*(1.-p(2)*p(2))
+ +(1.-p(2)*p(2))*(1.-p(0)*p(0)));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template<int dim>
+ Tensor<1,dim>
+ PillowFunction<dim>::gradient (const Point<dim> &p,
const unsigned int) const
-{
- return 0.;
-}
-
-
-template<int dim>
-void
-Q1WedgeFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+ {
+ Tensor<1,dim> result;
+ switch(dim)
+ {
+ case 1:
+ result[0] = -2.*p(0);
+ break;
+ case 2:
+ result[0] = -2.*p(0)*(1.-p(1)*p(1));
+ result[1] = -2.*p(1)*(1.-p(0)*p(0));
+ break;
+ case 3:
+ result[0] = -2.*p(0)*(1.-p(1)*p(1))*(1.-p(2)*p(2));
+ result[1] = -2.*p(1)*(1.-p(0)*p(0))*(1.-p(2)*p(2));
+ result[2] = -2.*p(2)*(1.-p(0)*p(0))*(1.-p(1)*p(1));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return result;
+ }
+
+ template<int dim>
+ void
+ PillowFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- values[i] = 0.;
-}
-
-
-
-template<int dim>
-Tensor<1,dim>
-Q1WedgeFunction<dim>::gradient (const Point<dim> &p,
- const unsigned int) const
-{
- Tensor<1,dim> erg;
- erg[0] = p(1);
- erg[1] = p(0);
- return erg;
-}
-
-
-
-template<int dim>
-void
-Q1WedgeFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0; i<points.size(); ++i)
- {
- gradients[i][0] = points[i](1);
- gradients[i][1] = points[i](0);
- }
-}
-
-
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ gradients[i][0] = -2.*p(0);
+ break;
+ case 2:
+ gradients[i][0] = -2.*p(0)*(1.-p(1)*p(1));
+ gradients[i][1] = -2.*p(1)*(1.-p(0)*p(0));
+ break;
+ case 3:
+ gradients[i][0] = -2.*p(0)*(1.-p(1)*p(1))*(1.-p(2)*p(2));
+ gradients[i][1] = -2.*p(1)*(1.-p(0)*p(0))*(1.-p(2)*p(2));
+ gradients[i][2] = -2.*p(2)*(1.-p(0)*p(0))*(1.-p(1)*p(1));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
//////////////////////////////////////////////////////////////////////
-
-
-template<int dim>
-PillowFunction<dim>::PillowFunction (const double offset)
- :
- offset(offset)
-{}
-
-
-template<int dim>
-double
-PillowFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return 1.-p(0)*p(0)+offset;
- case 2:
- return (1.-p(0)*p(0))*(1.-p(1)*p(1))+offset;
- case 3:
- return (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2))+offset;
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-PillowFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+
+ template<int dim>
+ double
+ CosineFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ switch(dim)
+ {
+ case 1:
+ return std::cos(M_PI_2*p(0));
+ case 2:
+ return std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ case 3:
+ return std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ CosineFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ values[i] = std::cos(M_PI_2*p(0));
+ break;
+ case 2:
+ values[i] = std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ break;
+ case 3:
+ values[i] = std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template<int dim>
+ double
+ CosineFunction<dim>::laplacian (const Point<dim> &p,
+ const unsigned int) const
+ {
+ switch(dim)
+ {
+ case 1:
+ return -M_PI_2*M_PI_2* std::cos(M_PI_2*p(0));
+ case 2:
+ return -2*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ case 3:
+ return -3*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ CosineFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ values[i] = -M_PI_2*M_PI_2* std::cos(M_PI_2*p(0));
+ break;
+ case 2:
+ values[i] = -2*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ break;
+ case 3:
+ values[i] = -3*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template<int dim>
+ Tensor<1,dim>
+ CosineFunction<dim>::gradient (const Point<dim> &p,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = 1.-p(0)*p(0)+offset;
- break;
- case 2:
- values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1))+offset;
- break;
- case 3:
- values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2))+offset;
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-
-
-template<int dim>
-double
-PillowFunction<dim>::laplacian (const Point<dim> &p,
+ {
+ Tensor<1,dim> result;
+ switch(dim)
+ {
+ case 1:
+ result[0] = -M_PI_2* std::sin(M_PI_2*p(0));
+ break;
+ case 2:
+ result[0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ result[1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
+ break;
+ case 3:
+ result[0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ result[1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ result[2] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return result;
+ }
+
+ template<int dim>
+ void
+ CosineFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0));
+ break;
+ case 2:
+ gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ gradients[i][1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
+ break;
+ case 3:
+ gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ gradients[i][1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ gradients[i][2] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template<int dim>
+ Tensor<2,dim>
+ CosineFunction<dim>::hessian (const Point<dim> &p,
const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return -2.;
- case 2:
- return -2.*((1.-p(0)*p(0))+(1.-p(1)*p(1)));
- case 3:
- return -2.*((1.-p(0)*p(0))*(1.-p(1)*p(1))
- +(1.-p(1)*p(1))*(1.-p(2)*p(2))
- +(1.-p(2)*p(2))*(1.-p(0)*p(0)));
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-PillowFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+ {
+ const double pi2 = M_PI_2*M_PI_2;
+
+ Tensor<2,dim> result;
+ switch(dim)
+ {
+ case 1:
+ result[0][0] = -pi2* std::cos(M_PI_2*p(0));
+ break;
+ case 2:
+ if (true)
+ {
+ const double coco = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
+ result[0][0] = coco;
+ result[1][1] = coco;
+ result[0][1] = sisi;
+ result[1][0] = sisi;
+ }
+ break;
+ case 3:
+ if (true)
+ {
+ const double cococo = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ const double sisico = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ const double sicosi = pi2*std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+ const double cosisi = pi2*std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+
+ result[0][0] = cococo;
+ result[1][1] = cococo;
+ result[2][2] = cococo;
+ result[0][1] = sisico;
+ result[1][0] = sisico;
+ result[0][2] = sicosi;
+ result[2][0] = sicosi;
+ result[1][2] = cosisi;
+ result[2][1] = cosisi;
+ }
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return result;
+ }
+
+ template<int dim>
+ void
+ CosineFunction<dim>::hessian_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<2,dim> > &hessians,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = -2.;
- break;
- case 2:
- values[i] = -2.*((1.-p(0)*p(0))+(1.-p(1)*p(1)));
- break;
- case 3:
- values[i] = -2.*((1.-p(0)*p(0))*(1.-p(1)*p(1))
- +(1.-p(1)*p(1))*(1.-p(2)*p(2))
- +(1.-p(2)*p(2))*(1.-p(0)*p(0)));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-Tensor<1,dim>
-PillowFunction<dim>::gradient (const Point<dim> &p,
- const unsigned int) const
-{
- Tensor<1,dim> result;
- switch(dim)
- {
- case 1:
- result[0] = -2.*p(0);
- break;
- case 2:
- result[0] = -2.*p(0)*(1.-p(1)*p(1));
- result[1] = -2.*p(1)*(1.-p(0)*p(0));
- break;
- case 3:
- result[0] = -2.*p(0)*(1.-p(1)*p(1))*(1.-p(2)*p(2));
- result[1] = -2.*p(1)*(1.-p(0)*p(0))*(1.-p(2)*p(2));
- result[2] = -2.*p(2)*(1.-p(0)*p(0))*(1.-p(1)*p(1));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- return result;
-}
-
-template<int dim>
-void
-PillowFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- gradients[i][0] = -2.*p(0);
- break;
- case 2:
- gradients[i][0] = -2.*p(0)*(1.-p(1)*p(1));
- gradients[i][1] = -2.*p(1)*(1.-p(0)*p(0));
- break;
- case 3:
- gradients[i][0] = -2.*p(0)*(1.-p(1)*p(1))*(1.-p(2)*p(2));
- gradients[i][1] = -2.*p(1)*(1.-p(0)*p(0))*(1.-p(2)*p(2));
- gradients[i][2] = -2.*p(2)*(1.-p(0)*p(0))*(1.-p(1)*p(1));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
+ {
+ Assert (hessians.size() == points.size(),
+ ExcDimensionMismatch(hessians.size(), points.size()));
+
+ const double pi2 = M_PI_2*M_PI_2;
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ hessians[i][0][0] = -pi2* std::cos(M_PI_2*p(0));
+ break;
+ case 2:
+ if (true)
+ {
+ const double coco = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
+ const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
+ hessians[i][0][0] = coco;
+ hessians[i][1][1] = coco;
+ hessians[i][0][1] = sisi;
+ hessians[i][1][0] = sisi;
+ }
+ break;
+ case 3:
+ if (true)
+ {
+ const double cococo = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ const double sisico = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
+ const double sicosi = pi2*std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+ const double cosisi = pi2*std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
+
+ hessians[i][0][0] = cococo;
+ hessians[i][1][1] = cococo;
+ hessians[i][2][2] = cococo;
+ hessians[i][0][1] = sisico;
+ hessians[i][1][0] = sisico;
+ hessians[i][0][2] = sicosi;
+ hessians[i][2][0] = sicosi;
+ hessians[i][1][2] = cosisi;
+ hessians[i][2][1] = cosisi;
+ }
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
//////////////////////////////////////////////////////////////////////
-
-template<int dim>
-double
-CosineFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return std::cos(M_PI_2*p(0));
- case 2:
- return std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- case 3:
- return std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-CosineFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = std::cos(M_PI_2*p(0));
- break;
- case 2:
- values[i] = std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- break;
- case 3:
- values[i] = std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-double
-CosineFunction<dim>::laplacian (const Point<dim> &p,
+
+ template<int dim>
+ double
+ ExpFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ switch(dim)
+ {
+ case 1:
+ return std::exp(p(0));
+ case 2:
+ return std::exp(p(0)) * std::exp(p(1));
+ case 3:
+ return std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ ExpFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return -M_PI_2*M_PI_2* std::cos(M_PI_2*p(0));
- case 2:
- return -2*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- case 3:
- return -3*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-CosineFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = -M_PI_2*M_PI_2* std::cos(M_PI_2*p(0));
- break;
- case 2:
- values[i] = -2*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- break;
- case 3:
- values[i] = -3*M_PI_2*M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-Tensor<1,dim>
-CosineFunction<dim>::gradient (const Point<dim> &p,
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
+ {
+ case 1:
+ values[i] = std::exp(p(0));
+ break;
+ case 2:
+ values[i] = std::exp(p(0)) * std::exp(p(1));
+ break;
+ case 3:
+ values[i] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ }
+
+ template<int dim>
+ double
+ ExpFunction<dim>::laplacian (const Point<dim> &p,
const unsigned int) const
-{
- Tensor<1,dim> result;
- switch(dim)
- {
- case 1:
- result[0] = -M_PI_2* std::sin(M_PI_2*p(0));
- break;
- case 2:
- result[0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- result[1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
- break;
- case 3:
- result[0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- result[1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- result[2] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- return result;
-}
-
-template<int dim>
-void
-CosineFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
+ {
+ switch(dim)
+ {
+ case 1:
+ return std::exp(p(0));
+ case 2:
+ return 2 * std::exp(p(0)) * std::exp(p(1));
+ case 3:
+ return 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+ }
+
+ template<int dim>
+ void
+ ExpFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
+ std::vector<double> &values,
const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0));
- break;
- case 2:
- gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- gradients[i][1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
- break;
- case 3:
- gradients[i][0] = -M_PI_2* std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- gradients[i][1] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- gradients[i][2] = -M_PI_2* std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-Tensor<2,dim>
-CosineFunction<dim>::hessian (const Point<dim> &p,
- const unsigned int) const
-{
- const double pi2 = M_PI_2*M_PI_2;
-
- Tensor<2,dim> result;
- switch(dim)
- {
- case 1:
- result[0][0] = -pi2* std::cos(M_PI_2*p(0));
- break;
- case 2:
- if (true)
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
{
- const double coco = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
- result[0][0] = coco;
- result[1][1] = coco;
- result[0][1] = sisi;
- result[1][0] = sisi;
+ case 1:
+ values[i] = std::exp(p(0));
+ break;
+ case 2:
+ values[i] = 2 * std::exp(p(0)) * std::exp(p(1));
+ break;
+ case 3:
+ values[i] = 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
}
- break;
- case 3:
- if (true)
+ }
+ }
+
+ template<int dim>
+ Tensor<1,dim>
+ ExpFunction<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> result;
+ switch(dim)
+ {
+ case 1:
+ result[0] = std::exp(p(0));
+ break;
+ case 2:
+ result[0] = std::exp(p(0)) * std::exp(p(1));
+ result[1] = std::exp(p(0)) * std::exp(p(1));
+ break;
+ case 3:
+ result[0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ result[1] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ result[2] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return result;
+ }
+
+ template<int dim>
+ void
+ ExpFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<dim>& p = points[i];
+ switch(dim)
{
- const double cococo = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- const double sisico = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- const double sicosi = pi2*std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
- const double cosisi = pi2*std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
-
- result[0][0] = cococo;
- result[1][1] = cococo;
- result[2][2] = cococo;
- result[0][1] = sisico;
- result[1][0] = sisico;
- result[0][2] = sicosi;
- result[2][0] = sicosi;
- result[1][2] = cosisi;
- result[2][1] = cosisi;
+ case 1:
+ gradients[i][0] = std::exp(p(0));
+ break;
+ case 2:
+ gradients[i][0] = std::exp(p(0)) * std::exp(p(1));
+ gradients[i][1] = std::exp(p(0)) * std::exp(p(1));
+ break;
+ case 3:
+ gradients[i][0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ gradients[i][1] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ gradients[i][2] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
}
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- return result;
-}
-
-template<int dim>
-void
-CosineFunction<dim>::hessian_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<2,dim> > &hessians,
- const unsigned int) const
-{
- Assert (hessians.size() == points.size(),
- ExcDimensionMismatch(hessians.size(), points.size()));
-
- const double pi2 = M_PI_2*M_PI_2;
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- hessians[i][0][0] = -pi2* std::cos(M_PI_2*p(0));
- break;
- case 2:
- if (true)
- {
- const double coco = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1));
- const double sisi = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1));
- hessians[i][0][0] = coco;
- hessians[i][1][1] = coco;
- hessians[i][0][1] = sisi;
- hessians[i][1][0] = sisi;
- }
- break;
- case 3:
- if (true)
- {
- const double cococo = -pi2*std::cos(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- const double sisico = pi2*std::sin(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::cos(M_PI_2*p(2));
- const double sicosi = pi2*std::sin(M_PI_2*p(0)) * std::cos(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
- const double cosisi = pi2*std::cos(M_PI_2*p(0)) * std::sin(M_PI_2*p(1)) * std::sin(M_PI_2*p(2));
-
- hessians[i][0][0] = cococo;
- hessians[i][1][1] = cococo;
- hessians[i][2][2] = cococo;
- hessians[i][0][1] = sisico;
- hessians[i][1][0] = sisico;
- hessians[i][0][2] = sicosi;
- hessians[i][2][0] = sicosi;
- hessians[i][1][2] = cosisi;
- hessians[i][2][1] = cosisi;
- }
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
+ }
+ }
+
//////////////////////////////////////////////////////////////////////
-
-template<int dim>
-double
-ExpFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return std::exp(p(0));
- case 2:
- return std::exp(p(0)) * std::exp(p(1));
- case 3:
- return std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-ExpFunction<dim>::value_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = std::exp(p(0));
- break;
- case 2:
- values[i] = std::exp(p(0)) * std::exp(p(1));
- break;
- case 3:
- values[i] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-double
-ExpFunction<dim>::laplacian (const Point<dim> &p,
- const unsigned int) const
-{
- switch(dim)
- {
- case 1:
- return std::exp(p(0));
- case 2:
- return 2 * std::exp(p(0)) * std::exp(p(1));
- case 3:
- return 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- default:
- Assert(false, ExcNotImplemented());
- }
- return 0.;
-}
-
-template<int dim>
-void
-ExpFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &points,
- std::vector<double> &values,
+
+
+ double
+ LSingularityFunction::value (const Point<2> &p,
+ const unsigned int) const
+ {
+ double x = p(0);
+ double y = p(1);
+
+ if ((x>=0) && (y>=0))
+ return 0.;
+
+ double phi = std::atan2(y,-x)+M_PI;
+ double r2 = x*x+y*y;
+
+ return std::pow(r2,1./3.) * std::sin(2./3.*phi);
+ }
+
+
+ void
+ LSingularityFunction::value_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ double x = points[i](0);
+ double y = points[i](1);
+
+ if ((x>=0) && (y>=0))
+ values[i] = 0.;
+ else
+ {
+ double phi = std::atan2(y,-x)+M_PI;
+ double r2 = x*x+y*y;
+
+ values[i] = std::pow(r2,1./3.) * std::sin(2./3.*phi);
+ }
+ }
+ }
+
+
+ double
+ LSingularityFunction::laplacian (const Point<2> &,
+ const unsigned int) const
+ {
+ return 0.;
+ }
+
+
+ void
+ LSingularityFunction::laplacian_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ values[i] = 0.;
+ }
+
+
+ Tensor<1,2>
+ LSingularityFunction::gradient (const Point<2> &p,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- values[i] = std::exp(p(0));
- break;
- case 2:
- values[i] = 2 * std::exp(p(0)) * std::exp(p(1));
- break;
- case 3:
- values[i] = 3 * std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
-template<int dim>
-Tensor<1,dim>
-ExpFunction<dim>::gradient (const Point<dim> &p,
- const unsigned int) const
-{
- Tensor<1,dim> result;
- switch(dim)
- {
- case 1:
- result[0] = std::exp(p(0));
- break;
- case 2:
- result[0] = std::exp(p(0)) * std::exp(p(1));
- result[1] = std::exp(p(0)) * std::exp(p(1));
- break;
- case 3:
- result[0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- result[1] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- result[2] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- return result;
-}
-
-template<int dim>
-void
-ExpFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &points,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<dim>& p = points[i];
- switch(dim)
- {
- case 1:
- gradients[i][0] = std::exp(p(0));
- break;
- case 2:
- gradients[i][0] = std::exp(p(0)) * std::exp(p(1));
- gradients[i][1] = std::exp(p(0)) * std::exp(p(1));
- break;
- case 3:
- gradients[i][0] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- gradients[i][1] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- gradients[i][2] = std::exp(p(0)) * std::exp(p(1)) * std::exp(p(2));
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- }
-}
-
+ {
+ double x = p(0);
+ double y = p(1);
+ double phi = std::atan2(y,-x)+M_PI;
+ double r43 = std::pow(x*x+y*y,2./3.);
+
+ Tensor<1,2> result;
+ result[0] = 2./3.*(std::sin(2./3.*phi)*x + std::cos(2./3.*phi)*y)/r43;
+ result[1] = 2./3.*(std::sin(2./3.*phi)*y - std::cos(2./3.*phi)*x)/r43;
+ return result;
+ }
+
+
+ void
+ LSingularityFunction::gradient_list (const std::vector<Point<2> > &points,
+ std::vector<Tensor<1,2> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<2>& p = points[i];
+ double x = p(0);
+ double y = p(1);
+ double phi = std::atan2(y,-x)+M_PI;
+ double r43 = std::pow(x*x+y*y,2./3.);
+
+ gradients[i][0] = 2./3.*(std::sin(2./3.*phi)*x + std::cos(2./3.*phi)*y)/r43;
+ gradients[i][1] = 2./3.*(std::sin(2./3.*phi)*y - std::cos(2./3.*phi)*x)/r43;
+ }
+ }
+
//////////////////////////////////////////////////////////////////////
-
-
-double
-LSingularityFunction::value (const Point<2> &p,
- const unsigned int) const
-{
- double x = p(0);
- double y = p(1);
-
- if ((x>=0) && (y>=0))
- return 0.;
- double phi = std::atan2(y,-x)+M_PI;
- double r2 = x*x+y*y;
-
- return std::pow(r2,1./3.) * std::sin(2./3.*phi);
-}
-
-
-void
-LSingularityFunction::value_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
+
+ double
+ SlitSingularityFunction::value (const Point<2> &p,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- double x = points[i](0);
- double y = points[i](1);
-
- if ((x>=0) && (y>=0))
- values[i] = 0.;
- else
- {
- double phi = std::atan2(y,-x)+M_PI;
- double r2 = x*x+y*y;
-
- values[i] = std::pow(r2,1./3.) * std::sin(2./3.*phi);
- }
- }
-}
-
-
-double
-LSingularityFunction::laplacian (const Point<2> &,
- const unsigned int) const
-{
- return 0.;
-}
-
-
-void
-LSingularityFunction::laplacian_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
+ {
+ double x = p(0);
+ double y = p(1);
+
+ double phi = std::atan2(x,y)+M_PI;
+ double r2 = x*x+y*y;
+
+ return std::pow(r2,.25) * std::sin(.5*phi);
+ }
+
+
+ void
+ SlitSingularityFunction::value_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ double x = points[i](0);
+ double y = points[i](1);
+
+ double phi = std::atan2(x,y)+M_PI;
+ double r2 = x*x+y*y;
+
+ values[i] = std::pow(r2,.25) * std::sin(.5*phi);
+ }
+ }
+
+
+ double
+ SlitSingularityFunction::laplacian (const Point<2> &,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- values[i] = 0.;
-}
-
-
-Tensor<1,2>
-LSingularityFunction::gradient (const Point<2> &p,
- const unsigned int) const
-{
- double x = p(0);
- double y = p(1);
- double phi = std::atan2(y,-x)+M_PI;
- double r43 = std::pow(x*x+y*y,2./3.);
-
- Tensor<1,2> result;
- result[0] = 2./3.*(std::sin(2./3.*phi)*x + std::cos(2./3.*phi)*y)/r43;
- result[1] = 2./3.*(std::sin(2./3.*phi)*y - std::cos(2./3.*phi)*x)/r43;
- return result;
-}
-
-
-void
-LSingularityFunction::gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
+ {
+ return 0.;
+ }
+
+
+ void
+ SlitSingularityFunction::laplacian_list (const std::vector<Point<2> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch(values.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ values[i] = 0.;
+ }
+
+
+ Tensor<1,2>
+ SlitSingularityFunction::gradient (const Point<2> &p,
const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<2>& p = points[i];
- double x = p(0);
- double y = p(1);
- double phi = std::atan2(y,-x)+M_PI;
- double r43 = std::pow(x*x+y*y,2./3.);
-
- gradients[i][0] = 2./3.*(std::sin(2./3.*phi)*x + std::cos(2./3.*phi)*y)/r43;
- gradients[i][1] = 2./3.*(std::sin(2./3.*phi)*y - std::cos(2./3.*phi)*x)/r43;
- }
-}
-
+ {
+ double x = p(0);
+ double y = p(1);
+ double phi = std::atan2(x,y)+M_PI;
+ double r64 = std::pow(x*x+y*y,3./4.);
+
+
+ Tensor<1,2> result;
+ result[0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64;
+ result[1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64;
+ return result;
+ }
+
+
+ void
+ SlitSingularityFunction::gradient_list (const std::vector<Point<2> > &points,
+ std::vector<Tensor<1,2> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == points.size(),
+ ExcDimensionMismatch(gradients.size(), points.size()));
+
+ for (unsigned int i=0;i<points.size();++i)
+ {
+ const Point<2>& p = points[i];
+ double x = p(0);
+ double y = p(1);
+ double phi = std::atan2(x,y)+M_PI;
+ double r64 = std::pow(x*x+y*y,3./4.);
+
+ gradients[i][0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64;
+ gradients[i][1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64;
+ }
+ }
+
//////////////////////////////////////////////////////////////////////
-
-
-double
-SlitSingularityFunction::value (const Point<2> &p,
+
+ template<int dim>
+ JumpFunction<dim>::JumpFunction(const Point<dim> &direction,
+ const double steepness)
+ :
+ direction(direction),
+ steepness(steepness)
+ {
+ switch (dim)
+ {
+ case 1:
+ angle = 0;
+ break;
+ case 2:
+ angle = std::atan2(direction(0),direction(1));
+ break;
+ case 3:
+ Assert(false, ExcNotImplemented());
+ }
+ sine = std::sin(angle);
+ cosine = std::cos(angle);
+ }
+
+
+
+ template<int dim>
+ double
+ JumpFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ double x = steepness*(-cosine*p(0)+sine*p(1));
+ return -std::atan(x);
+ }
+
+
+
+ template<int dim>
+ void
+ JumpFunction<dim>::value_list (const typename std::vector<Point<dim> > &p,
+ std::vector<double> &values,
+ const unsigned int) const
+ {
+ Assert (values.size() == p.size(),
+ ExcDimensionMismatch(values.size(), p.size()));
+
+ for (unsigned int i=0;i<p.size();++i)
+ {
+ double x = steepness*(-cosine*p[i](0)+sine*p[i](1));
+ values[i] = -std::atan(x);
+ }
+ }
+
+
+ template<int dim>
+ double
+ JumpFunction<dim>::laplacian (const Point<dim> &p,
const unsigned int) const
-{
- double x = p(0);
- double y = p(1);
-
- double phi = std::atan2(x,y)+M_PI;
- double r2 = x*x+y*y;
-
- return std::pow(r2,.25) * std::sin(.5*phi);
-}
-
-
-void
-SlitSingularityFunction::value_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
+ {
+ double x = steepness*(-cosine*p(0)+sine*p(1));
+ double r = 1+x*x;
+ return 2*steepness*steepness*x/(r*r);
+ }
+
+
+ template<int dim>
+ void
+ JumpFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &p,
+ std::vector<double> &values,
const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- double x = points[i](0);
- double y = points[i](1);
-
- double phi = std::atan2(x,y)+M_PI;
- double r2 = x*x+y*y;
-
- values[i] = std::pow(r2,.25) * std::sin(.5*phi);
- }
-}
-
-
-double
-SlitSingularityFunction::laplacian (const Point<2> &,
- const unsigned int) const
-{
- return 0.;
-}
-
-
-void
-SlitSingularityFunction::laplacian_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == points.size(),
- ExcDimensionMismatch(values.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- values[i] = 0.;
-}
-
-
-Tensor<1,2>
-SlitSingularityFunction::gradient (const Point<2> &p,
- const unsigned int) const
-{
- double x = p(0);
- double y = p(1);
- double phi = std::atan2(x,y)+M_PI;
- double r64 = std::pow(x*x+y*y,3./4.);
+ {
+ Assert (values.size() == p.size(),
+ ExcDimensionMismatch(values.size(), p.size()));
+
+ double f = 2*steepness*steepness;
+
+ for (unsigned int i=0;i<p.size();++i)
+ {
+ double x = steepness*(-cosine*p[i](0)+sine*p[i](1));
+ double r = 1+x*x;
+ values[i] = f*x/(r*r);
+ }
+ }
- Tensor<1,2> result;
- result[0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64;
- result[1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64;
- return result;
-}
-
-
-void
-SlitSingularityFunction::gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
- const unsigned int) const
-{
- Assert (gradients.size() == points.size(),
- ExcDimensionMismatch(gradients.size(), points.size()));
-
- for (unsigned int i=0;i<points.size();++i)
- {
- const Point<2>& p = points[i];
- double x = p(0);
- double y = p(1);
- double phi = std::atan2(x,y)+M_PI;
- double r64 = std::pow(x*x+y*y,3./4.);
-
- gradients[i][0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64;
- gradients[i][1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64;
- }
-}
-
-//////////////////////////////////////////////////////////////////////
-
-template<int dim>
-JumpFunction<dim>::JumpFunction(const Point<dim> &direction,
- const double steepness)
- :
- direction(direction),
- steepness(steepness)
-{
- switch (dim)
- {
- case 1:
- angle = 0;
- break;
- case 2:
- angle = std::atan2(direction(0),direction(1));
- break;
- case 3:
- Assert(false, ExcNotImplemented());
- }
- sine = std::sin(angle);
- cosine = std::cos(angle);
-}
-
-
-
-template<int dim>
-double
-JumpFunction<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- double x = steepness*(-cosine*p(0)+sine*p(1));
- return -std::atan(x);
-}
-
-
-
-template<int dim>
-void
-JumpFunction<dim>::value_list (const typename std::vector<Point<dim> > &p,
- std::vector<double> &values,
+
+ template<int dim>
+ Tensor<1,dim>
+ JumpFunction<dim>::gradient (const Point<dim> &p,
const unsigned int) const
-{
- Assert (values.size() == p.size(),
- ExcDimensionMismatch(values.size(), p.size()));
-
- for (unsigned int i=0;i<p.size();++i)
- {
- double x = steepness*(-cosine*p[i](0)+sine*p[i](1));
- values[i] = -std::atan(x);
- }
-}
-
-
-template<int dim>
-double
-JumpFunction<dim>::laplacian (const Point<dim> &p,
- const unsigned int) const
-{
- double x = steepness*(-cosine*p(0)+sine*p(1));
- double r = 1+x*x;
- return 2*steepness*steepness*x/(r*r);
-}
-
-
-template<int dim>
-void
-JumpFunction<dim>::laplacian_list (const typename std::vector<Point<dim> > &p,
- std::vector<double> &values,
- const unsigned int) const
-{
- Assert (values.size() == p.size(),
- ExcDimensionMismatch(values.size(), p.size()));
-
- double f = 2*steepness*steepness;
-
- for (unsigned int i=0;i<p.size();++i)
- {
- double x = steepness*(-cosine*p[i](0)+sine*p[i](1));
- double r = 1+x*x;
- values[i] = f*x/(r*r);
- }
-}
-
-
-
-template<int dim>
-Tensor<1,dim>
-JumpFunction<dim>::gradient (const Point<dim> &p,
- const unsigned int) const
-{
- double x = steepness*(-cosine*p(0)+sine*p(1));
- double r = -steepness*(1+x*x);
- Tensor<1,dim> erg;
- erg[0] = cosine*r;
- erg[1] = sine*r;
- return erg;
-}
-
-
-
-template<int dim>
-void
-JumpFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &p,
- typename std::vector<Tensor<1,dim> > &gradients,
- const unsigned int) const
-{
- Assert (gradients.size() == p.size(),
- ExcDimensionMismatch(gradients.size(), p.size()));
-
- for (unsigned int i=0; i<p.size(); ++i)
- {
- double x = steepness*(cosine*p[i](0)+sine*p[i](1));
- double r = -steepness*(1+x*x);
- gradients[i][0] = cosine*r;
- gradients[i][1] = sine*r;
- }
-}
-
-
-
-template <int dim>
-unsigned int
-JumpFunction<dim>::memory_consumption () const
-{
- // only simple data elements, so
- // use sizeof operator
- return sizeof (*this);
-};
-
-
-
-
-
+ {
+ double x = steepness*(-cosine*p(0)+sine*p(1));
+ double r = -steepness*(1+x*x);
+ Tensor<1,dim> erg;
+ erg[0] = cosine*r;
+ erg[1] = sine*r;
+ return erg;
+ }
+
+
+
+ template<int dim>
+ void
+ JumpFunction<dim>::gradient_list (const typename std::vector<Point<dim> > &p,
+ typename std::vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+ {
+ Assert (gradients.size() == p.size(),
+ ExcDimensionMismatch(gradients.size(), p.size()));
+
+ for (unsigned int i=0; i<p.size(); ++i)
+ {
+ double x = steepness*(cosine*p[i](0)+sine*p[i](1));
+ double r = -steepness*(1+x*x);
+ gradients[i][0] = cosine*r;
+ gradients[i][1] = sine*r;
+ }
+ }
+
+
+
+ template <int dim>
+ unsigned int
+ JumpFunction<dim>::memory_consumption () const
+ {
+ // only simple data elements, so
+ // use sizeof operator
+ return sizeof (*this);
+ };
+
+
+
+
+
/* ---------------------- FourierSineFunction ----------------------- */
-
-
-template <int dim>
-FourierCosineFunction<dim>::
-FourierCosineFunction (const Point<dim> &fourier_coefficients)
- :
- Function<dim> (1),
- fourier_coefficients (fourier_coefficients)
-{};
-
-
-
-template <int dim>
-double
-FourierCosineFunction<dim>::value (const Point<dim> &p,
+
+
+ template <int dim>
+ FourierCosineFunction<dim>::
+ FourierCosineFunction (const Point<dim> &fourier_coefficients)
+ :
+ Function<dim> (1),
+ fourier_coefficients (fourier_coefficients)
+ {};
+
+
+
+ template <int dim>
+ double
+ FourierCosineFunction<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ double val=1;
+ for (unsigned int i=0; i<dim; ++i)
+ val *= std::cos(fourier_coefficients[i] * p[i]);
+ return val;
+ };
+
+
+
+ template <int dim>
+ Tensor<1,dim>
+ FourierCosineFunction<dim>::gradient (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ Tensor<1,dim> grad;
+ for (unsigned int i=0; i<dim; ++i)
+ grad[i] = 1;
+
+ for (unsigned int i=0; i<dim; ++i)
+ {
+ const double cos_i = std::cos(fourier_coefficients[i] * p[i]);
+ const double sin_i = std::sin(fourier_coefficients[i] * p[i]);
+
+ for (unsigned int d=0; d<dim; ++d)
+ if (d==i)
+ grad[d] *= - fourier_coefficients[i] * sin_i;
+ else
+ grad[d] *= cos_i;
+ };
+
+ return grad;
+ };
+
+
+
+ template <int dim>
+ double
+ FourierCosineFunction<dim>::laplacian (const Point<dim> &p,
+ const unsigned int component) const
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ double val = -(fourier_coefficients*fourier_coefficients);
+ for (unsigned int i=0; i<dim; ++i)
+ val *= std::cos(fourier_coefficients[i] * p[i]);
+ return val;
+ };
+
+
+
+
+/* ---------------------- FourierSineFunction ----------------------- */
+
+
+
+ template <int dim>
+ FourierSineFunction<dim>::
+ FourierSineFunction (const Point<dim> &fourier_coefficients)
+ :
+ Function<dim> (1),
+ fourier_coefficients (fourier_coefficients)
+ {};
+
+
+
+ template <int dim>
+ double
+ FourierSineFunction<dim>::value (const Point<dim> &p,
const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- double val=1;
- for (unsigned int i=0; i<dim; ++i)
- val *= std::cos(fourier_coefficients[i] * p[i]);
- return val;
-};
-
-
-
-template <int dim>
-Tensor<1,dim>
-FourierCosineFunction<dim>::gradient (const Point<dim> &p,
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ double val=1;
+ for (unsigned int i=0; i<dim; ++i)
+ val *= std::sin(fourier_coefficients[i] * p[i]);
+ return val;
+ };
+
+
+
+ template <int dim>
+ Tensor<1,dim>
+ FourierSineFunction<dim>::gradient (const Point<dim> &p,
const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- Tensor<1,dim> grad;
- for (unsigned int i=0; i<dim; ++i)
- grad[i] = 1;
-
- for (unsigned int i=0; i<dim; ++i)
- {
- const double cos_i = std::cos(fourier_coefficients[i] * p[i]);
- const double sin_i = std::sin(fourier_coefficients[i] * p[i]);
-
- for (unsigned int d=0; d<dim; ++d)
- if (d==i)
- grad[d] *= - fourier_coefficients[i] * sin_i;
- else
- grad[d] *= cos_i;
- };
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ Tensor<1,dim> grad;
+ for (unsigned int i=0; i<dim; ++i)
+ grad[i] = 1;
+
+ for (unsigned int i=0; i<dim; ++i)
+ {
+ const double cos_i = std::cos(fourier_coefficients[i] * p[i]);
+ const double sin_i = std::sin(fourier_coefficients[i] * p[i]);
+
+ for (unsigned int d=0; d<dim; ++d)
+ if (d==i)
+ grad[d] *= fourier_coefficients[i] * cos_i;
+ else
+ grad[d] *= sin_i;
+ };
+
+ return grad;
+ };
- return grad;
-};
-
-
-
-template <int dim>
-double
-FourierCosineFunction<dim>::laplacian (const Point<dim> &p,
+
+
+ template <int dim>
+ double
+ FourierSineFunction<dim>::laplacian (const Point<dim> &p,
const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- double val = -(fourier_coefficients*fourier_coefficients);
- for (unsigned int i=0; i<dim; ++i)
- val *= std::cos(fourier_coefficients[i] * p[i]);
- return val;
-};
-
-
-
-
-/* ---------------------- FourierSineFunction ----------------------- */
-
-
-
-template <int dim>
-FourierSineFunction<dim>::
-FourierSineFunction (const Point<dim> &fourier_coefficients)
- :
- Function<dim> (1),
- fourier_coefficients (fourier_coefficients)
-{};
-
-
-
-template <int dim>
-double
-FourierSineFunction<dim>::value (const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- double val=1;
- for (unsigned int i=0; i<dim; ++i)
- val *= std::sin(fourier_coefficients[i] * p[i]);
- return val;
-};
-
-
-
-template <int dim>
-Tensor<1,dim>
-FourierSineFunction<dim>::gradient (const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- Tensor<1,dim> grad;
- for (unsigned int i=0; i<dim; ++i)
- grad[i] = 1;
-
- for (unsigned int i=0; i<dim; ++i)
- {
- const double cos_i = std::cos(fourier_coefficients[i] * p[i]);
- const double sin_i = std::sin(fourier_coefficients[i] * p[i]);
-
- for (unsigned int d=0; d<dim; ++d)
- if (d==i)
- grad[d] *= fourier_coefficients[i] * cos_i;
- else
- grad[d] *= sin_i;
- };
+ {
+ Assert (component==0, ExcIndexRange(component,0,1));
+ double val = -(fourier_coefficients*fourier_coefficients);
+ for (unsigned int i=0; i<dim; ++i)
+ val *= std::sin(fourier_coefficients[i] * p[i]);
+ return val;
+ };
+
+
+
+
+ template class SquareFunction<1>;
+ template class SquareFunction<2>;
+ template class SquareFunction<3>;
+ template class Q1WedgeFunction<1>;
+ template class Q1WedgeFunction<2>;
+ template class Q1WedgeFunction<3>;
+ template class PillowFunction<1>;
+ template class PillowFunction<2>;
+ template class PillowFunction<3>;
+ template class CosineFunction<1>;
+ template class CosineFunction<2>;
+ template class CosineFunction<3>;
+ template class ExpFunction<1>;
+ template class ExpFunction<2>;
+ template class ExpFunction<3>;
+ template class JumpFunction<1>;
+ template class JumpFunction<2>;
+ template class JumpFunction<3>;
+ template class FourierCosineFunction<1>;
+ template class FourierCosineFunction<2>;
+ template class FourierCosineFunction<3>;
+ template class FourierSineFunction<1>;
+ template class FourierSineFunction<2>;
+ template class FourierSineFunction<3>;
+
- return grad;
-};
-
-
-
-template <int dim>
-double
-FourierSineFunction<dim>::laplacian (const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (component==0, ExcIndexRange(component,0,1));
- double val = -(fourier_coefficients*fourier_coefficients);
- for (unsigned int i=0; i<dim; ++i)
- val *= std::sin(fourier_coefficients[i] * p[i]);
- return val;
};
-
-
-
-
-template class SquareFunction<1>;
-template class SquareFunction<2>;
-template class SquareFunction<3>;
-template class Q1WedgeFunction<1>;
-template class Q1WedgeFunction<2>;
-template class Q1WedgeFunction<3>;
-template class PillowFunction<1>;
-template class PillowFunction<2>;
-template class PillowFunction<3>;
-template class CosineFunction<1>;
-template class CosineFunction<2>;
-template class CosineFunction<3>;
-template class ExpFunction<1>;
-template class ExpFunction<2>;
-template class ExpFunction<3>;
-template class JumpFunction<1>;
-template class JumpFunction<2>;
-template class JumpFunction<3>;
-template class FourierCosineFunction<1>;
-template class FourierCosineFunction<2>;
-template class FourierCosineFunction<3>;
-template class FourierSineFunction<1>;
-template class FourierSineFunction<2>;
-template class FourierSineFunction<3>;
-