37 #ifndef OPM_GEOMETRYHELPERS_HEADER 38 #define OPM_GEOMETRYHELPERS_HEADER 42 #include <opm/grid/utility/ErrorMacros.hpp> 43 #include "Volumes.hpp" 48 namespace GeometryHelpers
55 template <
class Po
int,
template <
class>
class Vector>
56 Point average(
const Vector<Point>& points)
58 int num_points = points.size();
59 assert(num_points > 0);
61 for (
int i = 1; i < num_points; ++i) {
64 pt /= double(num_points);
73 template <
class Po
int,
template <
class>
class Vector>
74 double polygonArea(
const Vector<Point>& points,
75 const Point& centroid)
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);
92 template <
class Po
int,
template <
class>
class Vector>
93 Point polygonCentroid(
const Vector<Point>& points,
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;
107 tot_centroid /= tot_area;
117 template <
class Po
int,
template <
class>
class Vector>
118 Point polygonNormal(
const Vector<Point>& points,
119 const Point& centroid)
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;
131 tot_normal /= tot_normal.two_norm();
141 template <
class Po
int,
template <
class>
class Vector>
142 double polygonCellVolume(
const Vector<Point>& points,
143 const Point& face_centroid,
144 const Point& cell_centroid)
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] };
151 assert(small_volume > 0);
152 tot_volume += small_volume;
154 assert(tot_volume>0);
164 template <
class Po
int,
template <
class>
class Vector>
165 Point polygonCellCentroid(
const Vector<Point>& points,
166 const Point& face_centroid,
167 const Point& cell_centroid)
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] };
175 assert(small_volume > 0);
176 Point small_centroid = tet[0];
177 for(
int j = 1; j < 4; ++j){
178 small_centroid += tet[j];
180 small_centroid *= small_volume/4.0;
181 centroid += small_centroid;
182 tot_volume += small_volume;
184 centroid /= tot_volume;
185 assert(tot_volume>0);
196 #endif // OPM_GEOMETRYHELPERS_HEADER Holds the implementation of the CpGrid as a pimple.
Definition: OpmParserIncludes.hpp:42
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