ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
triangle_functions.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19#ifndef UTILS_MATH_TRIANGLE_FUNCTIONS_HPP
20#define UTILS_MATH_TRIANGLE_FUNCTIONS_HPP
21
22#include "utils/Vector.hpp"
23
24#include <algorithm>
25#include <cmath>
26#include <numbers>
27
28namespace Utils {
29/**
30 * @brief Computes the normal vector of a triangle.
31 *
32 * The sign convention is such that P1P2, P1P3 and
33 * the normal form a right-handed system.
34 * The normal vector is not normalized, i.e. its length
35 * is arbitrary.
36 */
37inline Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2,
38 const Vector3d &P3) {
39 auto const u = P2 - P1;
40 auto const v = P3 - P1;
41
42 return vector_product(u, v);
43}
44
45/** Computes the area of triangle between vectors P1,P2,P3, by computing
46 * the cross product P1P2 x P1P3 and taking the half of its norm.
47 */
48inline double area_triangle(const Vector3d &P1, const Vector3d &P2,
49 const Vector3d &P3) {
50 return 0.5 * get_n_triangle(P1, P2, P3).norm();
51}
52
53/**
54 * @brief Computes the angle between two triangles in 3D space
55 *
56 * Returns the angle between two triangles in 3D space given by points P1P2P3
57 * and P2P3P4. Note, that the common edge is given as the second and the third
58 * argument. Here, the angle can have values from 0 to 2 * PI, depending on the
59 * orientation of the two triangles. So the angle can be convex or concave. For
60 * each triangle, an inward direction has been defined as the direction of one
61 * of the two normal vectors. Particularly, for triangle P1P2P3 it is the vector
62 * N1 = P2P1 x P2P3 and for triangle P2P3P4 it is N2 = P2P3 x P2P4. The method
63 * first computes the angle between N1 and N2, which gives always value between
64 * 0 and PI and then it checks whether this value must be corrected to a value
65 * between PI and 2 * PI.
66 *
67 * As an example, consider 4 points A,B,C,D in space given by coordinates
68 * A = [1,1,1], B = [2,1,1], C = [1,2,1], D = [1,1,2]. We want to determine
69 * the angle between triangles ABC and ACD. In case the orientation of the
70 * triangle ABC is [0,0,1] and orientation of ACD is [1,0,0], the resulting
71 * angle must be PI/2.0. To get correct results, note that the common edge is
72 * AC, and one must call the method as <tt>angle_btw_triangles(B,A,C,D)</tt>.
73 * With this call we have ensured that N1 = AB x AC (which coincides with
74 * [0,0,1]) and N2 = AC x AD (which coincides with [1,0,0]). Alternatively,
75 * if the orientations of the two triangles were the opposite, the correct
76 * call would be <tt>angle_btw_triangles(B,C,A,D)</tt> so that N1 = CB x CA
77 * and N2 = CA x CD.
78 */
79inline double angle_btw_triangles(const Vector3d &P1, const Vector3d &P2,
80 const Vector3d &P3, const Vector3d &P4) {
81 auto const normal1 = get_n_triangle(P2, P1, P3);
82 auto const normal2 = get_n_triangle(P2, P3, P4);
83 auto const denominator = std::sqrt(normal1.norm2() * normal2.norm2());
84 auto const cosine = std::clamp(normal1 * normal2 / denominator, -1., 1.);
85 // The angle between the faces (not considering
86 // the orientation, always less or equal to Pi)
87 // is equal to Pi minus angle between the normals
88 auto const phi = std::numbers::pi - std::acos(cosine);
89
90 // Now we need to determine, if the angle between two triangles is less than
91 // Pi or greater than Pi. To do this, we check if the point P4 lies in the
92 // halfspace given by triangle P1P2P3 and the normal to this triangle. If yes,
93 // we have an angle less than Pi, otherwise we have an angle greater than Pi.
94 // General equation of the plane is n_x*x + n_y*y + n_z*z + d = 0 where
95 // (n_x,n_y,n_z) is the normal to the plane.
96 // Point P1 lies in the plane, therefore d = -(n_x*P1_x + n_y*P1_y + n_z*P1_z)
97 // Point P4 lies in the halfspace given by normal iff n_x*P4_x + n_y*P4_y +
98 // n_z*P4_z + d >= 0
99 if (normal1 * P4 - normal1 * P1 < 0)
100 return 2. * std::numbers::pi - phi;
101
102 return phi;
103}
104} // namespace Utils
105
106#endif
Vector implementation and trait types for boost qvm interoperability.
T norm() const
Definition Vector.hpp:138
double angle_btw_triangles(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3, const Vector3d &P4)
Computes the angle between two triangles in 3D space.
Vector< T, 3 > vector_product(Vector< T, 3 > const &a, Vector< T, 3 > const &b)
Definition Vector.hpp:368
double area_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the area of triangle between vectors P1,P2,P3, by computing the cross product P1P2 x P1P3 an...
Vector3d get_n_triangle(const Vector3d &P1, const Vector3d &P2, const Vector3d &P3)
Computes the normal vector of a triangle.