My Project
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
45namespace 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
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
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 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