Incorporate crossProduct from Arithmetic

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

dune/matrix-vector/crossproduct.hh

+#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