preCICE v3.1.1
Searching...
No Matches
geometry.cpp
Go to the documentation of this file.
1#include "geometry.hpp"
2#include <Eigen/Core>
3#include <Eigen/Geometry>
4#include <algorithm>
5#include <cmath>
6#include <cstdlib>
7#include "logging/Logger.hpp"
8#include "math/math.hpp"
9#include "utils/Helpers.hpp"
10#include "utils/assertion.hpp"
11
13
14logging::Logger _log("math::geometry");
15
17 const Eigen::Ref<const Eigen::Vector2d> &a,
18 const Eigen::Ref<const Eigen::Vector2d> &b,
19 const Eigen::Ref<const Eigen::Vector2d> &c,
20 const Eigen::Ref<const Eigen::Vector2d> &d,
21 Eigen::Ref<Eigen::Vector2d> & intersectionPoint)
22{
23 // Compute denominator for solving 2x2 equation system
24 double D = a(0) * (d(1) - c(1)) +
25 b(0) * (c(1) - d(1)) +
26 d(0) * (b(1) - a(1)) -
27 c(0) * (a(1) - b(1));
28
29 // If D==0, the two lines are parallel
30 if (math::equals(D, 0.0)) {
31 return false;
32 }
33
34 // Compute line segment parameter s, which is in [0, 1], if the
35 // intersection point is within the segment.
36 double s = (a(0) * (d(1) - c(1)) + c(0) * (a(1) - d(1)) + d(0) * (c(1) - a(1))) / D;
37
38 intersectionPoint = a + s * (b - a);
39 return true;
40}
41
43 const Eigen::Vector3d &pointOnPlane,
44 const Eigen::Vector3d &planeNormal,
45 const Eigen::Vector3d &firstPointSegment,
46 const Eigen::Vector3d &secondPointSegment,
47 Eigen::Vector3d & intersectionPoint)
48{
49 // Methodology of "Computation Geometry in C", Joseph O'Rourke, Chapter 7.3.1
50
51 // The plane is represented by: p(x,y,z) dot planeNormal = d (1)
52 // We need to compute d by using pointOnPlane as p(x,y,z):
53 double d = pointOnPlane.dot(planeNormal);
54
55 // The segment is represented by:
56 // firstPointSegment + (secondPointSegment - firstPointSegment) * t (2)
57 // We insert (2) into (1) and reorder for t. The resulting equations reads
58 // t = (d - firstPointSegment dot planeNormal) /
59 // ((secondPointSegment - firstPointSegment) dot planeNormal). (3)
60 //
61 // If the denominator of (3) equals 0, the segment is parallel to the plane.
62 // If, in addition, the denominator of (3) equals 0, the segment lies in the
63 // plane. Otherwise we can compute a point of intersection.
64 Eigen::Vector3d segmentVec(secondPointSegment - firstPointSegment);
65 double nominator = d - firstPointSegment.dot(planeNormal);
66 double denominator = segmentVec.dot(planeNormal);
67 if (math::equals(denominator, 0.0)) {
68 if (math::equals(nominator, 0.0)) {
69 return CONTAINED;
70 } else {
71 return NO_INTERSECTION;
72 }
73 }
74 double t = nominator / denominator;
75
76 // If t is larger than 1 or smaller than zero, the intersection is not within
77 // the line segment
78 if (math::greater(t, 1.0) || math::greater(0.0, t)) {
79 return NO_INTERSECTION;
80 }
81
82 intersectionPoint = firstPointSegment + segmentVec * t;
83
84 // If t equals 1 or 0, the segment is just touching the plane, otherwise, a
85 // real intersection is happening.
86 if (math::equals(t, 0.0) || math::equals(t, 1.0)) {
88 }
89
90 return INTERSECTION;
91}
92
94 const Eigen::VectorXd &a,
95 const Eigen::VectorXd &b,
96 const Eigen::VectorXd &c)
97{
98 PRECICE_ASSERT(a.size() == b.size(), a.size(), b.size());
99 PRECICE_ASSERT(b.size() == c.size(), b.size(), c.size());
100 if (a.size() == 2) {
101 Eigen::Vector2d A = b;
102 A -= a;
103 Eigen::Vector2d B = c;
104 B -= a;
105 return 0.5 * std::fabs(A(0) * B(1) - A(1) * B(0));
106 } else {
107 PRECICE_ASSERT(a.size() == 3, a.size());
108 Eigen::Vector3d A = b;
109 A -= a;
110 Eigen::Vector3d B = c;
111 B -= a;
112 return 0.5 * A.cross(B).norm();
113 }
114}
115
117 const Eigen::Vector3d &a,
118 const Eigen::Vector3d &b,
119 const Eigen::Vector3d &c,
120 const Eigen::Vector3d &d)
121{
122 return std::abs((a - d).dot((b - d).cross(c - d))) / 6.0;
123}
124
125Eigen::Vector2d projectVector(
126 const Eigen::Vector3d &vector,
127 int indexDimensionToRemove)
128{
129 PRECICE_ASSERT(indexDimensionToRemove >= 0);
130 PRECICE_ASSERT(indexDimensionToRemove < 3);
131 Eigen::Vector2d projectedVector;
132 int reducedDim = 0;
133 for (int dim = 0; dim < 3; dim++) {
134 if (indexDimensionToRemove == dim) {
135 continue;
136 }
137 projectedVector(reducedDim) = vector(dim);
138 reducedDim++;
139 }
140 return projectedVector;
141}
142
144{
145 /*
146 All points need to be projected into a new plane with only 2 coordinates, x' and y'. These are used to check
147 the convexity of the quad. These new coordinates are stored in 'coords'.
148 */
149 PRECICE_ASSERT(std::all_of(coords.cbegin(), coords.cend(),
150 [](const auto &v) { return v.size() == 3; }),
151 "This only works in 3D.");
152
153 Eigen::Vector3d coordOrigin; // Origin point for the transformation of points onto the new plane
154 coordOrigin = coords[0];
155
156 // Normal of the plane of first three points in the list of vertices
157 Eigen::Vector3d e_1 = coords[1] - coordOrigin;
158 Eigen::Vector3d e_2 = coords[2] - coordOrigin;
159 Eigen::Vector3d normalVector = e_1.cross(e_2);
160
161 PRECICE_CHECK(math::equals(normalVector.dot(coords[3] - coordOrigin), 0.0),
162 "Non-planar quads are not supported. The vertex coordinates are: {}.", coords);
163
164 //Transform Coordinates - coord[0] is the origin
165 for (int i = 0; i < 4; i++) {
166 Eigen::Vector3d coordinateDifference = coords[i] - coordOrigin;
167 coords[i][0] = e_1.dot(coordinateDifference);
168 coords[i][1] = e_2.dot(coordinateDifference);
169 coords[i][2] = normalVector.dot(coordinateDifference);
170 }
171
172 /*
173 For the convex hull algorithm, the most left hand point regarding the x coordinate is chosen as the starting point.
174 The algorithm moves in an anti-clockwise position, finding the most right hand coordinate from the
175 previous most right hand point. The algorithm must find 3 other points in order to be a valid quad.
176 */
177
178 //First find point with smallest x coord. This point must be in the convex set then and is the starting point of gift wrapping algorithm
179 int idLowestPoint = 0;
180 for (int i = 1; i < 4; i++) {
181 if (coords[i][0] < coords[idLowestPoint][0]) {
182 idLowestPoint = i;
183 }
184 }
185
186 // Found starting point. Add this as the first vertex in the convex hull.
187 // current is the origin point => hull[0]
188 int validVertexIDCounter = 0; // Counts number of times a valid vertex is found. Must be 4 for a valid quad.
189 int currentVertex = idLowestPoint; // current valid vertex
190 ConvexityResult result{};
191 do {
192 // Add current point to result
193 result.vertexOrder[validVertexIDCounter] = currentVertex;
194
195 // Next potential valid vertex
196 int nextVertex = (currentVertex + 1) % 4; // remainder resets loop through vector of points
197 for (int i = 0; i < 4; i++) {
198 double y1 = coords[currentVertex][1] - coords[nextVertex][1];
199 double y2 = coords[currentVertex][1] - coords[i][1];
200 double x1 = coords[currentVertex][0] - coords[nextVertex][0];
201 double x2 = coords[currentVertex][0] - coords[i][0];
202 double val = y2 * x1 - y1 * x2;
203
204 if (val > 0) {
205 nextVertex = i;
206 }
207 }
208 // Now nextVertex is the most anti-clockwise with respect to current
209 // Set current as nextVertex for next iteration, so that nextVertex is added to
210 // result 'hull'
211 currentVertex = nextVertex;
212 validVertexIDCounter++;
213 } while (currentVertex != idLowestPoint); // While we don't come to first point
214
215 //Ordering of quad is hull 0-1-2-3-0
216 result.convex = (validVertexIDCounter == 4);
217
218 return result;
219}
220
221} // namespace precice::math::geometry
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
T all_of(T... args)
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T cbegin(T... args)
This class provides a lightweight logger.
Definition Logger.hpp:16
T cend(T... args)
T fabs(T... args)
Provides computational geometry operations.
Definition geometry.cpp:12
double triangleArea(const Eigen::VectorXd &a, const Eigen::VectorXd &b, const Eigen::VectorXd &c)
Computes the signed area of a triangle in 2D.
Definition geometry.cpp:93
double tetraVolume(const Eigen::Vector3d &a, const Eigen::Vector3d &b, const Eigen::Vector3d &c, const Eigen::Vector3d &d)
Computes the (unsigned) area of a triangle in 3D.
Definition geometry.cpp:116
Eigen::Vector2d projectVector(const Eigen::Vector3d &vector, int indexDimensionToRemove)
Projects a 3D vector to a 2D one by removing one dimension.
Definition geometry.cpp:125
bool lineIntersection(const Eigen::Ref< const Eigen::Vector2d > &a, const Eigen::Ref< const Eigen::Vector2d > &b, const Eigen::Ref< const Eigen::Vector2d > &c, const Eigen::Ref< const Eigen::Vector2d > &d, Eigen::Ref< Eigen::Vector2d > &intersectionPoint)
Determines the intersection point of two lines, if one exists.
Definition geometry.cpp:16
ResultConstants segmentPlaneIntersection(const Eigen::Vector3d &pointOnPlane, const Eigen::Vector3d &planeNormal, const Eigen::Vector3d &firstPointSegment, const Eigen::Vector3d &secondPointSegment, Eigen::Vector3d &intersectionPoint)
Determines the intersection point of a segment with a plane in 3D.
Definition geometry.cpp:42
logging::Logger _log("math::geometry")
ConvexityResult isConvexQuad(std::array< Eigen::VectorXd, 4 > coords)
Definition geometry.cpp:143
constexpr bool equals(const Eigen::MatrixBase< DerivedA > &A, const Eigen::MatrixBase< DerivedB > &B, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares two Eigen::MatrixBase for equality up to tolerance.
std::enable_if< std::is_arithmetic< Scalar >::value, bool >::type greater(Scalar A, Scalar B, Scalar tolerance=NUMERICAL_ZERO_DIFFERENCE)