Volumes.hpp
1 //===========================================================================
2 //
3 // File: Volumes.hpp
4 //
5 // Created: Mon Jun 22 15:46:32 2009
6 //
7 // Author(s): Jan B Thomassen <jbt@sintef.no>
8 //
9 // $Date$
10 //
11 // $Revision$
12 //
13 //===========================================================================
14 
15 /*
16  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
17  Copyright 2009, 2010 Statoil ASA.
18 
19  This file is part of The Open Porous Media project (OPM).
20 
21  OPM is free software: you can redistribute it and/or modify
22  it under the terms of the GNU General Public License as published by
23  the Free Software Foundation, either version 3 of the License, or
24  (at your option) any later version.
25 
26  OPM is distributed in the hope that it will be useful,
27  but WITHOUT ANY WARRANTY; without even the implied warranty of
28  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29  GNU General Public License for more details.
30 
31  You should have received a copy of the GNU General Public License
32  along with OPM. If not, see <http://www.gnu.org/licenses/>.
33 */
34 
35 #ifndef OPM_VOLUMES_HEADER
36 #define OPM_VOLUMES_HEADER
37 
38 #include <numeric>
39 
40 // Warning suppression for Dune includes.
41 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
42 
43 #include <dune/common/version.hh>
44 #include <dune/common/math.hh>
45 #include <dune/common/fvector.hh>
46 
47 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
48 
49 namespace Dune
50 {
51 
57  template <typename T>
58  FieldVector<T, 3> cross(const FieldVector<T, 3>& a, const FieldVector<T, 3>& b)
59  {
60  FieldVector<T, 3> res;
61  res[0] = a[1]*b[2] - a[2]*b[1];
62  res[1] = a[2]*b[0] - a[0]*b[2];
63  res[2] = a[0]*b[1] - a[1]*b[0];
64  return res;
65  }
66 
72  template <class Vector>
73  typename Vector::field_type inner(const Vector& a, const Vector& b)
74  {
75  return std::inner_product(a.begin(), a.end(), b.begin(), typename Vector::field_type());
76  }
77 
81  template<typename T, template <typename, int> class Point>
82  inline T determinantOf(const Point<T, 2>* a)
83  {
84  return a[0][0] * a[1][1] - a[1][0] * a[0][1];
85  }
86 
87 
91  template<typename T, template <typename, int> class Point>
92  inline T determinantOf(const Point<T, 3>* a)
93  {
94  return
95  a[0][0] * (a[1][1] * a[2][2] - a[2][1] * a[1][2]) -
96  a[0][1] * (a[1][0] * a[2][2] - a[2][0] * a[1][2]) +
97  a[0][2] * (a[1][0] * a[2][1] - a[2][0] * a[1][1]);
98  }
99 
100 
103  template<typename T, template <typename, int> class Point, int Dim>
104  inline T simplex_volume(const Point<T, Dim>* a)
105  {
106  Point<T, Dim> tmp[Dim];
107  for (int i = 0; i < Dim; ++i) {
108  tmp[i] = a[i+1] - a[i];
109  }
110  return determinantOf(tmp) / double(Factorial<Dim>::factorial);
111  // determinant / factorial
112  }
113 
114 
117  template <typename T, template <typename, int> class Point>
118  inline T area(const Point<T, 2>* c)
119  { return simplex_volume(c); }
120 
121 
124  template <typename T, template <typename, int> class Point>
125  inline T area(const Point<T, 3>* c)
126  {
127  // Using the one-half cross product rule
128  Point<T, 3> d0 = c[1] - c[0];
129  Point<T, 3> d1 = c[2] - c[0];
130  Point<T, 3> crossprod = cross(d0,d1);
131  return 0.5 * crossprod.two_norm();
132  }
133 
134 
136  template <typename T, template <typename, int> class Point>
137  inline T volume(const Point<T, 3>* c)
138  { return simplex_volume(c); }
139 
140 
143  template <typename T, template <typename, int> class Point>
144  T signed_area(const Point<T, 3>* c, const Point<T, 3>& normal)
145  {
146  // Using the one-half cross product rule
147  Point<T, 3> d0 = c[1] - c[0];
148  Point<T, 3> d1 = c[2] - c[0];
149  Point<T, 3> crossprod = cross(d0, d1);
150  if (inner(crossprod, normal) > 0) {
151  return 0.5 * crossprod.two_norm();
152  } else {
153  return -0.5 * crossprod.two_norm();
154  }
155  }
156 
157 
158 } // namespace Dune
159 
160 
161 
162 #endif // OPM_VOLUMES_HEADER
Holds the implementation of the CpGrid as a pimple.
Definition: OpmParserIncludes.hpp:42
T determinantOf(const Point< T, 2 > *a)
Calculates the determinant of a 2 x 2 matrix, represented in memory as an array of two-dimensional po...
Definition: Volumes.hpp:82
Vector::field_type inner(const Vector &a, const Vector &b)
Definition: Volumes.hpp:73
T simplex_volume(const Point< T, Dim > *a)
Computes the volume of a simplex consisting of (Dim+1) vertices embedded in Euclidean space of dimens...
Definition: Volumes.hpp:104
T volume(const Point< T, 3 > *c)
Computes the volume of a 3D simplex (embedded i 3D space).
Definition: Volumes.hpp:137
T area(const Point< T, 2 > *c)
Computes the area of a 2-dimensional triangle.
Definition: Volumes.hpp:118
FieldVector< T, 3 > cross(const FieldVector< T, 3 > &a, const FieldVector< T, 3 > &b)
Definition: Volumes.hpp:58
T signed_area(const Point< T, 3 > *c, const Point< T, 3 > &normal)
Computes the signed area of a triangle embedded in 3D space.
Definition: Volumes.hpp:144