preCICE v3.1.1
Loading...
Searching...
No Matches
barycenter.cpp
Go to the documentation of this file.
1#include "math/barycenter.hpp"
2#include <Eigen/Core>
3#include <Eigen/Geometry>
4#include <utility>
6#include "math/geometry.hpp"
7#include "utils/assertion.hpp"
8
10
12 const Eigen::VectorXd &a,
13 const Eigen::VectorXd &b,
14 const Eigen::VectorXd &u)
15{
16 using Eigen::Vector2d;
17 using Eigen::VectorXd;
18
19 const int dimensions = a.size();
20 PRECICE_ASSERT(dimensions == b.size(), "A and B need to have the same dimensions.", dimensions, b.size());
21 PRECICE_ASSERT(dimensions == u.size(), "A and the point need to have the same dimensions.", dimensions, u.size());
22 PRECICE_ASSERT((dimensions == 2) || (dimensions == 3), dimensions);
23
24 Vector2d barycentricCoords;
25 VectorXd ab, au;
26
27 // Let AB be the edge and U the input point. We compute the projection P of U on the edge.
28 // To find P, start from A and move from dot(AU, AB) / |AB| along the AB direction.
29 // Divide again by |AB| to find the barycentric coordinate of P relative to U
30 // This means we just need to compute dot(AU, AB) / dot(AB, AB)
31 ab = b - a;
32 au = u - a;
33
34 barycentricCoords(1) = au.dot(ab) / ab.dot(ab);
35 barycentricCoords(0) = 1 - barycentricCoords(1);
36
37 return barycentricCoords;
38}
39
40static double crossProduct2D(const Eigen::Vector2d &u, const Eigen::Vector2d &v)
41{
42 return u(0) * v(1) - u(1) * v(0);
43}
44
46 const Eigen::VectorXd &a,
47 const Eigen::VectorXd &b,
48 const Eigen::VectorXd &c,
49 const Eigen::VectorXd &u)
50{
51 using Eigen::Vector2d;
52 using Eigen::Vector3d;
53
54 const int dimensions = a.size();
55 PRECICE_ASSERT(dimensions == b.size(), "A and B need to have the same dimensions.", dimensions, b.size());
56 PRECICE_ASSERT(dimensions == c.size(), "A and C need to have the same dimensions.", dimensions, c.size());
57 PRECICE_ASSERT(dimensions == u.size(), "A and the point need to have the same dimensions.", dimensions, u.size());
58
59 Vector3d barycentricCoords;
60 double scaleFactor;
61
62 // constant per triangle
63 if (dimensions == 3) {
64 Vector3d ab, ac, au, n;
65
66 ab = b - a;
67 ac = c - a;
68 n = ab.cross(ac);
69 auto nDotN = n.dot(n);
70 PRECICE_ASSERT(nDotN != 0, "It seems a degenerate triangle was sent.");
71 scaleFactor = 1.0 / nDotN;
72
73 // varying per point
74 au = u - a;
75
76 barycentricCoords(2) = n.dot(ab.cross(au)) * scaleFactor;
77 barycentricCoords(1) = n.dot(au.cross(ac)) * scaleFactor;
78 barycentricCoords(0) = 1 - barycentricCoords(1) - barycentricCoords(2);
79
80 } else {
81 // Dimension == 2 => No cross product
82 Vector2d ab, ac, ub, uc, ua;
83 ab = b - a;
84 ac = c - a;
85 ub = b - u;
86 uc = c - u;
87 ua = a - u;
88
89 auto twiceArea = crossProduct2D(ab, ac);
90 PRECICE_ASSERT(twiceArea != 0, "It seems a degenerate triangle was sent.");
91 scaleFactor = 1.0 / twiceArea;
92
93 barycentricCoords(0) = crossProduct2D(ub, uc) * scaleFactor;
94 barycentricCoords(1) = crossProduct2D(uc, ua) * scaleFactor;
95 barycentricCoords(2) = 1 - barycentricCoords(0) - barycentricCoords(1);
96 }
97
98 return barycentricCoords;
99}
100
102 const Eigen::VectorXd &a,
103 const Eigen::VectorXd &b,
104 const Eigen::VectorXd &c,
105 const Eigen::VectorXd &d,
106 const Eigen::VectorXd &u)
107{
108 using Eigen::Vector3d;
109 using Eigen::Vector4d;
110
111 const int dimensions = a.size();
112
113 PRECICE_ASSERT(dimensions == 3, dimensions);
114 PRECICE_ASSERT(dimensions == b.size(), "A and B need to have the same dimensions.", dimensions, b.size());
115 PRECICE_ASSERT(dimensions == c.size(), "A and C need to have the same dimensions.", dimensions, c.size());
116 PRECICE_ASSERT(dimensions == d.size(), "A and D need to have the same dimensions.", dimensions, d.size());
117 PRECICE_ASSERT(dimensions == u.size(), "A and the point need to have the same dimensions.", dimensions, u.size());
118
119 Vector4d barycentricCoords;
120
121 // Varying per point
122 Vector3d au = u - a;
123 Vector3d du = u - d;
124 Vector3d cu = u - c;
125
126 // Necessary to compute triangles
127 Vector3d ab = b - a;
128 Vector3d ac = c - a;
129 Vector3d ad = d - a;
130 Vector3d bc = c - b;
131
132 // Triangles
133 Vector3d abc = ab.cross(bc);
134 Vector3d abd = ab.cross(-ad);
135 Vector3d acd = ac.cross(ad);
136
137 auto volume = abc.dot(ad);
138
139 barycentricCoords(3) = abc.dot(au) / volume;
140 barycentricCoords(2) = abd.dot(du) / volume;
141 barycentricCoords(1) = acd.dot(cu) / volume;
142 barycentricCoords(0) = 1 - barycentricCoords(3) - barycentricCoords(2) - barycentricCoords(1);
143
144 return barycentricCoords;
145}
146
147} // namespace precice::math::barycenter
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
Provides operations to calculate barycentric coordinates for a point's projection onto a primitive.
Definition barycenter.cpp:9
static double crossProduct2D(const Eigen::Vector2d &u, const Eigen::Vector2d &v)
Eigen::Vector3d calcBarycentricCoordsForTriangle(const Eigen::VectorXd &a, const Eigen::VectorXd &b, const Eigen::VectorXd &c, const Eigen::VectorXd &u)
Eigen::Vector4d calcBarycentricCoordsForTetrahedron(const Eigen::VectorXd &a, const Eigen::VectorXd &b, const Eigen::VectorXd &c, const Eigen::VectorXd &d, const Eigen::VectorXd &u)
Eigen::Vector2d calcBarycentricCoordsForEdge(const Eigen::VectorXd &a, const Eigen::VectorXd &b, const Eigen::VectorXd &u)