Think about the determinant function in the tensor package. Is it
useful, can it be generalized?
-
+> There are only non-zero alternating real valued n-forms, so a determinant
+> exists only for Tensor<dim,dim>
+> Alternating exterior forms of lower index return tensors.
+> Guido
Put Tensor<1,dim>::array_size directly into the array type. At
present, egcs 1.1.2 chokes on that.
Fill in docs for the timer class. Hopefully finally find a way to
let it measure times larger than half an hour.
-
+> Operating system dependence. Grrr. Guido
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.
+> Remove !? Guido
--- /dev/null
+// $Id$
+// Copyright Guido Kanschat, University of Minnesota, 1999
+
+#ifndef __base__function_lib_h
+#define __base__function_lib_h
+
+#include <base/function.h>
+
+/**
+ * n-quadratic pillow on the unit square.
+ *
+ * This is a function for testing the implementation. It has zero
+ * boundary values on the domain $(-1,1)^d$. In the inside, it is the
+ * product of $x_i^2-1$.
+ *
+ * Together with the function, its derivatives and Laplacian are defined.
+ * @author: Guido Kanschat, 1999
+ */
+template<int dim>
+class PillowFunction : 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;
+};
+
+#endif
--- /dev/null
+// $Id$
+// (c) Guido Kanschat, 1999
+
+#include <base/point.h>
+#include <base/function_lib.h>
+
+//TODO: Derivatives in 3d wrong
+
+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);
+ case 2:
+ return (1.-p(0)*p(0))*(1.-p(1)*p(1));
+ case 3:
+ return (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2));
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return 0.;
+}
+
+template<int dim>
+void
+PillowFunction<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] = 1.-p(0)*p(0);
+ break;
+ case 2:
+ values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1));
+ break;
+ case 3:
+ values[i] = (1.-p(0)*p(0))*(1.-p(1)*p(1))*(1.-p(2)*p(2));
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ }
+}
+
+template<int dim>
+double
+PillowFunction<dim>::laplacian (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 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] = 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(2)*p(2)));
+ 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 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] = 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));
+ return;
+ 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 PillowFunction<1>;
+template PillowFunction<2>;
+template PillowFunction<3>;
+