All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Groups Pages
GeometryHelpers.hpp
1 //===========================================================================
2 //
3 // File: GeometryHelpers.hpp
4 //
5 // Created: Mon Jun 22 13:43:23 2009
6 //
7 // Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8 // Halvor M Nilsen <hnil@sintef.no>
9 // Bjørn Spjelkavik <bsp@sintef.no>
10 //
11 // $Date$
12 //
13 // $Revision$
14 //
15 //===========================================================================
16 
17 /*
18  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
19  Copyright 2009, 2010 Statoil ASA.
20 
21  This file is part of The Open Porous Media project (OPM).
22 
23  OPM is free software: you can redistribute it and/or modify
24  it under the terms of the GNU General Public License as published by
25  the Free Software Foundation, either version 3 of the License, or
26  (at your option) any later version.
27 
28  OPM is distributed in the hope that it will be useful,
29  but WITHOUT ANY WARRANTY; without even the implied warranty of
30  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31  GNU General Public License for more details.
32 
33  You should have received a copy of the GNU General Public License
34  along with OPM. If not, see <http://www.gnu.org/licenses/>.
35 */
36 
37 #ifndef OPM_GEOMETRYHELPERS_HEADER
38 #define OPM_GEOMETRYHELPERS_HEADER
39 
40 #include <cmath>
41 
42 #include <opm/grid/utility/ErrorMacros.hpp>
43 #include "Volumes.hpp"
44 
45 namespace Dune
46 {
47 
48  namespace GeometryHelpers
49  {
55  template <class Point, template <class> class Vector>
56  Point average(const Vector<Point>& points)
57  {
58  int num_points = points.size();
59  assert(num_points > 0);
60  Point pt = points[0];
61  for (int i = 1; i < num_points; ++i) {
62  pt += points[i];
63  }
64  pt /= double(num_points);
65  return pt;
66  }
67 
73  template <class Point, template <class> class Vector>
74  double polygonArea(const Vector<Point>& points,
75  const Point& centroid)
76  {
77  double tot_area = 0.0;
78  int num_points = points.size();
79  for (int i = 0; i < num_points; ++i) {
80  Point tri[3] = { centroid, points[i], points[(i+1)%num_points] };
81  tot_area += area(tri);
82  }
83  return tot_area;
84  }
85 
86 
92  template <class Point, template <class> class Vector>
93  Point polygonCentroid(const Vector<Point>& points,
94  const Point& inpoint)
95  {
96  double tot_area = 0.0;
97  Point tot_centroid(0.0);
98  int num_points = points.size();
99  for (int i = 0; i < num_points; ++i) {
100  Point tri[3] = { inpoint, points[i], points[(i+1)%num_points] };
101  double tri_area = area(tri);
102  Point tri_w_mid = (tri[0] + tri[1] + tri[2]);
103  tri_w_mid *= tri_area/3.0;
104  tot_area += tri_area;
105  tot_centroid += tri_w_mid;
106  }
107  tot_centroid /= tot_area;
108  return tot_centroid;
109  }
110 
111 
117  template <class Point, template <class> class Vector>
118  Point polygonNormal(const Vector<Point>& points,
119  const Point& centroid)
120  {
121  Point tot_normal(0.0);
122  int num_points = points.size();
123  for (int i = 0; i < num_points; ++i) {
124  Point tri[3] = { centroid, points[i], points[(i+1)%num_points] };
125  Point d0 = tri[1] - tri[0];
126  Point d1 = tri[2] - tri[0];
127  Point w_normal = cross(d0, d1);
128  w_normal *= area(tri);
129  tot_normal += w_normal;
130  }
131  tot_normal /= tot_normal.two_norm();
132  return tot_normal;
133  }
134 
135 
141  template <class Point, template <class> class Vector>
142  double polygonCellVolume(const Vector<Point>& points,
143  const Point& face_centroid,
144  const Point& cell_centroid)
145  {
146  double tot_volume = 0.0;
147  int num_points = points.size();
148  for (int i = 0; i < num_points; ++i) {
149  Point tet[4] = { cell_centroid, face_centroid, points[i], points[(i+1)%num_points] };
150  double small_volume = std::fabs(simplex_volume(tet));
151  assert(small_volume > 0);
152  tot_volume += small_volume;
153  }
154  assert(tot_volume>0);
155  return tot_volume;
156  }
157 
158 
164  template <class Point, template <class> class Vector>
165  Point polygonCellCentroid(const Vector<Point>& points,
166  const Point& face_centroid,
167  const Point& cell_centroid)
168  {
169  Point centroid(0.0);
170  double tot_volume = 0.0;
171  int num_points = points.size();
172  for (int i = 0; i < num_points; ++i) {
173  Point tet[4] = { cell_centroid, face_centroid, points[i], points[(i+1)%num_points] };
174  double small_volume = std::fabs(simplex_volume(tet));
175  assert(small_volume > 0);
176  Point small_centroid = tet[0];
177  for(int j = 1; j < 4; ++j){
178  small_centroid += tet[j];
179  }
180  small_centroid *= small_volume/4.0;
181  centroid += small_centroid;
182  tot_volume += small_volume;
183  }
184  centroid /= tot_volume;
185  assert(tot_volume>0);
186  return centroid;
187 
188  }
189 
190 
191  } // namespace GeometryHelpers
192 
193 } // namespace Dune
194 
195 
196 #endif // OPM_GEOMETRYHELPERS_HEADER
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 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