* to zero.
*/
explicit Point ();
+
+ /**
+ * Convert a tensor to a point. Since no
+ * additional data is inside a point,
+ * this is ok.
+ */
+ Point (const Tensor<1,dim> &);
/**
* Constructor for one dimensional points. This
*/
Point<dim> operator * (const double) const;
+ /**
+ * Returns the scalar product of two vectors.
+ */
+ double operator * (const Point<dim> &) const;
+
/**
* Divide by a factor. If possible, use
* #operator /=# instead since this does not
*/
Point<dim> operator / (const double) const;
- /**
- * Add another vector, i.e. move this point by
- * the given offset.
- */
- Point<dim> & operator += (const Point<dim> &);
- /**
- * Subtract another vector.
- */
- Point<dim> & operator -= (const Point<dim> &);
-
- /**
- * Scale the vector by #factor#, i.e. multiply
- * all coordinates by #factor#.
- */
- Point<dim> & operator *= (const double &factor);
-
- /**
- * Scale the vector by #1/factor#.
- */
- Point<dim> & operator /= (const double &factor);
-
- /**
- * Returns the scalar product of two vectors.
- */
- double operator * (const Point<dim> &) const;
-
/**
* Returns the scalar product of this point
* vector with itself, i.e. the square, or
Tensor<1,dim>() {};
+template <int dim>
+inline
+Point<dim>::Point (const Tensor<1,dim> &t) :
+ Tensor<1,dim>(t) {};
+
+
template <>
inline
template <int dim>
inline
-Point<dim> operator * (const double factor, const Point<dim> &p) {
- return p*factor;
-};
-
-
-
-template <int dim>
-inline
-Point<dim> Point<dim>::operator / (const double factor) const {
- return (Point<dim>(*this) /= factor);
-};
-
-
-
-template <int dim>
-inline
-Point<dim> & Point<dim>::operator += (const Point<dim> &p) {
- for (unsigned int i=0; i<dim; ++i)
- values[i] += p.values[i];
- return *this;
-};
-
-
-
-template <int dim>
-inline
-Point<dim> & Point<dim>::operator -= (const Point<dim> &p) {
- for (unsigned int i=0; i<dim; ++i)
- values[i] -= p.values[i];
- return *this;
-};
-
-
-
-template <int dim>
-inline
-Point<dim> & Point<dim>::operator *= (const double &s) {
- for (unsigned int i=0; i<dim; ++i)
- values[i] *= s;
- return *this;
+double Point<dim>::operator * (const Point<dim> &p) const {
+ // simply pass down
+ return Tensor<1,dim>::operator * (p);
};
template <int dim>
inline
-Point<dim> & Point<dim>::operator /= (const double &s) {
- for (unsigned int i=0; i<dim; ++i)
- values[i] /= s;
- return *this;
+Point<dim> operator * (const double factor, const Point<dim> &p) {
+ return p*factor;
};
template <int dim>
inline
-double Point<dim>::operator * (const Point<dim> &p) const {
- double q=0;
- for (unsigned int i=0; i<dim; ++i)
- q += values[i] * p.values[i];
- return q;
+Point<dim> Point<dim>::operator / (const double factor) const {
+ return (Point<dim>(*this) /= factor);
};
-
+
template <int dim>
inline
double Point<dim>::square () const {
Tensor & operator = (const Tensor<rank_,dim> &);
/**
- * Test for equality of two points.
+ * Test for equality of two tensors.
*/
bool operator == (const Tensor<rank_,dim> &) const;
/**
- * Test for inequality of two points.
+ * Test for inequality of two tensors.
*/
bool operator != (const Tensor<rank_,dim> &) const;
+ /**
+ * Add another vector, i.e. move this tensor by
+ * the given offset.
+ */
+ Tensor<rank_,dim> & operator += (const Tensor<rank_,dim> &);
+
+ /**
+ * Subtract another vector.
+ */
+ Tensor<rank_,dim> & operator -= (const Tensor<rank_,dim> &);
+
+ /**
+ * Scale the vector by #factor#, i.e. multiply
+ * all coordinates by #factor#.
+ */
+ Tensor<rank_,dim> & operator *= (const double &factor);
+
+ /**
+ * Scale the vector by #1/factor#.
+ */
+ Tensor<rank_,dim> & operator /= (const double &factor);
+
/**
* Reset all values to zero.
*/
inline
Tensor<rank_,dim> & Tensor<rank_,dim>::operator = (const Tensor<rank_,dim> &t) {
for (unsigned int i=0; i<dim; ++i)
- values[i] = t.values[i];
+ subtensor[i] = t.subtensor[i];
return *this;
};
inline
bool Tensor<rank_,dim>::operator == (const Tensor<rank_,dim> &p) const {
for (unsigned int i=0; i<dim; ++i)
- if (values[i] != p.values[i]) return false;
+ if (subtensor[i] != p.subtensor[i]) return false;
return true;
};
+template <int rank_, int dim>
+inline
+Tensor<rank_,dim> & Tensor<rank_,dim>::operator += (const Tensor<rank_,dim> &p) {
+ for (unsigned int i=0; i<dim; ++i)
+ subtensor[i] += p.subtensor[i];
+ return *this;
+};
+
+
+
+template <int rank_, int dim>
+inline
+Tensor<rank_,dim> & Tensor<rank_,dim>::operator -= (const Tensor<rank_,dim> &p) {
+ for (unsigned int i=0; i<dim; ++i)
+ subtensor[i] -= p.subtensor[i];
+ return *this;
+};
+
+
+
+template <int rank_, int dim>
+inline
+Tensor<rank_,dim> & Tensor<rank_,dim>::operator *= (const double &s) {
+ for (unsigned int i=0; i<dim; ++i)
+ subtensor[i] *= s;
+ return *this;
+};
+
+
+
+template <int rank_, int dim>
+inline
+Tensor<rank_,dim> & Tensor<rank_,dim>::operator /= (const double &s) {
+ for (unsigned int i=0; i<dim; ++i)
+ subtensor[i] /= s;
+ return *this;
+};
+
+
+
template <int rank_, int dim>
inline
void Tensor<rank_,dim>::clear () {
+/* Exception class. This is certainly not the best possible place for its
+ declaration, but at present, local classes to any of Tensor<> can't
+ be properly accessed (haven't investigated why). If anyone has a better
+ idea, realize it!
+*/
+DeclException1 (ExcInvalidTensorIndex,
+ int,
+ << "Invalid tensor index " << arg1);
+
+
+
+template <int dim>
+inline
+void contract (Tensor<1,dim> &dest,
+ const Tensor<2,dim> &src1,
+ const Tensor<1,dim> &src2) {
+ dest.clear ();
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ dest[i] += src1[i][j] * src2[j];
+};
+
+
template <int dim>
inline
+template <int dim>
+inline
+void contract (Tensor<2,dim> &dest,
+ const Tensor<2,dim> &src1, const unsigned int index1,
+ const Tensor<2,dim> &src2, const unsigned int index2) {
+ dest.clear ();
+
+ switch (index1)
+ {
+ case 1:
+ switch (index2)
+ {
+ case 1:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[k][i] * src2[k][j];
+ break;
+ case 2:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[k][i] * src2[j][k];
+ break;
+
+ default:
+ Assert (false, ExcInvalidTensorIndex (index2));
+ };
+ break;
+ case 2:
+ switch (index2)
+ {
+ case 1:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[i][k] * src2[k][j];
+ break;
+ case 2:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[i][k] * src2[j][k];
+ break;
+
+ default:
+ Assert (false, ExcInvalidTensorIndex (index2));
+ };
+ break;
+
+ default:
+ Assert (false, ExcInvalidTensorIndex (index1));
+ };
+};
+
+
+
+template <int dim>
+inline
+void contract (Tensor<2,dim> &dest,
+ const Tensor<3,dim> &src1, const unsigned int index1,
+ const Tensor<1,dim> &src2) {
+ dest.clear ();
+
+ switch (index1)
+ {
+ case 1:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[k][i][j] * src2[k];
+ break;
+
+ case 2:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[i][k][j] * src2[k];
+ break;
+
+ case 3:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ for (unsigned int k=0; k<dim; ++k)
+ dest[i][j] += src1[i][j][k] * src2[k];
+ break;
+
+ default:
+ Assert (false, ExcInvalidTensorIndex (index1));
+ };
+};
+
+
+
template <int dim>
inline
void contract (Tensor<3,dim> &dest,
dest[i][j][k][l] += src1[i][j][m] * src2[m][k][l];
};
-
+
+
+inline
+double determinant (const Tensor<2,2> &t) {
+ return ((t[0][0] * t[1][1]) -
+ (t[1][0] * t[0][1]));
+};
+
*/
bool operator != (const Tensor<1,dim> &) const;
+ /**
+ * Add another vector, i.e. move this point by
+ * the given offset.
+ */
+ Tensor<1,dim> & operator += (const Tensor<1,dim> &);
+ /**
+ * Subtract another vector.
+ */
+ Tensor<1,dim> & operator -= (const Tensor<1,dim> &);
+
+ /**
+ * Scale the vector by #factor#, i.e. multiply
+ * all coordinates by #factor#.
+ */
+ Tensor<1,dim> & operator *= (const double &factor);
+
+ /**
+ * Scale the vector by #1/factor#.
+ */
+ Tensor<1,dim> & operator /= (const double &factor);
+
+ /**
+ * Returns the scalar product of two vectors.
+ */
+ double operator * (const Tensor<1,dim> &) const;
+
/**
* Reset all values to zero.
*/
+template <int dim>
+inline
+Tensor<1,dim> & Tensor<1,dim>::operator += (const Tensor<1,dim> &p) {
+ for (unsigned int i=0; i<dim; ++i)
+ values[i] += p.values[i];
+ return *this;
+};
+
+
+
+template <int dim>
+inline
+Tensor<1,dim> & Tensor<1,dim>::operator -= (const Tensor<1,dim> &p) {
+ for (unsigned int i=0; i<dim; ++i)
+ values[i] -= p.values[i];
+ return *this;
+};
+
+
+
+template <int dim>
+inline
+Tensor<1,dim> & Tensor<1,dim>::operator *= (const double &s) {
+ for (unsigned int i=0; i<dim; ++i)
+ values[i] *= s;
+ return *this;
+};
+
+
+
+template <int dim>
+inline
+Tensor<1,dim> & Tensor<1,dim>::operator /= (const double &s) {
+ for (unsigned int i=0; i<dim; ++i)
+ values[i] /= s;
+ return *this;
+};
+
+
+
+template <int dim>
+inline
+double Tensor<1,dim>::operator * (const Tensor<1,dim> &p) const {
+ double q=0;
+ for (unsigned int i=0; i<dim; ++i)
+ q += values[i] * p.values[i];
+ return q;
+};
+
+
+
template <int dim>
inline
void Tensor<1,dim>::clear () {