Volumes.hpp

//===========================================================================
//
// File: Volumes.hpp
//
// Created: Mon Jun 22 15:46:32 2009
//
// Author(s): Jan B Thomassen <jbt@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================

/*
  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
  Copyright 2009, 2010 Statoil ASA.

  This file is part of The Open Porous Media project  (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_VOLUMES_HEADER
#define OPM_VOLUMES_HEADER

#include <numeric>

// Warning suppression for Dune includes.
#include <opm/grid/utility/platform_dependent/disable_warnings.h>

#include <dune/common/version.hh>
#include <dune/common/math.hh>
#include <dune/common/fvector.hh>

#include <opm/grid/utility/platform_dependent/reenable_warnings.h>

namespace Dune
{

    template <typename T>
    FieldVector<T, 3> cross(const FieldVector<T, 3>& a, const FieldVector<T, 3>& b)
    {
        FieldVector<T, 3> res;
        res[0] = a[1]*b[2] - a[2]*b[1];
        res[1] = a[2]*b[0] - a[0]*b[2];
        res[2] = a[0]*b[1] - a[1]*b[0];
        return res;
    }

    template <class Vector>
    typename Vector::field_type inner(const Vector& a, const Vector& b)
    {
        return std::inner_product(a.begin(), a.end(), b.begin(), typename Vector::field_type());
    }

    template<typename T, template <typename, int> class Point>
    inline T determinantOf(const Point<T, 2>* a)
    {
        return a[0][0] * a[1][1] - a[1][0] * a[0][1];
    }


    template<typename T, template <typename, int> class Point>
    inline T determinantOf(const Point<T, 3>* a)
    {
        return
            a[0][0] * (a[1][1] * a[2][2] - a[2][1] * a[1][2]) -
            a[0][1] * (a[1][0] * a[2][2] - a[2][0] * a[1][2]) +
            a[0][2] * (a[1][0] * a[2][1] - a[2][0] * a[1][1]);
    }


    template<typename T, template <typename, int> class Point, int Dim>
    inline T simplex_volume(const Point<T, Dim>* a)
    {
        Point<T, Dim> tmp[Dim];
        for (int i = 0; i < Dim; ++i) {
            tmp[i] = a[i+1] - a[i];
        }

        return determinantOf(tmp) / double(factorial(Dim));
        // determinant / factorial
    }


    template <typename T, template <typename, int> class Point>
    inline T area(const Point<T, 2>* c)
    { return simplex_volume(c); }


    template <typename T, template <typename, int> class Point>
    inline T area(const Point<T, 3>* c)
    {
        // Using the one-half cross product rule
        Point<T, 3> d0 = c[1] - c[0];
        Point<T, 3> d1 = c[2] - c[0];
        Point<T, 3> crossprod = cross(d0,d1);
        return 0.5 * crossprod.two_norm();
    }


    template <typename T, template <typename, int> class Point>
    inline T volume(const Point<T, 3>* c)
    { return simplex_volume(c); }


    template <typename T, template <typename, int> class Point>
    T signed_area(const Point<T, 3>* c, const Point<T, 3>& normal)
    {
        // Using the one-half cross product rule
        Point<T, 3> d0 = c[1] - c[0];
        Point<T, 3> d1 = c[2] - c[0];
        Point<T, 3> crossprod = cross(d0, d1);
        if (inner(crossprod, normal) > 0) {
            return 0.5 * crossprod.two_norm();
        } else {
            return -0.5 * crossprod.two_norm();
        }
    }


} // namespace Dune



#endif // OPM_VOLUMES_HEADER