diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt index 9a4b2c4464fb6a3ab38889dca1aeafa0c408d7cb..090d12d36b758e655f6127fbaaa269dc4d099b30 100644 --- a/dune/matrix-vector/CMakeLists.txt +++ b/dune/matrix-vector/CMakeLists.txt @@ -3,6 +3,7 @@ install(FILES addtodiagonal.hh axpy.hh componentwisematrixmap.hh + crossproduct.hh genericvectortools.hh ldlt.hh matrixtraits.hh diff --git a/dune/matrix-vector/crossproduct.hh b/dune/matrix-vector/crossproduct.hh new file mode 100644 index 0000000000000000000000000000000000000000..374375670984a6ce40f888b7e6b237e89190e3b2 --- /dev/null +++ b/dune/matrix-vector/crossproduct.hh @@ -0,0 +1,41 @@ +#ifndef DUNE_MATRIX_VECTOR_CROSSPRODUCT_HH +#define DUNE_MATRIX_VECTOR_CROSSPRODUCT_HH + +#include <dune/common/fvector.hh> + +namespace Dune { +namespace MatrixVector { + + //! Helper class for computing the cross product + template <class T, int n> + struct CrossProductHelper { + static Dune::FieldVector<T, n> crossProduct( + const Dune::FieldVector<T, n>& a, const Dune::FieldVector<T, n>& b) { + DUNE_UNUSED_PARAMETER(a); + DUNE_UNUSED_PARAMETER(b); + DUNE_THROW(Dune::Exception, "You can only call crossProduct with dim==3"); + } + }; + + //! Specialisation for n=3 + template <class T> + struct CrossProductHelper<T, 3> { + static Dune::FieldVector<T, 3> crossProduct( + const Dune::FieldVector<T, 3>& a, const Dune::FieldVector<T, 3>& b) { + Dune::FieldVector<T, 3> r; + r[0] = a[1] * b[2] - a[2] * b[1]; + r[1] = a[2] * b[0] - a[0] * b[2]; + r[2] = a[0] * b[1] - a[1] * b[0]; + return r; + } + }; + + //! Compute the cross product of two vectors. Only works for n==3 + template <class T, int n> + Dune::FieldVector<T, n> crossProduct(const Dune::FieldVector<T, n>& a, + const Dune::FieldVector<T, n>& b) { + return CrossProductHelper<T, n>::crossProduct(a, b); + } +} +} +#endif