crossproduct.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=8 sw=2 sts=2:
#ifndef DUNE_PDELAB_COMMON_CROSSPRODUCT_HH
#define DUNE_PDELAB_COMMON_CROSSPRODUCT_HH
 
#include <dune/common/fvector.hh>
 
namespace Dune {
  namespace PDELab {
 
    //
    //  Cross product
    //
 
 
    template<unsigned dimB_, unsigned dimC_>
    class CrossProduct {
      static_assert(AlwaysFalse<CrossProduct>::value,
                    "CrossProduct cannot be used unspecialized");
    public:
      static const unsigned dimA;
      static const unsigned dimB;
      static const unsigned dimC;
 
      template<typename AType, typename BType, typename CType>
      CrossProduct(AType& A, const BType& B, const CType& C);
    };
 
 
    template<>
    struct CrossProduct<3, 3> {
      static const unsigned dimA = 3;
      static const unsigned dimB = 3;
      static const unsigned dimC = 3;
 
      template<typename AType, typename BType, typename CType>
      CrossProduct(AType& A, const BType& B, const CType& C) {
        for(unsigned i = 0; i < 3; ++i) {
          unsigned j = (i+1)%3;
          unsigned k = (i+2)%3;
          A[i] = B[j]*C[k] - B[k]*C[j];
        }
      }
    };
 
 
    template<>
    struct CrossProduct<2, 2> {
      static const unsigned dimA = 1;
      static const unsigned dimB = 2;
      static const unsigned dimC = 2;
 
      template<typename AType, typename BType, typename CType>
      CrossProduct(AType& A, const BType& B, const CType& C) {
        A[0] = B[0]*C[1] - B[1]*C[0];
      }
    };
 
 
    template<>
    struct CrossProduct<2, 1> {
      static const unsigned dimA = 2;
      static const unsigned dimB = 2;
      static const unsigned dimC = 1;
 
      template<typename AType, typename BType, typename CType>
      CrossProduct(AType& A, const BType& B, const CType& C) {
        A[0] =   B[1]*C[0];
        A[1] = - B[0]*C[0];
      }
    };
 
 
    template<>
    struct CrossProduct<1, 2> {
      static const unsigned dimA = 2;
      static const unsigned dimB = 1;
      static const unsigned dimC = 2;
 
      template<typename AType, typename BType, typename CType>
      CrossProduct(AType& A, const BType& B, const CType& C) {
        A[0] = - B[0]*C[1];
        A[1] =   B[0]*C[0];
      }
    };
 
    //
    //  Free-standing cross product functions for FieldVectors
    //
 
 
    template<typename T, int dimB, int dimC>
    void crossproduct(const FieldVector<T, dimB>& b,
                      const FieldVector<T, dimC>& c,
                      FieldVector<T, (CrossProduct<dimB, dimC>::dimA)>& a) {
      CrossProduct<dimB,dimC>(a,b,c);
    }
 
 
    template<typename T, int dimB, int dimC>
    FieldVector<T, CrossProduct<dimB, dimC>::dimA>
    crossproduct(const FieldVector<T, dimB>& b,
                 const FieldVector<T, dimC>& c) {
      FieldVector<T, CrossProduct<dimB, dimC>::dimA> a;
      crossproduct(b,c,a);
      return a;
    }
 
  } // namespace PDELab
} //namespace Dune
 
#endif // DUNE_PDELAB_COMMON_CROSSPRODUCT_HH