YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
bnd_triangle.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#ifdef HAVE_CONFIG_H
6// Get the definition of the 'restrict' keyword.
7#include "config.h"
8#endif
9
10#include <math.h>
11#include <string.h>
12
13#include "utils_core.h"
14#include "geometry.h"
16
17static inline void compute_middle_point(
18 double (* restrict vertices)[3], size_t num_vertices,
19 double * restrict middle_point) {
20
21 middle_point[0] = 0;
22 middle_point[1] = 0;
23 middle_point[2] = 0;
24
25 for (size_t i = 0; i < num_vertices; ++i) {
26
27 middle_point[0] += vertices[i][0];
28 middle_point[1] += vertices[i][1];
29 middle_point[2] += vertices[i][2];
30 }
31
32 double scale = 1.0 / sqrt(middle_point[0] * middle_point[0] +
33 middle_point[1] * middle_point[1] +
34 middle_point[2] * middle_point[2]);
35
36 middle_point[0] *= scale;
37 middle_point[1] *= scale;
38 middle_point[2] *= scale;
39}
40
41static inline double * find_furthest_vertex(
42 double (* restrict vertices)[3], size_t num_vertices,
43 double * restrict middle_point) {
44
45 double * x_furthest = NULL;
46 double min_dot = 2.0; // actual co-domain is [-1;1]
47
48#ifdef __NEC__
49// vectorization of the following loop leads to failue of
50// test_interp_method_hcsbb_parallel
51// with NEC compiler when 'CFLAGS='-O1 -finline-functions'
52#pragma _NEC novector
53#endif
54 for (size_t i = 0; i < num_vertices; ++i) {
55
56 double curr_dot = vertices[i][0] * middle_point[0] +
57 vertices[i][1] * middle_point[1] +
58 vertices[i][2] * middle_point[2];
59
60 // the dot product is defined by dot = cos(a,b)*|a|*|b|
61 // since all our vectors have a length of 1, we directly get the cosine of
62 // the angle between the vectors
63 // we are only interested in determining which vertex is the furtherest,
64 // therefore we do not need the actual angle
65 if (curr_dot < min_dot) {
66 x_furthest = &(vertices[i][0]);
67 min_dot = curr_dot;
68 }
69 }
70
71 return x_furthest;
72}
73
74// computes the triangle
75static void compute_triangle(
76 double * a, double * x_middle, double * x_t, double (*x_triangle)[3]) {
77
78 // compute the three norm vectors for the planes defined by the edges and the
79 // origin
80 double edge_norm_vector[3][3];
81 for (size_t i = 0; i < 3; ++i) {
82 for (size_t j = 0; j < 3; ++j)
83 edge_norm_vector[i][j] =
84 -1.0 * a[1] * x_middle[j] + a[0] * x_t[3*i+j];
85 normalise_vector(edge_norm_vector[i]);
86 }
87
88 // the corners of the triangle are the intersections of planes
89 // they are perpendicular to the norm vectors of the adjacent edges, hence
90 // we can use the cross product to determine them
91 // we just have to be carefull with the ordering of the norm vector or else we
92 // might get a vector opposite of the actual corner of the triangle
93 crossproduct_d(edge_norm_vector[0], edge_norm_vector[1], &(x_triangle[0][0]));
94 crossproduct_d(edge_norm_vector[1], edge_norm_vector[2], &(x_triangle[1][0]));
95 crossproduct_d(edge_norm_vector[2], edge_norm_vector[0], &(x_triangle[2][0]));
96 for (size_t i = 0; i < 3; ++i) normalise_vector(&(x_triangle[i][0]));
97}
98
100 double * vertices, size_t num_vertices, double triangle[][3],
101 size_t num_tests) {
102
103 // the middle point of the triangle that is to contain all vertices
104 double x_middle[3];
105 compute_middle_point((double(*)[3])vertices, num_vertices, x_middle);
106
107 // find the vertex that is furthest away from the middle point
108 double * x_furthest =
109 find_furthest_vertex((double(*)[3])vertices, num_vertices, x_middle);
110
111 // a triangle whose edge goes through the furthest index and whose edge
112 // through the furthest index is perpendicular to the plane defined by
113 // x_middle, x_furthest and x_origin will always contain all points
114
115 // now we compute three vectors x_tj (j = 0..2), with the following properties
116 // I. all x_tj are orthogonal to x_middle
117 // II. x_t0 is in the plane (x_middle, x_furthest, x_origin)
118 // III. the angle between x_t0 and x_t1/x_t2 is 120�/-120�
119 double x_t[3][3];
120 {
121 double t[3];
122 crossproduct_d(x_middle, x_furthest, t);
123 crossproduct_d(t, x_middle, x_t[0]);
124 }
125 normalise_vector(x_t[0]);
126 rotate_vector(x_middle, (2.0 * M_PI) / 3.0, x_t[0], x_t[1]);
127 rotate_vector(x_middle, -(2.0 * M_PI) / 3.0, x_t[0], x_t[2]);
128
129 // we define the opening angle of the triangle to be the angle between
130 // x_middle and the vector to the middle of the edge of the triangle
131
132 double min_cos_angle = -2.0; // is an invalid value
133 double best_x_t[3][3] = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
134 double best_a[2] = {0.0, 0.0};
135
136 double * a = xmalloc(3 * num_vertices * 3 * sizeof(*a));
137
138 // we have one solution that should work. now we check whether different
139 // rotations of the current solution give better results
140 // we only need to check a rotation range of ]0�,60�[
141 double d_rot_angle = M_PI / (double)(3 * num_tests);
142 for (size_t test_idx = 0; test_idx < num_tests; ++test_idx) {
143
144 // current rotation angle of triangle to be checked
145 double curr_rot_angle = (double)test_idx * d_rot_angle;
146
147 double x_t_[3][3];
148
149 // for the three edges
150 for (size_t i = 0; i < 3; ++i) {
151
152 memcpy(a + i*3*num_vertices, vertices, 3*num_vertices * sizeof(*a));
153
154 // rotation of x_ti gives x_ti_
155 rotate_vector(x_middle, curr_rot_angle, x_t[i], x_t_[i]);
156
157 double A[3][3];
158 // x_o is orthogonal to x_ti_ and x_middle
159 crossproduct_d(x_middle, x_t_[i], &A[2][0]/*x_o*/);
160
161 // the vectors x_middle, x_ti_ and x_o define a coordinate system
162 // each vector can be define as
163 // v = a0*x_middle + a1*x_ti_ + a2*x_o
164
165 // the sine and cosine of the opening angle required for the current
166 // triangle to contain a vertex can be computed by:
167 // cos_angle = a0 / sqrt(a0*a0 + a1*a1)
168 // sin_angle = a1 / sqrt(a0*a0 + a1*a1)
169
170 // to get a0 and a1 we have to solve the following linear system:
171 // A*a = v
172 // where: A is a matrix consisting of (x_middle, x_ti_, x_o) and
173 // a is the vector (a0, a1, a2)
174
175 lapack_int n = 3, nrhs = (lapack_int) num_vertices, lda = n, ldx = n, ipiv[3];
176 A[0][0] = x_middle[0], A[0][1] = x_middle[1], A[0][2] = x_middle[2];
177 A[1][0] = x_t_[i][0], A[1][1] = x_t_[i][1], A[1][2] = x_t_[i][2];
178
179 // we use LAPACK to solve the linear system for all vertices at once
180 // initially a contains all vertices, after the call to LAPACKE_dgesv it
181 // contains the results
183 !LAPACKE_dgesv(
184 LAPACK_COL_MAJOR, n, nrhs, &A[0][0], lda,
185 &ipiv[0], a + i * 3 * num_vertices, ldx),
186 "ERROR: internal error (could not solve linear 3x3 system)")
187 }
188
189 // for each test_idx we actually check two triangles, which have a rotation
190 // angle of 180� between themselves
191 // each edge of the two triangles has an opposite counterpart. a vertex will
192 // influence the required opening angle for one of the two edges. the
193 // effected edge is determined base on sign of the sine of the required
194 // opening angle (see flag)
195 double curr_min_cos_angle[2] = {2.0, 2.0};
196 double * curr_x_furthest[2] = {NULL, NULL};
197 double curr_best_a[2][2];
198
199 // for all three edges and each vertex
200 for (size_t i = 0; i < 3 * num_vertices; ++i) {
201
202 int flag = a[1+3*i] < 0.0;
203 double scale = 1.0 / sqrt(a[0+3*i]*a[0+3*i] + a[1+3*i]*a[1+3*i]);
204 double cos_angle = a[0+3*i] * scale;
205
206 if (cos_angle < curr_min_cos_angle[flag]) {
207 curr_min_cos_angle[flag] = cos_angle;
208 curr_x_furthest[flag] = vertices + 3 * (i%num_vertices);
209 curr_best_a[flag][0] = cos_angle;
210 curr_best_a[flag][1] = a[1+3*i]*scale;
211 }
212 }
213
214 int flag = curr_min_cos_angle[0] < curr_min_cos_angle[1];
215
216#ifdef DEBUG
217 // debug output
218 if (curr_x_furthest[flag] != NULL) {
219 double x_triangle[3][3];
220 compute_triangle(&curr_best_a[flag][0], &x_middle[0], &x_t_[0][0],
221 &x_triangle[0][0]);
222 print_bnd_triangle(vertices, num_vertices, &x_triangle[0][0], test_idx);
223 }
224#endif // DEBUG
225
226 // if the current triangle is smaller than the best solution till now
227 if ((curr_x_furthest[flag] != NULL) &&
228 (curr_min_cos_angle[flag] > min_cos_angle)) {
229 memcpy(&best_x_t[0][0], &x_t_[0][0], 9 * sizeof(x_t_[0][0]));
230 min_cos_angle = curr_min_cos_angle[flag];
231 x_furthest = curr_x_furthest[flag];
232 best_a[0] = curr_best_a[flag][0];
233 best_a[1] = curr_best_a[flag][1];
234 }
235 }
236
237 free(a);
238
240 min_cos_angle > -2.0,
241 "ERROR: internal error (could not fine any matching triangle)")
242
243 // to get the actual triangle we have to compute the intersections of the
244 // edges
245 compute_triangle(&best_a[0], &x_middle[0], &best_x_t[0][0], triangle);
246}
247
void yac_compute_bnd_triangle(double *vertices, size_t num_vertices, double triangle[][3], size_t num_tests)
static double * find_furthest_vertex(double(*restrict vertices)[3], size_t num_vertices, double *restrict middle_point)
static void compute_triangle(double *a, double *x_middle, double *x_t, double(*x_triangle)[3])
static void compute_middle_point(double(*restrict vertices)[3], size_t num_vertices, double *restrict middle_point)
static void crossproduct_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
static void rotate_vector(double axis[], double angle, double v_in[], double v_out[])
Definition geometry.h:732
static void normalise_vector(double v[])
Definition geometry.h:689
#define xmalloc(size)
Definition ppm_xfuncs.h:66
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16