37#ifndef OPM_GEOMETRYHELPERS_HEADER
38#define OPM_GEOMETRYHELPERS_HEADER
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>
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>
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;
108 if (std::abs(tot_area) > 0.0) {
109 tot_centroid /= tot_area;
112 tot_centroid = inpoint;
124 template <
class Po
int,
template <
class>
class Vector>
126 const Point& centroid)
128 Point tot_normal(0.0);
129 int num_points = points.size();
130 for (
int i = 0; i < num_points; ++i) {
131 Point tri[3] = { centroid, points[i], points[(i+1)%num_points] };
132 Point d0 = tri[1] - tri[0];
133 Point d1 = tri[2] - tri[0];
134 Point w_normal =
cross(d0, d1);
135 w_normal *=
area(tri);
136 tot_normal += w_normal;
139 if (
const auto length = tot_normal.two_norm(); length > 0.0) {
140 tot_normal /= length;
152 template <
class Po
int,
template <
class>
class Vector>
154 const Point& face_centroid,
155 const Point& cell_centroid)
157 double tot_volume = 0.0;
158 int num_points = points.size();
159 for (
int i = 0; i < num_points; ++i) {
160 Point tet[4] = { cell_centroid, face_centroid, points[i], points[(i+1)%num_points] };
162 assert(small_volume >= 0);
163 tot_volume += small_volume;
165 assert(tot_volume>=0);
175 template <
class Po
int,
template <
class>
class Vector>
177 const Point& face_centroid,
178 const Point& cell_centroid)
181 double tot_volume = 0.0;
182 int num_points = points.size();
183 for (
int i = 0; i < num_points; ++i) {
184 Point tet[4] = { cell_centroid, face_centroid, points[i], points[(i+1)%num_points] };
186 assert(small_volume >= 0);
187 Point small_centroid = tet[0];
188 for(
int j = 1; j < 4; ++j){
189 small_centroid += tet[j];
191 small_centroid *= small_volume/4.0;
192 centroid += small_centroid;
193 tot_volume += small_volume;
195 centroid /= tot_volume;
196 assert(tot_volume>=0);
Point polygonCentroid(const Vector< Point > &points, const Point &inpoint)
Definition: GeometryHelpers.hpp:93
double polygonArea(const Vector< Point > &points, const Point ¢roid)
Definition: GeometryHelpers.hpp:74
Point average(const Vector< Point > &points)
Definition: GeometryHelpers.hpp:56
Point polygonCellCentroid(const Vector< Point > &points, const Point &face_centroid, const Point &cell_centroid)
Definition: GeometryHelpers.hpp:176
double polygonCellVolume(const Vector< Point > &points, const Point &face_centroid, const Point &cell_centroid)
Definition: GeometryHelpers.hpp:153
Point polygonNormal(const Vector< Point > &points, const Point ¢roid)
Definition: GeometryHelpers.hpp:125
The namespace Dune is the main namespace for all Dune code.
Definition: common/CartesianIndexMapper.hpp:10
T area(const Point< T, 2 > *c)
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)
Definition: Volumes.hpp:104
Dune::FieldVector< double, 3 > Vector
Definition: cpgrid/GridHelpers.hpp:385