/*---------------------------- tensor.h ---------------------------*/
-/* $Id$ */
+/* $Id$ */
#ifndef __tensor_H
#define __tensor_H
/*---------------------------- tensor.h ---------------------------*/
--- /dev/null
+/*---------------------------- tensorfunction.h ---------------------------*/
+/* $Id$ */
+/* Copyright G. Kanschat, University of Heidelberg, 1999 */
+#ifndef __tensorfunction_H
+#define __tensorfunction_H
+/*---------------------------- tensorfunction.h ---------------------------*/
+
+
+#include <base/exceptions.h>
+#include <vector>
+#include <base/point.h>
+#include <base/functiontime.h>
+#include <base/tensorindex.h>
+
+
+template<int rank_, int dim>
+class Tensor;
+
+/**
+ * This class is a model for a tensor valued function.
+ * It returns the value
+ * at a given point through the #operator ()# member functions,
+ * which are virtual. It also has a function to return a whole list of
+ * values at different points to reduce the overhead of the virtual function
+ * calls; this function is preset to successively call the function returning
+ * one value at a time.
+ *
+ * There are other functions return the gradient of the function at one or
+ * several points. You only have to overload those functions you need; the
+ * functions returning several values at a time will call those returning
+ * only one value, while those ones will throw an exception when called but
+ * not overloaded.
+ *
+ * Usually, efficiency of your program increases if you overload the
+ * complex virtual functions, too.
+ *
+ * @author Guido Kanschat, 1999
+ */
+template <int rank_, int dim>
+class TensorFunction :
+ public FunctionTime
+{
+ public:
+ /**
+ * Constructor. May take an initial vakue
+ * for the time variable, which defaults
+ * to zero.
+ */
+ TensorFunction (const double initial_time = 0.0);
+
+ /**
+ * Virtual destructor; absolutely
+ * necessary in this case.
+ */
+ virtual ~TensorFunction ();
+
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual Tensor<rank_, dim> operator () (const Point<dim> &p) const;
+
+ /**
+ * Set #values# to the point values
+ * of the function at the #points#.
+ * It is assumed that #values#
+ * already has the right size, i.e.
+ * the same size as the #points#
+ * array.
+ */
+ virtual void value_list (const vector<Point<dim> > &points,
+ vector<Tensor<rank_,dim> > &values) const;
+
+ /**
+ * Return one component of the value.
+ */
+ virtual double operator() (TensorIndex<rank_> i, const Point<dim>& p) const;
+
+
+ /**
+ * Return the gradient of the function
+ * at the given point.
+ */
+ virtual Tensor<rank_+1,dim> gradient (const Point<dim> &p) const;
+
+ /**
+ * Set #gradients# to the gradients of
+ * the function at the #points#.
+ * It is assumed that #values#
+ * already has the right size, i.e.
+ * the same size as the #points# array.
+ */
+ virtual void gradient_list (const vector<Point<dim> > &points,
+ vector<Tensor<rank_+1,dim> > &gradients) const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcPureFunctionCalled);
+ /**
+ * Exception
+ */
+ DeclException2 (ExcVectorHasWrongSize,
+ int, int,
+ << "The vector has size " << arg1 << " but should have "
+ << arg2 << " elements.");
+
+};
+
+
+
+
+
+
+/*---------------------------- tensorfunction.h ---------------------------*/
+/* end of #ifndef __tensorfunction_H */
+#endif
+/*---------------------------- tensorfunction.h ---------------------------*/
--- /dev/null
+/*---------------------------- tensorindex.h ---------------------------*/
+/* $Id$ */
+#ifndef __tensorindex_H
+#define __tensorindex_H
+/*---------------------------- tensorindex.h ---------------------------*/
+
+#include <base/exceptions.h>
+
+/**
+ *
+ * Rank-independent access to elements of #Tensor#. A little class
+ * template remembering #rank# integers.
+ *
+ * This template is implemented only for ranks between 1 and 4. If you
+ * really need higher ranks, please modify the source code accordingly
+ * or contact the developer.
+ *
+ * @author Guido Kanschat, 1999
+ *
+ */
+template<int rank>
+class TensorIndex
+{
+private:
+ /**
+ * Field of indices.
+ */
+ unsigned index[rank];
+public:
+ /**
+ * Constructor taking #rank# indices.
+ */
+ TensorIndex(...);
+
+ /**
+ *
+ * Access operator returning index
+ * in #n#th component
+ *
+ */
+ unsigned operator () (unsigned n) const;
+ DeclException1(ExcRank, unsigned,
+ << "Index " << arg1 << " higher than maximum " << rank-1);
+};
+
+
+template<>
+class TensorIndex<4>
+{
+private:
+ /**
+ * Field of indices.
+ */
+ unsigned index[4];
+public:
+ /**
+ * Constructor taking #rank# indices.
+ */
+ TensorIndex(unsigned i0, unsigned i1, unsigned i2, unsigned i3)
+ {
+ index[0] = i0;
+ index[1] = i1;
+ index[2] = i2;
+ index[3] = i3;
+ }
+
+
+ /**
+ *
+ * Access operator returning index
+ * in #n#th component
+ *
+ */
+ unsigned int operator () (unsigned int n) const
+ {
+ Assert(n<4, ExcRank(n));
+ return index[n];
+
+ }
+ DeclException1(ExcRank, unsigned,
+ << "Index " << arg1 << " higher than maximum 3");
+};
+
+template<>
+class TensorIndex<3>
+{
+private:
+ /**
+ * Field of indices.
+ */
+ unsigned index[3];
+public:
+ /**
+ * Constructor taking #rank# indices.
+ */
+ TensorIndex(unsigned i0, unsigned i1, unsigned i2)
+ {
+ index[0] = i0;
+ index[1] = i1;
+ index[2] = i2;
+ }
+
+
+ /**
+ *
+ * Access operator returning index
+ * in #n#th component
+ *
+ */
+ unsigned operator () (unsigned n) const
+ {
+ Assert(n<3, ExcRank(n));
+ return index[n];
+
+ }
+
+ DeclException1(ExcRank, unsigned,
+ << "Index " << arg1 << " higher than maximum 2");
+};
+
+template<>
+class TensorIndex<2>
+{
+private:
+ /**
+ * Field of indices.
+ */
+ unsigned index[4];
+public:
+ /**
+ * Constructor taking #rank# indices.
+ */
+ TensorIndex(unsigned i0, unsigned i1)
+ {
+ index[0] = i0;
+ index[1] = i1;
+ }
+
+
+ /**
+ *
+ * Access operator returning index
+ * in #n#th component
+ *
+ */
+ unsigned operator () (unsigned n) const
+ {
+ Assert(n<2, ExcRank(n));
+ return index[n];
+
+ }
+ DeclException1(ExcRank, unsigned,
+ << "Index " << arg1 << " higher than maximum 1");
+};
+
+template<>
+class TensorIndex<1>
+{
+private:
+ /**
+ * Field of indices.
+ */
+ unsigned index[1];
+public:
+ /**
+ * Constructor taking #rank# indices.
+ */
+ TensorIndex(unsigned i0)
+ {
+ index[0] = i0;
+ }
+
+
+ /**
+ *
+ * Access operator returning index
+ * in #n#th component
+ *
+ */
+ unsigned operator () (unsigned n) const
+ {
+ Assert(n<1, ExcRank(n));
+ return index[n];
+
+ }
+ DeclException1(ExcRank, unsigned,
+ << "Index " << arg1 << " higher than maximum 0");
+};
+
+template<int rank>
+inline unsigned
+TensorIndex<rank>::operator() (unsigned n) const
+{
+ Assert(n<rank, ExcRank(n));
+ return index[n];
+}
+
+
+
+/*---------------------------- tensorindex.h ---------------------------*/
+/* end of #ifndef __tensorindex_H */
+#endif
+/*---------------------------- tensorindex.h ---------------------------*/
--- /dev/null
+// $Id$
+
+
+#include <base/tensorfunction.h>
+#include <vector>
+#include <base/tensor.h>
+
+
+template <int rank_, int dim>
+TensorFunction<rank_, dim>::TensorFunction (const double initial_time) :
+ FunctionTime(initial_time)
+{};
+
+
+
+template <int rank_, int dim>
+TensorFunction<rank_, dim>::~TensorFunction ()
+{};
+
+
+
+template <int rank_, int dim>
+double
+TensorFunction<rank_, dim>::operator () (TensorIndex<rank_> i,
+ const Point<dim> &) const
+{
+ int k=i(0);
+ k++;
+
+ Assert (false, ExcPureFunctionCalled());
+ return 0;
+};
+
+
+
+template <int rank_, int dim>
+void
+TensorFunction<rank_, dim>::value_list (const vector<Point<dim> > &points,
+ vector<Tensor<rank_,dim> > &values) const
+{
+ Assert (values.size() == points.size(),
+ ExcVectorHasWrongSize(values.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ values[i] = this->operator() (points[i]);
+};
+
+
+
+
+template <int rank_, int dim>
+Tensor<rank_+1,dim>
+TensorFunction<rank_, dim>::gradient (const Point<dim> &) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return Tensor<rank_+1,dim>();
+};
+
+
+
+template <int rank_, int dim>
+void
+TensorFunction<rank_, dim>::gradient_list (const vector<Point<dim> > &points,
+ vector<Tensor<rank_+1,dim> > &gradients) const
+{
+ Assert (gradients.size() == points.size(),
+ ExcVectorHasWrongSize(gradients.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ gradients[i] = gradient(points[i]);
+};
+
+template class TensorFunction<1,1>;
+template class TensorFunction<2,1>;
+template class TensorFunction<3,1>;
+template class TensorFunction<4,1>;
+template class TensorFunction<1,2>;
+template class TensorFunction<2,2>;
+template class TensorFunction<3,2>;
+template class TensorFunction<4,2>;
+template class TensorFunction<1,3>;
+template class TensorFunction<2,3>;
+template class TensorFunction<3,3>;
+template class TensorFunction<4,3>;