const unsigned int component = 0) const;
};
+/**
+ * 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 vector<Point<dim> > &points,
+ 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 vector<Point<dim> > &points,
+ vector<Tensor<1,dim> > &gradients,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at a single point.
+ */
+ double laplacian(const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * Laplacian at multiple points.
+ */
+ void laplacian_list(const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int component = 0) const;
+};
+
+
/**
* Singularity on the L-shaped domain in 2D.
* @author Guido Kanschat, 1999
//////////////////////////////////////////////////////////////////////
+template<int dim>
+double
+CosineFunction<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ switch(dim)
+ {
+ case 1:
+ return cos(M_PI_2*p(0));
+ case 2:
+ return cos(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ case 3:
+ return cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+}
+
+template<int dim>
+void
+CosineFunction<dim>::value_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int) const
+{
+ Assert (values.size() == points.size(),
+ ExcVectorHasWrongSize(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] = cos(M_PI_2*p(0));
+ break;
+ case 2:
+ values[i] = cos(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ break;
+ case 3:
+ values[i] = cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * 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* cos(M_PI_2*p(0));
+ case 2:
+ return -2*M_PI_2*M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ case 3:
+ return -3*M_PI_2*M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+}
+
+template<int dim>
+void
+CosineFunction<dim>::laplacian_list (const vector<Point<dim> > &points,
+ vector<double> &values,
+ const unsigned int) const
+{
+ Assert (values.size() == points.size(),
+ ExcVectorHasWrongSize(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* cos(M_PI_2*p(0));
+ break;
+ case 2:
+ values[i] = -2*M_PI_2*M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ break;
+ case 3:
+ values[i] = -3*M_PI_2*M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * 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
+{
+ Tensor<1,dim> result;
+ switch(dim)
+ {
+ case 1:
+ result[0] = -M_PI_2* sin(M_PI_2*p(0));
+ break;
+ case 2:
+ result[0] = -M_PI_2* sin(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ result[1] = -M_PI_2* cos(M_PI_2*p(0)) * sin(M_PI_2*p(1));
+ break;
+ case 3:
+ result[0] = -M_PI_2* sin(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ result[1] = -M_PI_2* cos(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ result[2] = -M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * sin(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return result;
+}
+
+template<int dim>
+void
+CosineFunction<dim>::gradient_list (const vector<Point<dim> > &points,
+ vector<Tensor<1,dim> > &gradients,
+ const unsigned int) const
+{
+ Assert (gradients.size() == points.size(),
+ ExcVectorHasWrongSize(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* sin(M_PI_2*p(0));
+ break;
+ case 2:
+ gradients[i][0] = -M_PI_2* sin(M_PI_2*p(0)) * cos(M_PI_2*p(1));
+ gradients[i][1] = -M_PI_2* cos(M_PI_2*p(0)) * sin(M_PI_2*p(1));
+ return;
+ case 3:
+ gradients[i][0] = -M_PI_2* sin(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ gradients[i][1] = -M_PI_2* cos(M_PI_2*p(0)) * sin(M_PI_2*p(1)) * cos(M_PI_2*p(2));
+ gradients[i][2] = -M_PI_2* cos(M_PI_2*p(0)) * cos(M_PI_2*p(1)) * sin(M_PI_2*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+}
+
+//////////////////////////////////////////////////////////////////////
+
double
LSingularityFunction::value (const Point<2> &p,
template PillowFunction<1>;
template PillowFunction<2>;
template PillowFunction<3>;
+template CosineFunction<1>;
+template CosineFunction<2>;
+template CosineFunction<3>;