YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
interp_method_hcsbb.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 <string.h>
11
13#include "interp_method_hcsbb.h"
15#include "geometry.h"
16#include "ensure_array_size.h"
17
18// the default values for HCSBB_GAUSS_NNN, HCSBB_D_NNN, and HCSBB_GAUSS_ORDER
19// were determined empirically
20
21#ifndef HCSBB_GAUSS_NNN
22#define HCSBB_GAUSS_NNN (4)
23#endif
24#ifndef HCSBB_D_NNN
25#define HCSBB_D_NNN (7)
26#endif
27
28#ifndef HCSBB_GAUSS_ORDER
29#define HCSBB_GAUSS_ORDER (6)
30#endif
31// static size_t const num_gauss_points[] = {1,1,3,4,6,7,12,13,16,19,25,27,33};
32// -> HCSBB_NUM_GAUSS_POINTS = num_gauss_points[HCSBB_GAUSS_ORDER]
33#if HCSBB_GAUSS_ORDER == 0
34#define HCSBB_NUM_GAUSS_POINTS (1)
35#elif HCSBB_GAUSS_ORDER == 1
36#define HCSBB_NUM_GAUSS_POINTS (1)
37#elif HCSBB_GAUSS_ORDER == 2
38#define HCSBB_NUM_GAUSS_POINTS (3)
39#elif HCSBB_GAUSS_ORDER == 3
40#define HCSBB_NUM_GAUSS_POINTS (4)
41#elif HCSBB_GAUSS_ORDER == 4
42#define HCSBB_NUM_GAUSS_POINTS (6)
43#elif HCSBB_GAUSS_ORDER == 5
44#define HCSBB_NUM_GAUSS_POINTS (7)
45#elif HCSBB_GAUSS_ORDER == 6
46#define HCSBB_NUM_GAUSS_POINTS (12)
47#elif HCSBB_GAUSS_ORDER == 7
48#define HCSBB_NUM_GAUSS_POINTS (13)
49#elif HCSBB_GAUSS_ORDER == 8
50#define HCSBB_NUM_GAUSS_POINTS (16)
51#elif HCSBB_GAUSS_ORDER == 9
52#define HCSBB_NUM_GAUSS_POINTS (19)
53#else
54#define HCSBB_NUM_GAUSS_POINTS (-1)
55#error "interpolation_method_hcsbb: the specified gauss order is currently not supported."
56#endif
57
58// number spherical Bernstein-Bezier coefficients
59// n_c = = ((degree + 2) * (degree + 1)) / 2
60// sbb_degree = 3
61#define HCSBB_NUM_SBB_COEFFS (10)
62
63// This tolerence needs to be lower than the
64// one being used to find matching cells.
65#define SB_COORD_TOL (1.0e-6)
66
67static size_t do_search_hcsbb(struct interp_method * method,
68 struct yac_interp_grid * interp_grid,
69 size_t * tgt_points, size_t count,
70 struct yac_interp_weights * weights);
71static void delete_hcsbb(struct interp_method * method);
72
73static struct interp_method_vtable
77
83
88
89 double bnd_triangle[3][3];
90
91 size_t * src_points;
93
98 double * weights;
99};
100
101struct weight_vector_data {
102 double weight;
103 size_t point_id;
104};
105
110
112
114 size_t n;
115};
116
118 ON_VERTEX = 0,
121 ON_TRIANGLE = 2,
123};
124
135
137
144 struct weight_vector c[2]; // c_210, c_120;
145
148
150 struct {
152 double coordinate_xyz[3];
156
160
163 double g[3];
164};
165
167
170
173};
174
176
177 size_t tgt_point;
178 double coordinate_xyz[3];
179
181 size_t src_points[3]; // sorted by global ids
182 double sb_coords[3];
183
186 union {
189 size_t idx;
191};
192
194 double coords[3][3]; // this need to be the first element in this struct
196};
197
198// compute the spherical barycentric coordinates of the given vertices with
199// respect to the given triangle
201 double * sb_coords, size_t num_vertices, double triangle[3][3]) {
202
203 double A[3][3];
204 lapack_int n = 3, nrhs = (lapack_int) num_vertices, lda = n, ldx = n, ipiv[3];
205 memcpy(A, triangle, sizeof(A));
206
207 // for a vertex v the spherical barycentric coordinates b are defined as
208 // follows
209 // A * b = v
210 // where: A is the matrix consisting of the vertex coordinates of the three
211 // corners of the triangle
212 // we compute b by solving this linear system using LAPACK
214 !LAPACKE_dgesv(
215 LAPACK_COL_MAJOR, n, nrhs, &A[0][0], lda, ipiv, sb_coords, ldx),
216 "ERROR: internal error (could not solve linear 3x3 system)")
217}
218
219static inline int compare_size_t(const void * a, const void * b) {
220
221 size_t const * a_ = a, * b_ = b;
222
223 return (*a_ > *b_) - (*b_ > *a_);
224}
225
226static inline int compare_2_size_t(const void * a, const void * b) {
227
228 size_t const * a_ = a, * b_ = b;
229 int ret;
230 if ((ret = (a_[0] > b_[0]) - (b_[0] > a_[0]))) return ret;
231 else return (a_[1] > b_[1]) - (b_[1] > a_[1]);
232}
233
234static inline int compare_3_size_t(const void * a, const void * b) {
235
236 size_t const * a_ = a, * b_ = b;
237 int ret;
238 if ((ret = (a_[0] > b_[0]) - (b_[0] > a_[0]))) return ret;
239 if ((ret = (a_[1] > b_[1]) - (b_[1] > a_[1]))) return ret;
240 else return (a_[2] > b_[2]) - (b_[2] > a_[2]);
241}
242
243static inline int compare_HCSBB_D_NNN_size_t(const void * a, const void * b) {
244
245 size_t const * a_ = a, * b_ = b;
246 int ret;
247 for (size_t i = 0; i < HCSBB_D_NNN; ++i)
248 if ((ret = (a_[i] > b_[i]) - (b_[i] > a_[i])))
249 return ret;
250 return 0;
251}
252
253// Go through the source field triangles for all target points and check for
254// duplicated vertices, edges, and triangles. Some computation for each
255// require source field vertex, edge, and triangle can be done independent of
256// the actual target point, and hence needs to be done only once.
258 struct tgt_point_search_data * tgt_point_data, size_t num_edges,
259 size_t num_triangles,
260 size_t **unique_vertices, size_t * num_unique_vertices,
261 size_t (**unique_edge_to_unique_vertex)[2], size_t * num_unique_edges,
262 size_t (**unique_triangle_to_unique_edge)[3], size_t * num_unique_triangles) {
263
264 size_t (*edge_vertices)[2] =
265 xmalloc((num_edges + 3 * num_triangles) * sizeof(*edge_vertices));
266 size_t (*triangle_vertices)[3] =
267 xmalloc(num_triangles * sizeof(*triangle_vertices));
268
269 // get all edges and triangles
270 for (size_t i = 0; i < num_edges; ++i, ++tgt_point_data)
271 memcpy(edge_vertices[i], tgt_point_data->src_points, 2 * sizeof(size_t));
272 for (size_t i = 0; i < num_triangles; ++i, ++tgt_point_data) {
273 edge_vertices[3*i+num_edges+0][0] = tgt_point_data->src_points[0];
274 edge_vertices[3*i+num_edges+0][1] = tgt_point_data->src_points[1];
275 edge_vertices[3*i+num_edges+1][0] = tgt_point_data->src_points[0];
276 edge_vertices[3*i+num_edges+1][1] = tgt_point_data->src_points[2];
277 edge_vertices[3*i+num_edges+2][0] = tgt_point_data->src_points[1];
278 edge_vertices[3*i+num_edges+2][1] = tgt_point_data->src_points[2];
279 memcpy(
280 triangle_vertices[i], tgt_point_data->src_points, 3 * sizeof(size_t));
281 }
282
283 // sort edge end triangle data
284 qsort(edge_vertices, num_edges + 3 * num_triangles,
285 sizeof(*edge_vertices), compare_2_size_t);
286 qsort(triangle_vertices, num_triangles,
287 sizeof(*triangle_vertices), compare_3_size_t);
288
289 *num_unique_edges = num_edges + 3 * num_triangles;
290 *num_unique_triangles = num_triangles;
291
292 // remove duplicated edges and vertices
293 yac_remove_duplicates_size_t_2(edge_vertices, num_unique_edges);
294 yac_remove_duplicates_size_t_3(triangle_vertices, num_unique_triangles);
295
296 // get all unique vertices
297 *num_unique_vertices = 2 * (*num_unique_edges);
298 size_t * vertices = xmalloc((*num_unique_vertices) * sizeof(*vertices));
299 memcpy(
300 vertices, &(edge_vertices[0][0]), *num_unique_vertices * sizeof(*vertices));
301 qsort(vertices, *num_unique_vertices, sizeof(*vertices), compare_size_t);
302 yac_remove_duplicates_size_t(vertices, num_unique_vertices);
303 *unique_vertices =
304 xrealloc(vertices, *num_unique_vertices * sizeof(*vertices));
305
306 size_t (*triangle_edges)[3] =
307 xmalloc(3 * (*num_unique_triangles) * sizeof(*triangle_edges));
308 for (size_t i = 0; i < *num_unique_triangles; ++i) {
309 // edge across first vertex
310 triangle_edges[3*i+0][0] = triangle_vertices[i][1];
311 triangle_edges[3*i+0][1] = triangle_vertices[i][2];
312 triangle_edges[3*i+0][2] = 3*i+0;
313 // edge across second vertex
314 triangle_edges[3*i+1][0] = triangle_vertices[i][0];
315 triangle_edges[3*i+1][1] = triangle_vertices[i][2];
316 triangle_edges[3*i+1][2] = 3*i+1;
317 // edge across third vertex
318 triangle_edges[3*i+2][0] = triangle_vertices[i][0];
319 triangle_edges[3*i+2][1] = triangle_vertices[i][1];
320 triangle_edges[3*i+2][2] = 3*i+2;
321 }
322
323 // sort triangle edges
324 qsort(triangle_edges, 3 * (*num_unique_triangles),
325 sizeof(*triangle_edges), compare_2_size_t);
326
327 // look up triangle edges in the unique edges
328 *unique_triangle_to_unique_edge =
329 xmalloc(*num_unique_triangles * sizeof(**unique_triangle_to_unique_edge));
330 for (size_t i = 0, j = 0; i < 3 * (*num_unique_triangles); ++i) {
331 size_t * curr_edge = triangle_edges[i];
332 while (compare_2_size_t(curr_edge, edge_vertices[j])) ++j;
333 (&(*unique_triangle_to_unique_edge)[0][0])[curr_edge[2]] = j;
334 }
335 free(triangle_edges);
336
337 tgt_point_data -= num_edges + num_triangles;
338 for (size_t i = 0, j = 0; i < num_edges; ++i, ++tgt_point_data) {
339 while (compare_2_size_t(edge_vertices[j], tgt_point_data->src_points)) ++j;
341 j < *num_unique_edges, "ERROR(get_unique_interp_data): edge ids missmatch")
342 tgt_point_data->interp_data.idx = j;
343 }
344 for (size_t i = 0, j = 0; i < num_triangles; ++i, ++tgt_point_data) {
345 while (compare_3_size_t(triangle_vertices[j], tgt_point_data->src_points))
346 ++j;
348 j < *num_unique_triangles,
349 "ERROR(get_unique_interp_data): triangle ids missmatch")
350 tgt_point_data->interp_data.idx = j;
351 }
352 free(triangle_vertices);
353
354 // generate unique_edge_to_unique_vertex
355 *unique_edge_to_unique_vertex =
356 xmalloc(*num_unique_edges * sizeof(**unique_edge_to_unique_vertex));
357 size_t * reorder_idx =
358 xmalloc(2 * (*num_unique_edges) * sizeof(*reorder_idx));
359 for (size_t i = 0; i < 2 * (*num_unique_edges); ++i) reorder_idx[i] = i;
361 &(edge_vertices[0][0]), 2 * (*num_unique_edges), reorder_idx);
362 for (size_t i = 0, j = 0; i < 2 * (*num_unique_edges); ++i) {
363 size_t curr_vertex = (&(edge_vertices[0][0]))[i];
364 while (curr_vertex != (*unique_vertices)[j]) ++j;
365 (&(*unique_edge_to_unique_vertex[0][0]))[reorder_idx[i]] = j;
366 }
367 free(edge_vertices);
368 free(reorder_idx);
369}
370
373 double (*gauss_points)[3], size_t * nnn_search_results,
374 double w_nnn[][HCSBB_NUM_GAUSS_POINTS], size_t * src_point_buffer,
375 yac_int * global_id_buffer, yac_const_coordinate_pointer src_point_coords,
376 const_yac_int_pointer src_point_global_ids) {
377
378 // get all unique nnn_search result points
379 memcpy(src_point_buffer, nnn_search_results,
380 HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS * sizeof(*src_point_buffer));
381 qsort(src_point_buffer, HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS,
382 sizeof(*src_point_buffer), compare_size_t);
383 size_t num_unique_result_points = HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS;
384 yac_remove_duplicates_size_t(src_point_buffer, &num_unique_result_points);
385
386 // sort result points by their global ids
387 for (size_t i = 0; i < num_unique_result_points; ++i)
388 global_id_buffer[i] = src_point_global_ids[src_point_buffer[i]];
390 global_id_buffer, num_unique_result_points, src_point_buffer);
391
393 xmalloc(num_unique_result_points * sizeof(*src_point_buffer));
394 memcpy(gauss_nnn_patch->src_points, src_point_buffer,
395 num_unique_result_points * sizeof(*src_point_buffer));
396 gauss_nnn_patch->num_src_points = num_unique_result_points;
397
398 // compute w_nnn
399 {
400 for (size_t i = 0; i < num_unique_result_points; ++i)
401 for (size_t j = 0; j < HCSBB_NUM_GAUSS_POINTS; ++j)
402 w_nnn[i][j] = 0.0;
403
404 double inv_distances[HCSBB_GAUSS_NNN];
405 size_t reorder_idx[HCSBB_GAUSS_NNN];
406 for (size_t i = 0; i < HCSBB_NUM_GAUSS_POINTS; ++i) {
407
408 size_t * curr_result_points = nnn_search_results + i * HCSBB_GAUSS_NNN;
409 double * curr_gauss_points = gauss_points[i];
410 double distance_sum = 0.0;
411 size_t j;
412
413 // for all results of the current request point
414 for (j = 0; j < HCSBB_GAUSS_NNN; ++j) {
415
416 double distance =
418 curr_gauss_points,
419 (double*)(src_point_coords[curr_result_points[j]]));
420
421 // if src and tgt point are at the same location
422 if (distance < yac_angle_tol) break;
423
424 double inv_distance = 1.0 / (distance * distance);
425 inv_distances[j] = inv_distance;
426 distance_sum += inv_distance;
427 }
428
429 // if there is no source point that exactly matches the target point
430 if (j == HCSBB_GAUSS_NNN) {
431
432 yac_int curr_global_ids[HCSBB_GAUSS_NNN];
433 for (size_t k = 0; k < HCSBB_GAUSS_NNN; ++k)
434 curr_global_ids[k] = src_point_global_ids[curr_result_points[k]];
435 for (size_t l = 0; l < HCSBB_GAUSS_NNN; ++l) reorder_idx[l] = l;
436
437 // sort current results by their global ids
439 curr_global_ids, HCSBB_GAUSS_NNN, reorder_idx);
440
441 distance_sum = 1.0 / distance_sum;
442 for (size_t k = 0, l = 0; k < HCSBB_GAUSS_NNN; ++k) {
443 yac_int curr_global_id = curr_global_ids[k];
444 while (curr_global_id != global_id_buffer[l]) ++l;
445 w_nnn[l][i] = inv_distances[reorder_idx[k]] * distance_sum;
446 }
447
448 } else {
449
450 yac_int curr_global_id = src_point_global_ids[curr_result_points[j]];
451 size_t l = 0;
452 while (curr_global_id != global_id_buffer[l]) ++l;
453 w_nnn[l][i] = 1.0;
454 }
455 }
456 }
457
458 double (*w_g)[HCSBB_NUM_GAUSS_POINTS] =
460
461 // compute and store w_g * w_nnn
462 double (*weights)[HCSBB_NUM_SBB_COEFFS] =
463 xmalloc(num_unique_result_points * sizeof(*weights));
464 gauss_nnn_patch->weights = (double*)weights;
465 for (size_t i = 0; i < num_unique_result_points; ++i) {
466 for (size_t j = 0; j < HCSBB_NUM_SBB_COEFFS; ++j) {
467 double accu = 0.0;
468 for (size_t l = 0; l < HCSBB_NUM_GAUSS_POINTS; ++l)
469 accu += w_g[j][l] * w_nnn[i][l];
470 weights[i][j] = accu;
471 }
472 }
473
474 free(w_g);
475}
476
477// computes the directional derivatives of the spherical Bernstein basis
478// polynomials for degree = 3 in the given direction
479// the polynomials B_ijk are ordered in the following order:
480// B_003, B_012, B_021, B_030, B_102, B_111, B_120, B_201, B_210, B_300
482 double bnd_triangle[3][3], double direction[3], double coordinate_xyz[3],
483 double d_sbb_polynomials[HCSBB_NUM_SBB_COEFFS]) {
484
485 double sb_coord_p[3] = {
486 coordinate_xyz[0], coordinate_xyz[1], coordinate_xyz[2]};
487 double sb_coord_d[3] = {direction[0], direction[1], direction[2]};
488
489 // compute the spherical barycentric coordinates of the point and the
490 // direction
491 compute_sb_coords(&(sb_coord_p[0]), 1, bnd_triangle);
492 compute_sb_coords(&(sb_coord_d[0]), 1, bnd_triangle);
493
494#define POW_0(x) (1.0)
495#define POW_1(x) ((x))
496#define POW_2(x) ((x)*(x))
497#define POW_3(x) ((x)*(x)*(x))
498#define FAC_0 (1.0)
499#define FAC_1 (1.0)
500#define FAC_2 (2.0)
501#define FAC_3 (6.0)
502
503 d_sbb_polynomials[0] =
504 // sb_coord_d[0] * (0.0 * FAC_3 / (FAC_0 * FAC_0 * FAC_3)) *
505 // (POW_0(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_3(sb_coord_p[2])) +
506 // sb_coord_d[1] * (0.0 * FAC_3 / (FAC_0 * FAC_0 * FAC_3)) *
507 // (POW_0(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_3(sb_coord_p[2])) +
508 sb_coord_d[2] * (3.0 * FAC_3 / (FAC_0 * FAC_0 * FAC_3)) *
509 (POW_0(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_2(sb_coord_p[2]));
510 d_sbb_polynomials[1] =
511 // sb_coord_d[0] * (0.0 * FAC_3 / (FAC_0 * FAC_1 * FAC_2)) *
512 // (POW_0(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_2(sb_coord_p[2])) +
513 sb_coord_d[1] * (1.0 * FAC_3 / (FAC_0 * FAC_1 * FAC_2)) *
514 (POW_0(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_2(sb_coord_p[2])) +
515 sb_coord_d[2] * (2.0 * FAC_3 / (FAC_0 * FAC_1 * FAC_2)) *
516 (POW_0(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_1(sb_coord_p[2]));
517 d_sbb_polynomials[2] =
518 // sb_coord_d[0] * (0.0 * FAC_3 / (FAC_0 * FAC_2 * FAC_1)) *
519 // (POW_0(sb_coord_p[0]) * POW_2(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
520 sb_coord_d[1] * (2.0 * FAC_3 / (FAC_0 * FAC_2 * FAC_1)) *
521 (POW_0(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
522 sb_coord_d[2] * (1.0 * FAC_3 / (FAC_0 * FAC_2 * FAC_1)) *
523 (POW_0(sb_coord_p[0]) * POW_2(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
524 d_sbb_polynomials[3] =
525 // sb_coord_d[0] * (0.0 * FAC_3 / (FAC_0 * FAC_3 * FAC_0)) *
526 // (POW_0(sb_coord_p[0]) * POW_3(sb_coord_p[1]) * POW_0(sb_coord_p[2])) +
527 sb_coord_d[1] * (3.0 * FAC_3 / (FAC_0 * FAC_3 * FAC_0)) *
528 (POW_0(sb_coord_p[0]) * POW_2(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
529 // sb_coord_d[2] * (0.0 * FAC_3 / (FAC_0 * FAC_3 * FAC_0)) *
530 // (POW_0(sb_coord_p[0]) * POW_3(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
531 d_sbb_polynomials[4] =
532 sb_coord_d[0] * (1.0 * FAC_3 / (FAC_1 * FAC_0 * FAC_2)) *
533 (POW_0(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_2(sb_coord_p[2])) +
534 // sb_coord_d[1] * (0.0 * FAC_3 / (FAC_1 * FAC_0 * FAC_2)) *
535 // (POW_1(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_2(sb_coord_p[2])) +
536 sb_coord_d[2] * (2.0 * FAC_3 / (FAC_1 * FAC_0 * FAC_2)) *
537 (POW_1(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_1(sb_coord_p[2]));
538 d_sbb_polynomials[5] =
539 sb_coord_d[0] * (1.0 * FAC_3 / (FAC_1 * FAC_1 * FAC_1)) *
540 (POW_0(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
541 sb_coord_d[1] * (1.0 * FAC_3 / (FAC_1 * FAC_1 * FAC_1)) *
542 (POW_1(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
543 sb_coord_d[2] * (1.0 * FAC_3 / (FAC_1 * FAC_1 * FAC_1)) *
544 (POW_1(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
545 d_sbb_polynomials[6] =
546 sb_coord_d[0] * (1.0 * FAC_3 / (FAC_1 * FAC_2 * FAC_0)) *
547 (POW_0(sb_coord_p[0]) * POW_2(sb_coord_p[1]) * POW_0(sb_coord_p[2])) +
548 sb_coord_d[1] * (2.0 * FAC_3 / (FAC_1 * FAC_2 * FAC_0)) *
549 (POW_1(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
550 // sb_coord_d[2] * (0.0 * FAC_3 / (FAC_1 * FAC_2 * FAC_0)) *
551 // (POW_1(sb_coord_p[0]) * POW_2(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
552 d_sbb_polynomials[7] =
553 sb_coord_d[0] * (2.0 * FAC_3 / (FAC_2 * FAC_0 * FAC_1)) *
554 (POW_1(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
555 // sb_coord_d[1] * (0.0 * FAC_3 / (FAC_2 * FAC_0 * FAC_1)) *
556 // (POW_2(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_1(sb_coord_p[2])) +
557 sb_coord_d[2] * (1.0 * FAC_3 / (FAC_2 * FAC_0 * FAC_1)) *
558 (POW_2(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
559 d_sbb_polynomials[8] =
560 sb_coord_d[0] * (2.0 * FAC_3 / (FAC_2 * FAC_1 * FAC_0)) *
561 (POW_1(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_0(sb_coord_p[2])) +
562 sb_coord_d[1] * (1.0 * FAC_3 / (FAC_2 * FAC_1 * FAC_0)) *
563 (POW_2(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
564 // sb_coord_d[2] * (0.0 * FAC_3 / (FAC_2 * FAC_1 * FAC_0)) *
565 // (POW_2(sb_coord_p[0]) * POW_1(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
566 d_sbb_polynomials[9] =
567 sb_coord_d[0] * (3.0 * FAC_3 / (FAC_3 * FAC_0 * FAC_0)) *
568 (POW_2(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
569 // sb_coord_d[1] * (0.0 * FAC_3 / (FAC_3 * FAC_0 * FAC_0)) *
570 // (POW_3(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_0(sb_coord_p[2])) +
571 // sb_coord_d[2] * (0.0 * FAC_3 / (FAC_3 * FAC_0 * FAC_0)) *
572 // (POW_3(sb_coord_p[0]) * POW_0(sb_coord_p[1]) * POW_0(sb_coord_p[2]));
573
574#undef POW_0
575#undef POW_1
576#undef POW_2
577#undef POW_3
578#undef FAC_0
579#undef FAC_1
580#undef FAC_2
581#undef FAC_3
582}
583
585 struct gauss_nnn_patch * gauss_nnn_patch, double coordinate_xyz[3],
586 double direction[3], struct weight_vector_data * weights,
587 double mult_factor) {
588
589 // compute w_p
590 double d_sbb_vertex_polynomials[HCSBB_NUM_SBB_COEFFS];
592 gauss_nnn_patch->bnd_triangle, direction, coordinate_xyz,
593 d_sbb_vertex_polynomials);
594
595 // compute w_p * (w_g * w_nnn)
596 size_t * src_points = gauss_nnn_patch->src_points;
597 size_t num_src_points = gauss_nnn_patch->num_src_points;
598 double * d_data_weights = gauss_nnn_patch->weights;
599 for (size_t i = 0, k = 0; i < num_src_points; ++i) {
600
601 double weight = 0.0;
602 for (size_t j = 0; j < HCSBB_NUM_SBB_COEFFS; ++j, ++k)
603 weight += d_sbb_vertex_polynomials[j] * d_data_weights[k];
604
605 weights[i].weight = weight * mult_factor;
606 weights[i].point_id = src_points[i];
607 }
608}
609
611 struct edge_interp_data * edge_data, size_t num_edges) {
612
613 size_t edge_weight_vector_data_buffer_size = 0;
614 // count total number of weight_vector_data elements required
615 for (size_t i = 0; i < num_edges; ++i)
616 edge_weight_vector_data_buffer_size +=
617 2 + edge_data[i].vertex_d_data[0]->gauss_nnn_patch->num_src_points +
619
620 // 2 -> two corner weights
621 // 2*(HCSBB_NUM_GAUSS_POINTS*HCSBB_GAUSS_NNN+1)
622 // -> maximum number of weights per edge coefficient
623 struct weight_vector_data * curr_edge_weight_vector_data =
624 (edge_weight_vector_data_buffer_size > 0)?
625 xmalloc(edge_weight_vector_data_buffer_size *
626 sizeof(*curr_edge_weight_vector_data)):NULL;
627
628 for (size_t i = 0; i < num_edges; ++i) {
629
630 struct edge_interp_data * curr_edge = edge_data + i;
631 struct vertex_interp_data * curr_vertex_d_data[2] = {
632 curr_edge->vertex_d_data[0], curr_edge->vertex_d_data[1]};
633
634 double direction[3] = {
635 curr_vertex_d_data[1]->coordinate_xyz[0] -
636 curr_vertex_d_data[0]->coordinate_xyz[0],
637 curr_vertex_d_data[1]->coordinate_xyz[1] -
638 curr_vertex_d_data[0]->coordinate_xyz[1],
639 curr_vertex_d_data[1]->coordinate_xyz[2] -
640 curr_vertex_d_data[0]->coordinate_xyz[2]};
641
642 // c_210 and c_120
643 for (size_t j = 0; j < 2; ++j) {
644
645 size_t n = curr_vertex_d_data[j]->gauss_nnn_patch->num_src_points + 1;
646
648 curr_vertex_d_data[j]->gauss_nnn_patch,
649 curr_vertex_d_data[j]->coordinate_xyz, direction,
650 curr_edge_weight_vector_data, 1.0/3.0);
651
652 curr_edge_weight_vector_data[n-1].weight = 1.0;
653 curr_edge_weight_vector_data[n-1].point_id =
654 curr_edge->vertex_d_data[j]->src_point;
655 curr_edge->c[j].data = curr_edge_weight_vector_data;
656 curr_edge->c[j].n = n;
657
658 direction[0] *= -1, direction[1] *= -1, direction[2] *= -1;
659
660 curr_edge_weight_vector_data += n;
661 }
662 }
663}
664
666 void const * a, void const * b) {
667
668 struct weight_vector_data const * weight_a =
669 ((struct weight_vector_data_pos const *)a)->weight;
670 struct weight_vector_data const * weight_b =
671 ((struct weight_vector_data_pos const *)b)->weight;
672
673 int ret = (weight_a->point_id > weight_b->point_id) -
674 (weight_a->point_id < weight_b->point_id);
675 if (ret) return ret;
676 double abs_weight_a = fabs(weight_a->weight);
677 double abs_weight_b = fabs(weight_b->weight);
678 ret = (abs_weight_a > abs_weight_b) - (abs_weight_a < abs_weight_b);
679 if (ret) return ret;
680 return (weight_a->weight > weight_b->weight) -
681 (weight_a->weight < weight_b->weight);
682}
683
685 void const * a, void const * b) {
686
687 struct weight_vector_data const * weight_a =
688 (struct weight_vector_data const *)a;
689 struct weight_vector_data const * weight_b =
690 (struct weight_vector_data const *)b;
691
692 return (weight_a->point_id > weight_b->point_id) -
693 (weight_a->point_id < weight_b->point_id);
694}
695
697 void const * a, void const * b) {
698
699 return ((struct weight_vector_data_pos const *)a)->pos -
700 ((struct weight_vector_data_pos const *)b)->pos;
701}
702
704 struct weight_vector_data * weights, size_t * n,
705 struct weight_vector_data_pos * buffer) {
706
707 size_t n_ = *n;
708
709 if (n_ <= 1) return;
710
711 for (size_t i = 0; i < n_; ++i) {
712 buffer[i].weight = weights + i;
713 buffer[i].pos = i;
714 }
715
716 qsort(buffer, n_, sizeof(*buffer), compare_weight_vector_data_pos_weight);
717
718 size_t new_n = 1;
719 struct weight_vector_data_pos * prev_weight_data_pos = buffer;
720 struct weight_vector_data * prev_weight_data = prev_weight_data_pos->weight;
721 struct weight_vector_data_pos * curr_weight_data_pos = buffer + 1;
722 for (size_t i = 1; i < n_; ++i, ++curr_weight_data_pos) {
723
724 struct weight_vector_data * curr_weight_data = curr_weight_data_pos->weight;
725 if (!compare_weight_vector_data_point(prev_weight_data, curr_weight_data)) {
726 prev_weight_data->weight += curr_weight_data->weight;
727 if (prev_weight_data_pos->pos > curr_weight_data_pos->pos)
728 prev_weight_data_pos->pos = curr_weight_data_pos->pos;
729 } else {
730 ++new_n;
731 ++prev_weight_data_pos;
732 *prev_weight_data_pos = *curr_weight_data_pos;
733 prev_weight_data = prev_weight_data_pos->weight;
734 }
735 }
736
737 if (n_ == new_n) return;
738
739 qsort(buffer, new_n, sizeof(*buffer), compare_weight_vector_data_pos_pos);
740
741 for (size_t i = 0; i < new_n; ++i) weights[i] = *(buffer[i].weight);
742
743 *n = new_n;
744}
745
747 struct triangle_interp_data * triangle_data, size_t num_triangles) {
748
749 // maximum number of coefficients per c_111 alpha value
750 // 6 edge coefficients (HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS + 1)
751 // 2 corner coefficients (1)
752 // 1 middle point coefficient (HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS)
753 struct weight_vector_data * weight_vector_data_buffer =
754 xmalloc((6 * (HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS + 1) + 2 * 1 +
756 sizeof(*weight_vector_data_buffer));
757 struct weight_vector_data_pos * weight_vector_data_pos_buffer =
758 xmalloc((6 * (HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS + 1) + 2 * 1 +
760 sizeof(*weight_vector_data_pos_buffer));
761
762 struct triangle_interp_data * curr_triangle = triangle_data;
763 for (size_t i = 0; i < num_triangles; ++i, ++curr_triangle) {
764
765 // z_w0/3 = B_020(w0) * (b_0(g0)*c_120 + b_1(g0)*c_030 + b_2(g0)*c_021) +
766 // B_011(w0) * (b_0(g0)*a_1 b_1(g0)*c_021 + b_2(g0)*c_012) +
767 // B_002(w0) * (b_0(g0)*c_102 + b_1(g0)*c_012 + b_2(g0)*c_003)
768 // w0 = ||v_1+v_2|| -> middle point of edge 0
769 // g0 = v_2 x v_1 -> vector, which is a tangent of the sphere in w0 and
770 // perpendicular to w0
771 // z_w0 = D_g0 f(w0) -> directional derivative of f in w0 in the direction
772 // of g0
773 // =>
774 // a_0 = (D_g f(w0)/3 -
775 // B_020(w0) * (b_0(g0)*c_120 + b_1(g0)*c_030 + b_2(g0)*c_021) -
776 // B_011(w0) * ( b_1(g0)*c_021 + b_2(g0)*c_012) -
777 // B_002(w0) * (b_0(g0)*c_102 + b_1(g0)*c_012 + b_2(g0)*c_003)) /
778 // (B_011(w0) * b_0(g0))
779 // a_1 = (D_g f(w1)/3 -
780 // B_200(w1) * (b_0(g1)*c_300 + b_1(g1)*c_210 + b_2(g1)*c_201) -
781 // B_101(w1) * (b_0(g1)*c_201 + b_2(g1)*c_102) -
782 // B_002(w1) * (b_0(g1)*c_102 + b_1(g1)*c_012 + b_2(g1)*c_003)) /
783 // (B_101(w1) * b_1(g1))
784 // a_2 = (D_g f(w2)/3 -
785 // B_200(w2) * (b_0(g2)*c_300 + b_1(g2)*c_210 + b_2(g2)*c_201) -
786 // B_110(w2) * (b_0(g2)*c_210 + b_1(g2)*c_120 ) -
787 // B_020(w2) * (b_0(g2)*c_120 + b_1(g2)*c_030 + b_2(g2)*c_021)) /
788 // (B_110(w2) * b_2(g2))
789
790 struct edge_interp_data ** curr_edges = curr_triangle->edges;
791 struct vertex_interp_data * vertex_d_data[3] = {
792 curr_edges[2]->vertex_d_data[0],
793 curr_edges[2]->vertex_d_data[1],
794 curr_edges[1]->vertex_d_data[1]};
795 size_t src_points[3] = {
796 vertex_d_data[0]->src_point,
797 vertex_d_data[1]->src_point,
798 vertex_d_data[2]->src_point};
799
800 // remark the corners of the triangle are sorted by global ids, the same is
801 // true for the edge data
802 struct weight_vector_data c_3_data[3] = {
803 {.weight = 1.0, .point_id = src_points[0]},
804 {.weight = 1.0, .point_id = src_points[1]},
805 {.weight = 1.0, .point_id = src_points[2]}};
806 struct weight_vector c_3[3] = {{.data = &(c_3_data[0]), .n = 1},
807 {.data = &(c_3_data[1]), .n = 1},
808 {.data = &(c_3_data[2]), .n = 1}};
809 struct weight_vector * c_300 = &(c_3[0]);
810 struct weight_vector * c_030 = &(c_3[1]);
811 struct weight_vector * c_003 = &(c_3[2]);
812 struct weight_vector * c_021 = &(curr_edges[0]->c[0]);
813 struct weight_vector * c_012 = &(curr_edges[0]->c[1]);
814 struct weight_vector * c_201 = &(curr_edges[1]->c[0]);
815 struct weight_vector * c_102 = &(curr_edges[1]->c[1]);
816 struct weight_vector * c_210 = &(curr_edges[2]->c[0]);
817 struct weight_vector * c_120 = &(curr_edges[2]->c[1]);
818
819 double corner_coordinate_xyz[3][3] = {
820 {vertex_d_data[0]->coordinate_xyz[0],
821 vertex_d_data[0]->coordinate_xyz[1],
822 vertex_d_data[0]->coordinate_xyz[2]},
823 {vertex_d_data[1]->coordinate_xyz[0],
824 vertex_d_data[1]->coordinate_xyz[1],
825 vertex_d_data[1]->coordinate_xyz[2]},
826 {vertex_d_data[2]->coordinate_xyz[0],
827 vertex_d_data[2]->coordinate_xyz[1],
828 vertex_d_data[2]->coordinate_xyz[2]}};
829
830 // vector which are tangents of the sphere in the middle points of the edges
831 // and are perpendicular to the edges
832 double * g[3] = {
833 &(curr_edges[0]->g[0]), &(curr_edges[1]->g[0]), &(curr_edges[2]->g[0])};
834
835 // spherical barycentric coordinates of g
836 double b_g[3][3] = {{g[0][0], g[0][1], g[0][2]},
837 {g[1][0], g[1][1], g[1][2]},
838 {g[2][0], g[2][1], g[2][2]}};
839 compute_sb_coords(&(b_g[0][0]), 3, corner_coordinate_xyz);
840
841 // computing quadratic spherical Bernstein-Bezier polynomials for the edge
842 // middle points (e.g. B_020(w0)) using w0 as an example
843 // I. w0 is in the middle of the edge opposing corner 0
844 // II. for the barycentric coordinates if w0 we can say the following
845 // a) b_0(w0) = 0 -> because it is on the edge opposing v0
846 // b) b_1(w0) = b_2(w0) -> because it is in the middle between v1 and v2
847 // III. B_ijk(p) = d!/(i!*j!*k!) * b_0(p)^i * b_1(p)^j * b_2(p)^k
848 // -> for all (i + j + k) = d
849 // -> p = w0 and d = 2 in our case
850 //
851 // from II.a) it follows: B_200 = B_110 = B_101 = 0
852 // from II.b) it follows: B_020 = B_002 = 2 * B_011 = b_1(w0) = b_2(w0) = bw
853
854 double bw[3] = {
855 curr_edges[0]->sb_coord_middle_point *
856 curr_edges[0]->sb_coord_middle_point,
857 curr_edges[1]->sb_coord_middle_point *
858 curr_edges[1]->sb_coord_middle_point,
859 curr_edges[2]->sb_coord_middle_point *
860 curr_edges[2]->sb_coord_middle_point};
861
862#define COPY_D_DATA(edge_idx) \
863 { \
864 get_derivative_weights( \
865 curr_edges[edge_idx]->middle_d_data.gauss_nnn_patch, \
866 curr_edges[edge_idx]->middle_d_data.coordinate_xyz, g[edge_idx], \
867 weight_vector_data_buffer + n, 1.0 / 3.0); \
868 n += curr_edges[edge_idx]->middle_d_data.gauss_nnn_patch->num_src_points; \
869 }
870#define COPY_COEFF(c_ijk, edge_idx, b_idx, bw_factor) \
871 { \
872 struct weight_vector_data * curr_buffer = weight_vector_data_buffer + n; \
873 size_t curr_n = c_ijk->n; \
874 memcpy(curr_buffer, c_ijk->data, curr_n * sizeof(*curr_buffer)); \
875 n += curr_n; \
876 double factor = - bw[edge_idx] * b_g[edge_idx][b_idx] * bw_factor; \
877 for (size_t i_ = 0; i_ < curr_n; ++i_) curr_buffer[i_].weight *= factor; \
878 }
879#define MULT_COEFF(edge_idx) \
880 { \
881 double factor = 1.0 / (2.0 * bw[edge_idx] * b_g[edge_idx][edge_idx]); \
882 for (size_t i_ = 0; i_ < n; ++i_) \
883 weight_vector_data_buffer[i_].weight *= factor; \
884 }
885#define COMPACT_WEIGHTS(edge_idx) \
886 { \
887 compact_weight_vector_data( \
888 weight_vector_data_buffer, &n, weight_vector_data_pos_buffer); \
889 curr_triangle->c_111_a[edge_idx].data = \
890 xmalloc(n * sizeof(*(curr_triangle->c_111_a[edge_idx].data))); \
891 memcpy(curr_triangle->c_111_a[edge_idx].data, weight_vector_data_buffer, \
892 n * sizeof(*weight_vector_data_buffer)); \
893 curr_triangle->c_111_a[edge_idx].n = n; \
894 }
895
896 // a_0 = (D_g f(w0)/3 -
897 // B_020(w0) * (b_0(g0)*c_120 + b_1(g0)*c_030 + b_2(g0)*c_021) -
898 // B_011(w0) * ( b_1(g0)*c_021 + b_2(g0)*c_012) -
899 // B_002(w0) * (b_0(g0)*c_102 + b_1(g0)*c_012 + b_2(g0)*c_003)) /
900 // (B_011(w0) * b_0(g0))
901 size_t n = 0;
902 COPY_D_DATA(0); // D_g f(w0)/3
903 COPY_COEFF(c_120, 0, 0, 1.0); // - B_020(w0) * b_0(g0) * c_120
904 COPY_COEFF(c_030, 0, 1, 1.0); // - B_020(w0) * b_1(g0) * c_030
905 COPY_COEFF(c_021, 0, 2, 1.0); // - B_020(w0) * b_2(g0) * c_021
906 COPY_COEFF(c_021, 0, 1, 2.0); // - B_011(w0) * b_1(g0) * c_021
907 COPY_COEFF(c_012, 0, 2, 2.0); // - B_011(w0) * b_2(g0) * c_012
908 COPY_COEFF(c_102, 0, 0, 1.0); // - B_002(w0) * b_0(g0) * c_102
909 COPY_COEFF(c_012, 0, 1, 1.0); // - B_002(w0) * b_1(g0) * c_012
910 COPY_COEFF(c_003, 0, 2, 1.0); // - B_002(w0) * b_2(g0) * c_003
911 MULT_COEFF(0); // (...) / (B_011(w0) * b_0(g0))
912 COMPACT_WEIGHTS(0); // a_0 = ...
913
914 // a_1 = (D_g f(w1)/3 -
915 // B_200(w1) * (b_0(g1)*c_300 + b_1(g1)*c_210 + b_2(g1)*c_201) -
916 // B_101(w1) * (b_0(g1)*c_201 + b_2(g1)*c_102) -
917 // B_002(w1) * (b_0(g1)*c_102 + b_1(g1)*c_012 + b_2(g1)*c_003)) /
918 // (B_101(w1) * b_1(g1))
919 n = 0;
920 COPY_D_DATA(1); // D_g f(w1)/3
921 COPY_COEFF(c_300, 1, 0, 1.0); // - B_200(w1) * b_0(g1) * c_300
922 COPY_COEFF(c_210, 1, 1, 1.0); // - B_200(w1) * b_1(g1) * c_210
923 COPY_COEFF(c_201, 1, 2, 1.0); // - B_200(w1) * b_2(g1) * c_201
924 COPY_COEFF(c_201, 1, 0, 2.0); // - B_101(w1) * b_0(g1) * c_201
925 COPY_COEFF(c_102, 1, 2, 2.0); // - B_101(w1) * b_2(g1) * c_102
926 COPY_COEFF(c_102, 1, 0, 1.0); // - B_002(w1) * b_0(g1) * c_102
927 COPY_COEFF(c_012, 1, 1, 1.0); // - B_002(w1) * b_1(g1) * c_012
928 COPY_COEFF(c_003, 1, 2, 1.0); // - B_002(w1) * b_2(g1) * c_003
929 MULT_COEFF(1); // (...) / (B_101(w1) * b_1(g1))
930 COMPACT_WEIGHTS(1); // a_1 = ...
931
932 // a_2 = (D_g f(w2)/3 -
933 // B_200(w2) * (b_0(g2)*c_300 + b_1(g2)*c_210 + b_2(g2)*c_201) -
934 // B_110(w2) * (b_0(g2)*c_210 + b_1(g2)*c_120 ) -
935 // B_020(w2) * (b_0(g2)*c_120 + b_1(g2)*c_030 + b_2(g2)*c_021)) /
936 // (B_110(w2) * b_2(g2))
937 n = 0;
938 COPY_D_DATA(2); // D_g f(w2)/3
939 COPY_COEFF(c_300, 2, 0, 1.0); // - B_200(w2) * b_0(g2) * c_300
940 COPY_COEFF(c_210, 2, 1, 1.0); // - B_200(w2) * b_1(g2) * c_210
941 COPY_COEFF(c_201, 2, 2, 1.0); // - B_200(w2) * b_2(g2) * c_201
942 COPY_COEFF(c_210, 2, 0, 2.0); // - B_110(w2) * b_0(g2) * c_210
943 COPY_COEFF(c_120, 2, 1, 2.0); // - B_110(w2) * b_1(g2) * c_120
944 COPY_COEFF(c_120, 2, 0, 1.0); // - B_020(w2) * b_0(g2) * c_120
945 COPY_COEFF(c_030, 2, 1, 1.0); // - B_020(w2) * b_1(g2) * c_030
946 COPY_COEFF(c_021, 2, 2, 1.0); // - B_020(w2) * b_2(g2) * c_021
947 MULT_COEFF(2); // (...) / (B_110(w2) * b_2(g2))
948 COMPACT_WEIGHTS(2); // a_2 = ...
949
950#undef COMPACT_WEIGHTS
951#undef MULT_COEFF
952#undef COPY_COEFF
953#undef COPY_D_DATA
954 }
955 free(weight_vector_data_buffer);
956 free(weight_vector_data_pos_buffer);
957}
958
960 struct tgt_point_search_data * tgt_point_data) {
961
963 (tgt_point_data->type_flag == ON_VERTEX) ||
964 (tgt_point_data->type_flag == ON_EDGE) ||
965 (tgt_point_data->type_flag == ON_TRIANGLE),
966 "ERROR(get_max_num_weights): invalid type_flag")
967 switch(tgt_point_data->type_flag) {
968
969 case (ON_VERTEX):
970 return 1;
971 case (ON_EDGE):
972 return 2 + tgt_point_data->interp_data.edge->c[0].n +
973 tgt_point_data->interp_data.edge->c[1].n;
974 default:
975 case (ON_TRIANGLE): {
976 struct triangle_interp_data * triangle_data =
977 tgt_point_data->interp_data.triangle;
978 return 3 +
979 triangle_data->edges[0]->c[0].n + triangle_data->edges[0]->c[1].n +
980 triangle_data->edges[1]->c[0].n + triangle_data->edges[1]->c[1].n +
981 triangle_data->edges[2]->c[0].n + triangle_data->edges[2]->c[1].n +
982 triangle_data->c_111_a[0].n +
983 triangle_data->c_111_a[1].n +
984 triangle_data->c_111_a[2].n;
985 }
986 };
987}
988
990 struct tgt_point_search_data * tgt_point_data,
991 struct weight_vector_data * weights, size_t * n) {
992
993 weights[0].weight = 1.0;
994 weights[0].point_id = tgt_point_data->src_points[0];
995 *n = 1;
996}
997
999 struct tgt_point_search_data * tgt_point_data,
1000 struct weight_vector_data * weights, size_t * n) {
1001
1002 // spherical barycentric coordinates of the target point
1003 double * sb_coords = tgt_point_data->sb_coords;
1004
1005 // compute the spherical Bernstein base polynomials for the given vertices
1006 // (since the target point is on the edge, we only need a subset of the
1007 // polynomials, because the rest are zero anyway)
1008 // B_ijk(p) = d!/(i!*j!*k!) * b_0(p)^i * b_1(p)^j * b_2(p)^k
1009 double B_300 = sb_coords[0] * sb_coords[0] * sb_coords[0];
1010 double B_030 = sb_coords[1] * sb_coords[1] * sb_coords[1];
1011 double B_210 = 3.0 * sb_coords[0] * sb_coords[0] * sb_coords[1];
1012 double B_120 = 3.0 * sb_coords[0] * sb_coords[1] * sb_coords[1];
1013
1014 // c_300 = f(v1)
1015 // c_030 = f(v2)
1016 // c_210 and c_120
1017 struct weight_vector * c = &(tgt_point_data->interp_data.edge->c[0]);
1018
1019 // f(p) = SUM(c_ijk * B_ijk) // for 300, 030, 210, and 120
1020
1021 size_t * src_points = tgt_point_data->src_points;
1022
1023 // c_300 * B_300
1024 // weights[0].weight = B_300 + B_210;
1025 weights[0].weight = B_300;
1026 weights[0].point_id = src_points[0];
1027 // c_030 * B_030
1028 // weights[1].weight = B_030 + B_120;
1029 weights[1].weight = B_030;
1030 weights[1].point_id = src_points[1];
1031 // c_210 * B_210
1032 struct weight_vector_data * curr_weights = weights + 2;
1033 memcpy(curr_weights, c[0].data, c[0].n * sizeof(*curr_weights));
1034 for (size_t i = 0; i < c[0].n; ++i) curr_weights[i].weight *= B_210;
1035 // c_120 * B_120
1036 curr_weights += c[0].n;
1037 memcpy(curr_weights, c[1].data, c[1].n * sizeof(*curr_weights));
1038 for (size_t i = 0; i < c[1].n; ++i) curr_weights[i].weight *= B_120;
1039
1040 *n = 2 + c[0].n + c[1].n;
1041}
1042
1044 double * sb_coords, double * A) {
1045
1046 double b[3] = {sb_coords[1]*sb_coords[1] * sb_coords[2]*sb_coords[2],
1047 sb_coords[0]*sb_coords[0] * sb_coords[2]*sb_coords[2],
1048 sb_coords[0]*sb_coords[0] * sb_coords[1]*sb_coords[1]};
1049 double temp = 1.0 / (b[0] + b[1] + b[2]);
1050
1051 for (size_t i = 0; i < 3; ++i) {
1052
1053 if (fabs(b[i]) < 1e-9) A[i] = 0.0;
1054 else A[i] = b[i] * temp;
1055 }
1056}
1057
1059 struct tgt_point_search_data * tgt_point_data,
1060 struct weight_vector_data * weights, size_t * n) {
1061
1062 // spherical barycentric coordinates of the target point
1063 double * sb_coords = tgt_point_data->sb_coords;
1064
1065 struct edge_interp_data * curr_edges[3] = {
1066 tgt_point_data->interp_data.triangle->edges[0],
1067 tgt_point_data->interp_data.triangle->edges[1],
1068 tgt_point_data->interp_data.triangle->edges[2]};
1069 struct vertex_interp_data * vertex_d_data[3] = {
1070 curr_edges[2]->vertex_d_data[0],
1071 curr_edges[2]->vertex_d_data[1],
1072 curr_edges[1]->vertex_d_data[1]};
1073
1074 // compute the spherical Bernstein base polynomials for the given vertices
1075 // B_ijk(p) = d!/(i!*j!*k!) * b_0(p)^i * b_1(p)^j * b_2(p)^k
1076 double B_vertex[3] =
1077 {sb_coords[0] * sb_coords[0] * sb_coords[0], // B_300
1078 sb_coords[1] * sb_coords[1] * sb_coords[1], // B_030
1079 sb_coords[2] * sb_coords[2] * sb_coords[2]}; // B_003
1080 double B_edge[3][2] =
1081 {{3.0 * sb_coords[1] * sb_coords[1] * sb_coords[2], // B_021
1082 3.0 * sb_coords[1] * sb_coords[2] * sb_coords[2]}, // B_012
1083 {3.0 * sb_coords[0] * sb_coords[0] * sb_coords[2], // B_201
1084 3.0 * sb_coords[0] * sb_coords[2] * sb_coords[2]}, // B_102
1085 {3.0 * sb_coords[0] * sb_coords[0] * sb_coords[1], // B_210
1086 3.0 * sb_coords[0] * sb_coords[1] * sb_coords[1]}}; // B_120
1087 double B_111 = 6.0 * sb_coords[0] * sb_coords[1] * sb_coords[2];
1088
1089 // f(p) = SUM(c_ijk * B_ijk) // for all (i+j+k) = 3
1090
1091 // the three vertices coefficients (300, 030, and 003)
1092 for (size_t i = 0; i < 3; ++i) {
1093 weights[i].weight = B_vertex[i];
1094 weights[i].point_id = vertex_d_data[i]->src_point;
1095 }
1096
1097 size_t n_ = 3;
1098
1099 // the six edges coefficients
1100 struct weight_vector * c_ijk[3][2] = {
1101 {&(curr_edges[0]->c[0]), // c_021
1102 &(curr_edges[0]->c[1])}, // c_012
1103 {&(curr_edges[1]->c[0]), // c_201
1104 &(curr_edges[1]->c[1])}, // c_102
1105 {&(curr_edges[2]->c[0]), // c_210
1106 &(curr_edges[2]->c[1])} // c_120
1107 };
1108 for (size_t edge_idx = 0; edge_idx < 3; ++edge_idx) {
1109 for (size_t i = 0; i < 2; ++i) {
1110 struct weight_vector_data * curr_weights = weights + n_;
1111 size_t curr_n = c_ijk[edge_idx][i]->n;
1112 n_ += curr_n;
1113 memcpy(curr_weights, c_ijk[edge_idx][i]->data,
1114 c_ijk[edge_idx][i]->n * sizeof(*curr_weights));
1115 double B = B_edge[edge_idx][i];
1116 for (size_t j = 0; j < curr_n; ++j)
1117 curr_weights[j].weight *= B;
1118 }
1119 }
1120
1121 // the middle vertex
1122 // c_111 = SUM(c_111_a[i] * A[i]) // for i = 0..2
1123 // where a is the blending function
1124 double A[3];
1125 evaluate_blending_function(sb_coords, &(A[0]));
1126
1127 struct weight_vector * c_111_a =
1128 tgt_point_data->interp_data.triangle->c_111_a;
1129 for (size_t i = 0; i < 3; ++i) {
1130 struct weight_vector_data * curr_weights = weights + n_;
1131 size_t curr_n = c_111_a[i].n;
1132 n_ += curr_n;
1133 memcpy(curr_weights, c_111_a[i].data, curr_n * sizeof(*curr_weights));
1134 double B_111_x_A_i = B_111 * A[i];
1135 for (size_t j = 0; j < curr_n; ++j)
1136 curr_weights[j].weight *= B_111_x_A_i;
1137 }
1138
1139 *n = n_;
1140}
1141
1143 struct tgt_point_search_data * tgt_point_data,
1144 struct weight_vector_data * weights, size_t * n,
1145 struct weight_vector_data_pos * compact_buffer) {
1146
1147 YAC_ASSERT(
1148 (tgt_point_data->type_flag == ON_VERTEX) ||
1149 (tgt_point_data->type_flag == ON_EDGE) ||
1150 (tgt_point_data->type_flag == ON_TRIANGLE),
1151 "ERROR(compute_hcsbb_weights): invalid type_flag")
1152
1153 switch(tgt_point_data->type_flag) {
1154
1155 case (ON_VERTEX):
1156 compute_hcsbb_weights_vertex(tgt_point_data, weights, n);
1157 break;
1158 case (ON_EDGE):
1159 compute_hcsbb_weights_edge(tgt_point_data, weights, n);
1160 break;
1161 default:
1162 case (ON_TRIANGLE):
1163 compute_hcsbb_weights_triangle(tgt_point_data, weights, n);
1164 break;
1165 }
1166
1167 compact_weight_vector_data(weights, n, compact_buffer);
1168}
1169
1171 struct tgt_point_search_data * tgt_point_data, size_t num_tgt_points,
1172 struct edge_interp_data * edge_data, size_t num_edges,
1173 struct triangle_interp_data * triangle_data, size_t num_triangles,
1174 struct weight_vector_data ** weights, size_t ** num_weights_per_tgt,
1175 size_t * total_num_weights) {
1176
1177 //--------------------------------------------------------------------------
1178 // prepare patch data for the computation of the weights for each individual
1179 // target point
1180 //--------------------------------------------------------------------------
1181
1182 // compute edge coefficients
1183 compute_edge_coefficients(edge_data, num_edges);
1184
1185 // compute triangle coefficients
1186 compute_triangle_coefficients(triangle_data, num_triangles);
1187
1188 struct weight_vector_data * weight_vector_data = NULL;
1189 size_t weight_vector_data_array_size = 0;
1190 size_t total_num_weights_ = 0;
1191 size_t * num_weights_per_tgt_ =
1192 xmalloc(num_tgt_points * sizeof(*num_weights_per_tgt_));
1193 struct weight_vector_data_pos * compact_buffer =
1194 xmalloc((3 * 1 + // c_300, c_030, c_003
1195 3 * (1 + HCSBB_NUM_GAUSS_POINTS * HCSBB_GAUSS_NNN) + // c_012, c_021, c_102, ...
1197 6 * (1 + HCSBB_NUM_GAUSS_POINTS * HCSBB_GAUSS_NNN))) * // c_111_a[3]
1198 sizeof(*compact_buffer));
1199
1200 //--------------------------------------
1201 // compute weights for the target points
1202 //--------------------------------------
1203
1204 for (size_t i = 0; i < (num_tgt_points + 127) / 128; ++i) {
1205
1206 size_t curr_num_tgt_points = MIN(128, num_tgt_points - i * 128);
1207 struct tgt_point_search_data * curr_tgt_point_data =
1208 tgt_point_data + i * 128;
1209 size_t * curr_num_weights_per_tgt = num_weights_per_tgt_ + i * 128;
1210
1211 size_t curr_max_num_weights = 0;
1212 for (size_t j = 0; j < curr_num_tgt_points; ++j)
1213 curr_max_num_weights += get_max_num_weights(curr_tgt_point_data + j);
1214
1215 ENSURE_ARRAY_SIZE(weight_vector_data, weight_vector_data_array_size,
1216 total_num_weights_ + curr_max_num_weights);
1217
1218 for (size_t j = 0; j < curr_num_tgt_points; ++j) {
1220 curr_tgt_point_data + j, weight_vector_data + total_num_weights_,
1221 curr_num_weights_per_tgt + j, compact_buffer);
1222 total_num_weights_ += curr_num_weights_per_tgt[j];
1223 }
1224 }
1225
1226 free(compact_buffer);
1227
1228 *weights =
1229 xrealloc(
1230 weight_vector_data, total_num_weights_ * sizeof(*weight_vector_data));
1231 *num_weights_per_tgt = num_weights_per_tgt_;
1232 *total_num_weights = total_num_weights_;
1233}
1234
1236 const void * a, const void * b) {
1237
1238 struct tgt_point_search_data const * a_ = a;
1239 struct tgt_point_search_data const * b_ = b;
1240
1241 int ret;
1242 if ((ret = a_->is_masked - b_->is_masked)) return ret;
1243 if (a_->is_masked) return 0;
1244 if ((ret = (int)(a_->type_flag) - (int)(b_->type_flag))) return ret;
1245 for (int i = 0; i <= (int)a_->type_flag; ++i)
1246 if ((ret = (a_->src_points[i] > b_->src_points[i]) -
1247 (a_->src_points[i] < b_->src_points[i]))) return ret;
1248 return (a_->tgt_point > b_->tgt_point) - (a_->tgt_point < b_->tgt_point);
1249}
1250
1252 size_t * points, yac_const_coordinate_pointer coords,
1253 double bnd_triangle[][3]) {
1254
1255 double point_coords[HCSBB_D_NNN][3];
1256
1257 for (size_t i = 0; i < HCSBB_D_NNN; ++i)
1258 memcpy(point_coords[i], coords[points[i]], 3 * sizeof(double));
1259
1261 &(point_coords[0][0]), HCSBB_D_NNN, bnd_triangle, 16);
1262}
1263
1264static int compare_bnd_triangles(void const * a, void const * b) {
1265
1266 double const * tri_a = a, * tri_b = b;
1267 int ret = 0;
1268
1269 for (size_t i = 0; (!ret) && (i < 9); ++i) {
1270 double a_ = tri_a[i], b_ = tri_b[i];
1271 if (fabs(a_ - b_) > 1e-9) ret = (a_ > b_) - (a_ < b_);
1272 }
1273 return ret;
1274}
1275
1277 double gauss_vertices[HCSBB_NUM_GAUSS_POINTS][3],
1278 double gauss_sb_coords[HCSBB_NUM_GAUSS_POINTS][3],
1279 double bnd_triangle[3][3]) {
1280
1281#if HCSBB_GAUSS_ORDER == 3
1282 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1283 {1.0-0.333333333333333-0.333333333333333,0.333333333333333,0.333333333333333},
1284 {1.0-0.200000000000000-0.200000000000000,0.200000000000000,0.200000000000000},
1285 {1.0-0.600000000000000-0.200000000000000,0.600000000000000,0.200000000000000},
1286 {1.0-0.200000000000000-0.600000000000000,0.200000000000000,0.600000000000000}
1287 };
1288#elif HCSBB_GAUSS_ORDER == 4
1289 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1290 {1.0-0.091576213509771-0.091576213509771,0.091576213509771,0.091576213509771},
1291 {1.0-0.816847572980459-0.091576213509771,0.816847572980459,0.091576213509771},
1292 {1.0-0.091576213509771-0.816847572980459,0.091576213509771,0.816847572980459},
1293 {1.0-0.445948490915965-0.445948490915965,0.445948490915965,0.445948490915965},
1294 {1.0-0.108103018168070-0.445948490915965,0.108103018168070,0.445948490915965},
1295 {1.0-0.445948490915965-0.108103018168070,0.445948490915965,0.108103018168070}
1296 };
1297#elif HCSBB_GAUSS_ORDER == 5
1298 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1299 {1.0-0.333333333333333-0.333333333333333,0.333333333333333,0.333333333333333},
1300 {1.0-0.101286507323456-0.101286507323456,0.101286507323456,0.101286507323456},
1301 {1.0-0.797426985353087-0.101286507323456,0.797426985353087,0.101286507323456},
1302 {1.0-0.101286507323456-0.797426985353087,0.101286507323456,0.797426985353087},
1303 {1.0-0.470142064105115-0.470142064105115,0.470142064105115,0.470142064105115},
1304 {1.0-0.059715871789770-0.470142064105115,0.059715871789770,0.470142064105115},
1305 {1.0-0.470142064105115-0.059715871789770,0.470142064105115,0.059715871789770}
1306 };
1307#elif HCSBB_GAUSS_ORDER == 6
1308 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1309 {1.0-0.063089014491502-0.063089014491502,0.063089014491502,0.063089014491502},
1310 {1.0-0.873821971016996-0.063089014491502,0.873821971016996,0.063089014491502},
1311 {1.0-0.063089014491502-0.873821971016996,0.063089014491502,0.873821971016996},
1312 {1.0-0.249286745170910-0.249286745170910,0.249286745170910,0.249286745170910},
1313 {1.0-0.501426509658179-0.249286745170910,0.501426509658179,0.249286745170910},
1314 {1.0-0.249286745170910-0.501426509658179,0.249286745170910,0.501426509658179},
1315 {1.0-0.310352451033785-0.053145049844816,0.310352451033785,0.053145049844816},
1316 {1.0-0.053145049844816-0.310352451033785,0.053145049844816,0.310352451033785},
1317 {1.0-0.636502499121399-0.053145049844816,0.636502499121399,0.053145049844816},
1318 {1.0-0.053145049844816-0.636502499121399,0.053145049844816,0.636502499121399},
1319 {1.0-0.636502499121399-0.310352451033785,0.636502499121399,0.310352451033785},
1320 {1.0-0.310352451033785-0.636502499121399,0.310352451033785,0.636502499121399}
1321 };
1322#elif HCSBB_GAUSS_ORDER == 7
1323 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1324 {1.0-0.333333333333333-0.333333333333333,0.333333333333333,0.333333333333333},
1325 {1.0-0.260345966079038-0.260345966079038,0.260345966079038,0.260345966079038},
1326 {1.0-0.479308067841923-0.260345966079038,0.479308067841923,0.260345966079038},
1327 {1.0-0.260345966079038-0.479308067841923,0.260345966079038,0.479308067841923},
1328 {1.0-0.065130102902216-0.065130102902216,0.065130102902216,0.065130102902216},
1329 {1.0-0.869739794195568-0.065130102902216,0.869739794195568,0.065130102902216},
1330 {1.0-0.065130102902216-0.869739794195568,0.065130102902216,0.869739794195568},
1331 {1.0-0.312865496004874-0.048690315425316,0.312865496004874,0.048690315425316},
1332 {1.0-0.048690315425316-0.312865496004874,0.048690315425316,0.312865496004874},
1333 {1.0-0.638444188569809-0.048690315425316,0.638444188569809,0.048690315425316},
1334 {1.0-0.048690315425316-0.638444188569809,0.048690315425316,0.638444188569809},
1335 {1.0-0.638444188569809-0.312865496004874,0.638444188569809,0.312865496004874},
1336 {1.0-0.312865496004874-0.638444188569809,0.312865496004874,0.638444188569809}
1337 };
1338#elif HCSBB_GAUSS_ORDER == 8
1339 static double const abscissa[HCSBB_NUM_GAUSS_POINTS][3] = {
1340 {1.0-0.333333333333333-0.333333333333333,0.333333333333333,0.333333333333333},
1341 {1.0-0.081414823414554-0.459292588292723,0.081414823414554,0.459292588292723},
1342 {1.0-0.459292588292723-0.081414823414554,0.459292588292723,0.081414823414554},
1343 {1.0-0.459292588292723-0.459292588292723,0.459292588292723,0.459292588292723},
1344 {1.0-0.658861384496480-0.170569307751760,0.658861384496480,0.170569307751760},
1345 {1.0-0.170569307751760-0.658861384496480,0.170569307751760,0.658861384496480},
1346 {1.0-0.170569307751760-0.170569307751760,0.170569307751760,0.170569307751760},
1347 {1.0-0.898905543365938-0.050547228317031,0.898905543365938,0.050547228317031},
1348 {1.0-0.050547228317031-0.898905543365938,0.050547228317031,0.898905543365938},
1349 {1.0-0.050547228317031-0.050547228317031,0.050547228317031,0.050547228317031},
1350 {1.0-0.008394777409958-0.263112829634638,0.008394777409958,0.263112829634638},
1351 {1.0-0.263112829634638-0.008394777409958,0.263112829634638,0.008394777409958},
1352 {1.0-0.008394777409958-0.728492392955404,0.008394777409958,0.728492392955404},
1353 {1.0-0.728492392955404-0.008394777409958,0.728492392955404,0.008394777409958},
1354 {1.0-0.263112829634638-0.728492392955404,0.263112829634638,0.728492392955404},
1355 {1.0-0.728492392955404-0.263112829634638,0.728492392955404,0.263112829634638}
1356 };
1357#else
1358 static double const abscissa[1][3];
1359#error "interpolation_method_hcsbb: the specified gauss order is currently not supported."
1360#endif
1361
1362 for (size_t i = 0; i < HCSBB_NUM_GAUSS_POINTS; ++i) {
1363
1364 gauss_vertices[i][0] = bnd_triangle[0][0] * abscissa[i][0] +
1365 bnd_triangle[1][0] * abscissa[i][1] +
1366 bnd_triangle[2][0] * abscissa[i][2];
1367 gauss_vertices[i][1] = bnd_triangle[0][1] * abscissa[i][0] +
1368 bnd_triangle[1][1] * abscissa[i][1] +
1369 bnd_triangle[2][1] * abscissa[i][2];
1370 gauss_vertices[i][2] = bnd_triangle[0][2] * abscissa[i][0] +
1371 bnd_triangle[1][2] * abscissa[i][1] +
1372 bnd_triangle[2][2] * abscissa[i][2];
1373 double scale = 1.0 / sqrt(gauss_vertices[i][0] * gauss_vertices[i][0] +
1374 gauss_vertices[i][1] * gauss_vertices[i][1] +
1375 gauss_vertices[i][2] * gauss_vertices[i][2]);
1376 gauss_vertices[i][0] *= scale;
1377 gauss_vertices[i][1] *= scale;
1378 gauss_vertices[i][2] *= scale;
1379 gauss_sb_coords[i][0] = abscissa[i][0] * scale;
1380 gauss_sb_coords[i][1] = abscissa[i][1] * scale;
1381 gauss_sb_coords[i][2] = abscissa[i][2] * scale;
1382 }
1383}
1384
1385// computes the spherical Bernstein basis polynomials
1386// for degree = 3 the polynomials B_ijk are ordered in the following order:
1387// B_003, B_012, B_021, B_030, B_102, B_111, B_120, B_201, B_210, B_300
1390 double sbb_polynomials[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_GAUSS_POINTS]) {
1391
1392#define POW_0(x) (1.0)
1393#define POW_1(x) ((x))
1394#define POW_2(x) ((x)*(x))
1395#define POW_3(x) ((x)*(x)*(x))
1396#define FAC_0 (1.0)
1397#define FAC_1 (1.0)
1398#define FAC_2 (2.0)
1399#define FAC_3 (6.0)
1400
1401 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1402 sbb_polynomials[0][vertex_idx] =
1403 (FAC_3 / (FAC_0 * FAC_0 * FAC_3)) * (POW_0(sb_coords[vertex_idx][0]) *
1404 POW_0(sb_coords[vertex_idx][1]) *
1405 POW_3(sb_coords[vertex_idx][2]));
1406 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1407 sbb_polynomials[1][vertex_idx] =
1408 (FAC_3 / (FAC_0 * FAC_1 * FAC_2)) * (POW_0(sb_coords[vertex_idx][0]) *
1409 POW_1(sb_coords[vertex_idx][1]) *
1410 POW_2(sb_coords[vertex_idx][2]));
1411 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1412 sbb_polynomials[2][vertex_idx] =
1413 (FAC_3 / (FAC_0 * FAC_2 * FAC_1)) * (POW_0(sb_coords[vertex_idx][0]) *
1414 POW_2(sb_coords[vertex_idx][1]) *
1415 POW_1(sb_coords[vertex_idx][2]));
1416 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1417 sbb_polynomials[3][vertex_idx] =
1418 (FAC_3 / (FAC_0 * FAC_3 * FAC_0)) * (POW_0(sb_coords[vertex_idx][0]) *
1419 POW_3(sb_coords[vertex_idx][1]) *
1420 POW_0(sb_coords[vertex_idx][2]));
1421 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1422 sbb_polynomials[4][vertex_idx] =
1423 (FAC_3 / (FAC_1 * FAC_0 * FAC_2)) * (POW_1(sb_coords[vertex_idx][0]) *
1424 POW_0(sb_coords[vertex_idx][1]) *
1425 POW_2(sb_coords[vertex_idx][2]));
1426 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1427 sbb_polynomials[5][vertex_idx] =
1428 (FAC_3 / (FAC_1 * FAC_1 * FAC_1)) * (POW_1(sb_coords[vertex_idx][0]) *
1429 POW_1(sb_coords[vertex_idx][1]) *
1430 POW_1(sb_coords[vertex_idx][2]));
1431 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1432 sbb_polynomials[6][vertex_idx] =
1433 (FAC_3 / (FAC_1 * FAC_2 * FAC_0)) * (POW_1(sb_coords[vertex_idx][0]) *
1434 POW_2(sb_coords[vertex_idx][1]) *
1435 POW_0(sb_coords[vertex_idx][2]));
1436 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1437 sbb_polynomials[7][vertex_idx] =
1438 (FAC_3 / (FAC_2 * FAC_0 * FAC_1)) * (POW_2(sb_coords[vertex_idx][0]) *
1439 POW_0(sb_coords[vertex_idx][1]) *
1440 POW_1(sb_coords[vertex_idx][2]));
1441 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1442 sbb_polynomials[8][vertex_idx] =
1443 (FAC_3 / (FAC_2 * FAC_1 * FAC_0)) * (POW_2(sb_coords[vertex_idx][0]) *
1444 POW_1(sb_coords[vertex_idx][1]) *
1445 POW_0(sb_coords[vertex_idx][2]));
1446 for (size_t vertex_idx = 0; vertex_idx < HCSBB_NUM_GAUSS_POINTS; ++vertex_idx)
1447 sbb_polynomials[9][vertex_idx] =
1448 (FAC_3 / (FAC_3 * FAC_0 * FAC_0)) * (POW_3(sb_coords[vertex_idx][0]) *
1449 POW_0(sb_coords[vertex_idx][1]) *
1450 POW_0(sb_coords[vertex_idx][2]));
1451
1452#undef POW_0
1453#undef POW_1
1454#undef POW_2
1455#undef POW_3
1456#undef FAC_0
1457#undef FAC_1
1458#undef FAC_2
1459#undef FAC_3
1460}
1461
1463
1464// LAPACKE_dsytrf_work and LAPACKE_dsytri might not be available (see yac_lapack_interface.h).
1465#ifdef YAC_LAPACK_NO_DSYTR
1466 lapack_int ipiv[HCSBB_NUM_SBB_COEFFS + 1] = {0};
1467 double work[HCSBB_NUM_SBB_COEFFS*HCSBB_NUM_SBB_COEFFS] = {0};
1468
1469 YAC_ASSERT(
1470 !LAPACKE_dgetrf(
1471 LAPACK_COL_MAJOR, (lapack_int) HCSBB_NUM_SBB_COEFFS,
1472 (lapack_int) HCSBB_NUM_SBB_COEFFS, &A[0][0],
1473 (lapack_int) HCSBB_NUM_SBB_COEFFS, ipiv), "internal ERROR: dgetrf")
1474
1475 YAC_ASSERT(
1476 !LAPACKE_dgetri_work(
1477 LAPACK_COL_MAJOR, (lapack_int) HCSBB_NUM_SBB_COEFFS, &A[0][0],
1478 (lapack_int) HCSBB_NUM_SBB_COEFFS, ipiv, work,
1480 "internal ERROR: dgetri")
1481#else
1482 lapack_int ipiv[HCSBB_NUM_SBB_COEFFS] = {0};
1483 double work[HCSBB_NUM_SBB_COEFFS] = {0};
1484
1485 YAC_ASSERT(
1486 !LAPACKE_dsytrf_work(
1487 LAPACK_COL_MAJOR, 'L', (lapack_int) HCSBB_NUM_SBB_COEFFS , &A[0][0],
1488 (lapack_int) HCSBB_NUM_SBB_COEFFS, ipiv, work,
1489 (lapack_int) HCSBB_NUM_SBB_COEFFS), "internal ERROR: dsytrf")
1490
1491 YAC_ASSERT(
1492 !LAPACKE_dsytri_work(
1493 LAPACK_COL_MAJOR, 'L', (lapack_int) HCSBB_NUM_SBB_COEFFS , &A[0][0],
1494 (lapack_int) HCSBB_NUM_SBB_COEFFS, ipiv, work),
1495 "internal ERROR: dsytri")
1496
1497 for (size_t i = 0; i < HCSBB_NUM_SBB_COEFFS; ++i)
1498 for (size_t j = i + 1; j < HCSBB_NUM_SBB_COEFFS; ++j)
1499 A[j][i] = A[i][j];
1500#endif
1501}
1502
1504 struct gauss_nnn_patch * gauss_nnn_patch, double (*bnd_triangle)[3],
1505 double (*gauss_vertices)[3],
1506 double gauss_sbb_polynomials[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_GAUSS_POINTS],
1508
1509 // compute bounding triangle
1510 memcpy(&(gauss_nnn_patch->bnd_triangle[0][0]), bnd_triangle,
1511 3 * sizeof(*bnd_triangle));
1512
1513 double gauss_sb_coords[HCSBB_NUM_GAUSS_POINTS][3];
1514
1515 // compute gauss points
1516 generate_gauss_legendre_points(gauss_vertices, gauss_sb_coords, bnd_triangle);
1517
1518 double (*w_g)[HCSBB_NUM_GAUSS_POINTS] =
1519 xmalloc(HCSBB_NUM_SBB_COEFFS * sizeof(w_g[0]));
1520 gauss_nnn_patch->weights = (double*)w_g;
1521
1522 // compute w_g
1523 {
1524 // compute the spherical Bernstein basis polynomials for the
1525 // gauss points (P_g)
1526 compute_sbb_polynomials_gauss_3d(gauss_sb_coords, gauss_sbb_polynomials);
1527
1528 // P_g^T * P_g
1529 // (actually this is a triangular matrix -> computing halve would be enough)
1530 for (size_t i = 0; i < HCSBB_NUM_SBB_COEFFS; ++i) {
1531 for (size_t j = 0; j < HCSBB_NUM_SBB_COEFFS; ++j) {
1532 double accu = 0;
1533 for (size_t k = 0; k < HCSBB_NUM_GAUSS_POINTS; ++k) {
1534 accu += gauss_sbb_polynomials[i][k] * gauss_sbb_polynomials[j][k];
1535 }
1536 C[i][j] = accu;
1537 }
1538 }
1539
1540 // (P_p^T * P_p)^-1
1541 inverse(C);
1542
1543 // (P_p^T * P_p)^-1 * P_p^T
1544 for (size_t i = 0; i < HCSBB_NUM_SBB_COEFFS; ++i) {
1545 for (size_t j = 0; j < HCSBB_NUM_GAUSS_POINTS; ++j) {
1546 double accu = 0;
1547 for (size_t k = 0; k < HCSBB_NUM_SBB_COEFFS; ++k) {
1548 accu += C[i][k] * gauss_sbb_polynomials[k][j];
1549 }
1550 w_g[i][j] = accu;
1551 }
1552 }
1553 }
1554}
1555
1557 struct gauss_nnn_patch * gauss_nnn_patches, size_t num_gauss_nnn_patches) {
1558
1559 for (size_t i = 0; i < num_gauss_nnn_patches; ++i) {
1560
1561 free(gauss_nnn_patches[i].src_points);
1562 free(gauss_nnn_patches[i].weights);
1563 }
1564
1565 free(gauss_nnn_patches);
1566}
1567
1568// This routine generates a triangle around each coordinate. The generated
1569// triangle contains at around HCSBB_D_NNN source points, this ensures that the
1570// size of the triangle scales with the resolution of the source grid in the
1571// area around the coordinate.
1573 size_t count, yac_coordinate_pointer coords,
1574 struct yac_interp_grid * interp_grid, size_t * num_unique_triangles_,
1575 double (**unique_bnd_triangles_)[3][3], size_t ** coord_to_triangle) {
1576
1577 // for all points at which we need to estimate the derivative, do a nearest
1578 // neighbour search
1579 size_t * nnn_search_results =
1580 xmalloc((HCSBB_D_NNN + 1) * count * sizeof(*nnn_search_results));
1582 interp_grid, coords, count, HCSBB_D_NNN, nnn_search_results);
1583
1584 // sort the results for each search coordinate by the global ids of the
1585 // result ids (we do not actually care about the distance to the search
1586 // point (we also add a reorder index)
1587 const_yac_int_pointer src_global_ids =
1589 for (size_t i = 0, j = count - 1; i < count; ++i, --j) {
1590 // get global_ids for the current search results
1591 size_t * curr_nnn_search_results =
1592 nnn_search_results + j * HCSBB_D_NNN;
1593 yac_int nnn_search_result_global_ids[HCSBB_D_NNN];
1594 for (int k = 0; k < HCSBB_D_NNN; ++k)
1595 nnn_search_result_global_ids[k] =
1596 src_global_ids[curr_nnn_search_results[k]];
1597 // sort results of current point
1599 nnn_search_result_global_ids, HCSBB_D_NNN, curr_nnn_search_results);
1600 // make space for reorder index
1601 memmove(nnn_search_results + j * (HCSBB_D_NNN + 1),
1602 curr_nnn_search_results,
1603 HCSBB_D_NNN * sizeof(*nnn_search_results));
1604 // set reorder index
1605 nnn_search_results[j * (HCSBB_D_NNN + 1) + HCSBB_D_NNN] = j;
1606 }
1607
1608 // sort all results in order to be able to finde duplicated results
1609 qsort(nnn_search_results, count,
1610 (HCSBB_D_NNN + 1) * sizeof(*nnn_search_results),
1612
1613 size_t * coord_to_search_results =
1614 xmalloc(count * sizeof(*coord_to_search_results));
1615 size_t num_unique_search_results = 0;
1616 size_t dummy_search_results[HCSBB_D_NNN + 1];
1617 for (size_t i = 0; i < HCSBB_D_NNN + 1; ++i)
1618 dummy_search_results[i] = SIZE_MAX;
1619 for (size_t i = 0, * prev_search_result = dummy_search_results;
1620 i < count; ++i) {
1621
1622 size_t * curr_search_result = nnn_search_results + i * (HCSBB_D_NNN + 1);
1623 size_t coord_idx = curr_search_result[HCSBB_D_NNN];
1624 if (compare_HCSBB_D_NNN_size_t(prev_search_result, curr_search_result)) {
1625 prev_search_result = curr_search_result;
1626 ++num_unique_search_results;
1627 }
1628 coord_to_search_results[coord_idx] = num_unique_search_results - 1;
1629 }
1630
1631 // Remark: we have to call yac_interp_grid_get_src_field_coords after the call
1632 // to yac_interp_grid_do_nnn_search_src, because the source grid might have been
1633 // changed by the search
1634 yac_const_coordinate_pointer src_point_coords =
1636
1637 // generate bnd_triangles for all unique search results
1638 struct bnd_triangle_reorder * bnd_triangles_reorder =
1639 xmalloc(num_unique_search_results * sizeof(*bnd_triangles_reorder));
1640 for (size_t i = 0, j = 0, * prev_search_result = dummy_search_results;
1641 i < count; ++i) {
1642
1643 size_t * curr_search_result = nnn_search_results + i * (HCSBB_D_NNN + 1);
1644 if (compare_HCSBB_D_NNN_size_t(prev_search_result, curr_search_result)) {
1646 curr_search_result, src_point_coords, bnd_triangles_reorder[j].coords);
1647 bnd_triangles_reorder[j].reorder_idx = j;
1648 prev_search_result = curr_search_result;
1649 ++j;
1650 }
1651 }
1652 free(nnn_search_results);
1653
1654 // sort generated bounding triangles
1655 qsort(bnd_triangles_reorder, num_unique_search_results,
1656 sizeof(*bnd_triangles_reorder), compare_bnd_triangles);
1657
1658 // determine number of unique bnd_triangles
1659 size_t * search_result_to_bnd_triangle =
1660 xmalloc(num_unique_search_results * sizeof(*search_result_to_bnd_triangle));
1661 double (*unique_bnd_triangles)[3][3] = (double(*)[3][3])bnd_triangles_reorder;
1662 size_t num_unique_triangles = 0;
1663 double dummy_bnd_triangle[3][3] = {{-1,-1,-1},{-1,-1,-1},{-1,-1,-1}};
1664 double (*prev_bnd_triangle)[3] = dummy_bnd_triangle;
1665 for (size_t i = 0; i < num_unique_search_results; ++i) {
1666 double (*curr_bnd_triangle)[3] = bnd_triangles_reorder[i].coords;
1667 size_t reorder_idx = bnd_triangles_reorder[i].reorder_idx;
1668 if (compare_bnd_triangles(prev_bnd_triangle, curr_bnd_triangle)) {
1669 memmove(unique_bnd_triangles[num_unique_triangles], curr_bnd_triangle,
1670 sizeof(*unique_bnd_triangles));
1671 ++num_unique_triangles;
1672 prev_bnd_triangle = curr_bnd_triangle;
1673 }
1674 search_result_to_bnd_triangle[reorder_idx] = num_unique_triangles - 1;
1675 }
1676
1677 // convert coord_to_search_results to coord_to_triangle
1678 for (size_t i = 0; i < count; ++i)
1679 coord_to_search_results[i] =
1680 search_result_to_bnd_triangle[coord_to_search_results[i]];
1681 free(search_result_to_bnd_triangle);
1682
1683 *num_unique_triangles_ = num_unique_triangles;
1684 *unique_bnd_triangles_ =
1685 xrealloc(unique_bnd_triangles,
1686 num_unique_triangles * sizeof(*unique_bnd_triangles));
1687 *coord_to_triangle =
1688 xrealloc(coord_to_search_results, count * sizeof(*coord_to_search_results));
1689}
1690
1691// To estimate the derivate at a given point we first build triangle around the
1692// point. The size of this triangle scales with the source grid resolution.
1693// For each triangle we generate a number of supporting points, which are
1694// interpolated using a simple nearest neighbour approach. Using these
1695// supporting points we can interpolate the source field based on
1696// Bernstein basis polynomials. Based on these polynomials, we can generate the
1697// weights to compute derivatives of points.
1699 size_t num_vertices, struct vertex_interp_data * vertex_data,
1700 size_t num_edges, struct edge_interp_data * edge_data,
1701 struct yac_interp_grid * interp_grid, size_t * num_gauss_nnn_patches_,
1702 struct gauss_nnn_patch ** gauss_nnn_patches_) {
1703
1704 // coordinates at which we will need to estimate a derivative
1706 xmalloc((num_vertices + num_edges) * sizeof(*coords));
1707
1708 // get coordinates of all vertices and all edge middle points
1709 for (size_t i = 0; i < num_vertices; ++i)
1710 memcpy(coords[i], vertex_data[i].coordinate_xyz, sizeof(*coords));
1711 for (size_t i = 0; i < num_edges; ++i)
1712 memcpy(
1713 coords[i + num_vertices], edge_data[i].middle_d_data.coordinate_xyz,
1714 sizeof(*coords));
1715
1716 size_t num_gauss_nnn_patches;
1717 double (*bnd_triangles)[3][3];
1718 size_t * coord_to_gauss_nnn_patch;
1720 num_vertices + num_edges, coords, interp_grid,
1721 &num_gauss_nnn_patches, &bnd_triangles, &coord_to_gauss_nnn_patch);
1722 free(coords);
1723
1724 struct gauss_nnn_patch * gauss_nnn_patches =
1725 xmalloc(num_gauss_nnn_patches * sizeof(*gauss_nnn_patches));
1726 double (*gauss_points)[HCSBB_NUM_GAUSS_POINTS][3] =
1727 xmalloc(num_gauss_nnn_patches * sizeof(*gauss_points));
1728 double (*gauss_sbb_polynomials_buffer)[HCSBB_NUM_GAUSS_POINTS] =
1729 xmalloc(HCSBB_NUM_SBB_COEFFS * sizeof(*gauss_sbb_polynomials_buffer));
1730 double (*work_buffer)[HCSBB_NUM_SBB_COEFFS] =
1731 xmalloc(HCSBB_NUM_SBB_COEFFS * sizeof(*work_buffer));
1732
1733 *num_gauss_nnn_patches_ = num_gauss_nnn_patches;
1734 *gauss_nnn_patches_ = gauss_nnn_patches;
1735
1736 // initialise all gauss patches
1737 for (size_t i = 0; i < num_gauss_nnn_patches; ++i)
1739 gauss_nnn_patches + i, bnd_triangles[i],
1740 gauss_points[i], gauss_sbb_polynomials_buffer, work_buffer);
1741 free(work_buffer);
1742 free(gauss_sbb_polynomials_buffer);
1743 free(bnd_triangles);
1744
1745 for (size_t i = 0; i < num_vertices; ++i)
1746 vertex_data[i].gauss_nnn_patch =
1747 gauss_nnn_patches + coord_to_gauss_nnn_patch[i];
1748 for (size_t i = 0; i < num_edges; ++i)
1749 edge_data[i].middle_d_data.gauss_nnn_patch =
1750 gauss_nnn_patches + coord_to_gauss_nnn_patch[i + num_vertices];
1751 free(coord_to_gauss_nnn_patch);
1752
1753 // do an NNN search for all gauss integration point
1754 size_t * nnn_search_results =
1755 xmalloc(HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS * num_gauss_nnn_patches *
1756 sizeof(*nnn_search_results));
1758 interp_grid, (yac_coordinate_pointer)&(gauss_points[0][0][0]),
1759 HCSBB_NUM_GAUSS_POINTS * num_gauss_nnn_patches,
1760 HCSBB_GAUSS_NNN, nnn_search_results);
1761
1762 // generate the nnn weights for the interpolation of the gauss integration
1763 // points which are required for the derivative estimates
1764 double (*w_nnn_buffer)[HCSBB_NUM_GAUSS_POINTS] =
1765 xmalloc(HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS * sizeof(*w_nnn_buffer));
1766 size_t * src_point_buffer =
1767 xmalloc(
1768 HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS * sizeof(*src_point_buffer));
1769 yac_int * global_id_buffer =
1770 xmalloc(
1771 HCSBB_GAUSS_NNN * HCSBB_NUM_GAUSS_POINTS * sizeof(*global_id_buffer));
1772 {
1773 yac_const_coordinate_pointer src_point_coords =
1775 const_yac_int_pointer src_point_global_ids =
1777 for (size_t i = 0; i < num_gauss_nnn_patches; ++i)
1779 gauss_nnn_patches + i, gauss_points[i],
1780 nnn_search_results + i * HCSBB_NUM_GAUSS_POINTS * HCSBB_GAUSS_NNN,
1781 w_nnn_buffer, src_point_buffer, global_id_buffer, src_point_coords,
1782 src_point_global_ids);
1783 }
1784 free(nnn_search_results);
1785 free(global_id_buffer);
1786 free(src_point_buffer);
1787 free(w_nnn_buffer);
1788 free(gauss_points);
1789}
1790
1791static int check_polygon(
1792 double * check_coord, int num_vertices, size_t const * vertices,
1793 size_t vertex_start_idx, const_yac_int_pointer global_ids,
1794 yac_const_coordinate_pointer point_coords,
1795 size_t * triangle, double * sb_coords) {
1796
1797 // number of triangle the current source cell can be split into
1798 size_t num_triangles = (size_t)(num_vertices - 2);
1799 size_t triangle_indices[num_triangles][3];
1801 vertices, num_vertices, vertex_start_idx, &(triangle_indices[0]));
1802
1803 int ret = 0;
1804 for (size_t i = 0; (!ret) && (i < num_triangles); ++i) {
1805
1806 triangle[0] = triangle_indices[i][0];
1807 triangle[1] = triangle_indices[i][1];
1808 triangle[2] = triangle_indices[i][2];
1809 yac_int triangle_ids[3] =
1810 {global_ids[triangle[0]],
1811 global_ids[triangle[1]],
1812 global_ids[triangle[2]]};
1813
1814 // this is done to ensure bit-reproducibility
1815 yac_quicksort_index_yac_int_size_t(triangle_ids, 3, triangle);
1816
1817 double triangle_coords[3][3];
1818 for (int j = 0; j < 3; ++j)
1819 memcpy(
1820 triangle_coords[j], point_coords[triangle[j]], 3 * sizeof(double));
1821
1822 sb_coords[0] = check_coord[0];
1823 sb_coords[1] = check_coord[1];
1824 sb_coords[2] = check_coord[2];
1825 compute_sb_coords(sb_coords, 1, triangle_coords);
1826
1827 // if the point is in the current triangle
1828 ret = (sb_coords[0] > -SB_COORD_TOL) &&
1829 (sb_coords[1] > -SB_COORD_TOL) &&
1830 (sb_coords[2] > -SB_COORD_TOL);
1831 }
1832
1833 return ret;
1834}
1835
1837 double * tgt_coord, size_t src_cell,
1838 struct yac_const_basic_grid_data * src_grid_data,
1839 yac_const_coordinate_pointer src_point_coords,
1840 size_t * src_triangle_vertices, double * sb_coords) {
1841
1842 size_t const * src_vertices =
1843 src_grid_data->cell_to_vertex +
1844 src_grid_data->cell_to_vertex_offsets[src_cell];
1845 int num_src_vertices = src_grid_data->num_vertices_per_cell[src_cell];
1846 yac_int const * vertex_ids = src_grid_data->ids[YAC_LOC_CORNER];
1847
1848 // Find start point for triangulation of current source cell. This ensures
1849 // that the result is decomposition independent.
1850 size_t lowest_vertex_id_idx = 0;
1851 {
1852 yac_int lowest_vertex_id = vertex_ids[src_vertices[0]];
1853 for (int i = 1; i < num_src_vertices; ++i) {
1854 if (vertex_ids[src_vertices[i]] < lowest_vertex_id) {
1855 lowest_vertex_id = vertex_ids[src_vertices[i]];
1856 lowest_vertex_id_idx = i;
1857 }
1858 }
1859 }
1860
1861 return
1863 tgt_coord, num_src_vertices, src_vertices, lowest_vertex_id_idx,
1864 vertex_ids, src_point_coords, src_triangle_vertices, sb_coords);
1865}
1866
1868 double * tgt_coord, size_t src_cell,
1869 struct yac_const_basic_grid_data * src_grid_data,
1870 yac_const_coordinate_pointer src_point_coords,
1871 size_t * vertex_to_cell, size_t * vertex_to_cell_offsets,
1872 int * num_cells_per_vertex, size_t * src_triangle, double * sb_coords) {
1873
1874 size_t num_vertices = src_grid_data->num_vertices_per_cell[src_cell];
1875 size_t const * curr_vertices =
1876 src_grid_data->cell_to_vertex +
1877 src_grid_data->cell_to_vertex_offsets[src_cell];
1878 yac_int const * cell_ids = src_grid_data->ids[YAC_LOC_CELL];
1879
1880 int ret = 0;
1881
1882 // for all auxiliary grid cells
1883 for (size_t i = 0; (!ret) && (i < num_vertices); ++i) {
1884
1885 size_t curr_vertex = curr_vertices[i];
1886
1887 int curr_num_src_cells = num_cells_per_vertex[curr_vertex];
1888 if (curr_num_src_cells == 0) continue;
1889
1890 size_t * curr_src_cells =
1891 vertex_to_cell + vertex_to_cell_offsets[curr_vertex];
1892
1893 // the auxiliary grid cells are already constructed in a way that ensures
1894 // that they are decomposition independent => start index == 0
1895 ret =
1897 tgt_coord, curr_num_src_cells, curr_src_cells, 0, cell_ids,
1898 src_point_coords, src_triangle, sb_coords);
1899 }
1900
1901 return ret;
1902}
1903
1905 size_t num_tgt_points, size_t * selected_tgt, size_t * tgt_points,
1906 yac_coordinate_pointer tgt_coords, size_t * src_cells,
1907 struct yac_interp_grid * interp_grid) {
1908
1909 enum yac_location src_location =
1911
1912 struct tgt_point_search_data * tgt_point_data =
1913 xmalloc(num_tgt_points * sizeof(*tgt_point_data));
1914
1915 // if the source data is located at cells, we need to have information on the
1916 // auxiliary grid (cell points are the vertices of this grid)
1917 size_t * vertex_to_cell = NULL;
1918 size_t * vertex_to_cell_offsets = NULL;
1919 int * num_cells_per_vertex = NULL;
1920 if (src_location == YAC_LOC_CELL)
1922 interp_grid, src_cells, num_tgt_points,
1923 &vertex_to_cell, &vertex_to_cell_offsets, &num_cells_per_vertex);
1924
1925 // Remark: you have to get the pointers after the call to
1926 // yac_dist_grid_get_aux_grid_cells, because the routine might change
1927 // the underlying memory
1928 struct yac_const_basic_grid_data * src_grid_data =
1930 yac_const_coordinate_pointer src_point_coords =
1932 const_int_pointer src_mask =
1933 yac_interp_grid_get_src_field_mask(interp_grid, 0);
1934
1935 for (size_t i = 0; i < num_tgt_points; ++i, ++tgt_point_data) {
1936
1937 size_t tgt_idx = selected_tgt[i];
1938
1939 double * tgt_coord = tgt_coords[tgt_idx];
1940 tgt_point_data->tgt_point = tgt_points[tgt_idx];
1941 memcpy(
1942 tgt_point_data->coordinate_xyz, tgt_coord, 3 * sizeof(*tgt_coord));
1943
1944 // For the target point we need to find a triangle of source points connected
1945 // by great circles. The triangulation of the source cells needs to be
1946 // identical for all source cells on all processes
1947 size_t src_triangle_vertices[3];
1948 double sb_coords[3];
1949 int successfull;
1950 YAC_ASSERT(
1951 (src_location == YAC_LOC_CORNER) || (src_location == YAC_LOC_CELL),
1952 "ERROR(init_tgt_point_search_data): invalid source location")
1953 if (src_location == YAC_LOC_CORNER) {
1954 successfull =
1956 tgt_coord, src_cells[i], src_grid_data, src_point_coords,
1957 src_triangle_vertices, sb_coords);
1958 } else {
1959 successfull =
1961 tgt_coord, src_cells[i], src_grid_data, src_point_coords,
1962 vertex_to_cell, vertex_to_cell_offsets, num_cells_per_vertex,
1963 src_triangle_vertices, sb_coords);
1964 }
1965
1966 if (successfull) {
1967
1968 int tmp_type_flag = -1;
1969 int is_masked = 0;
1970 for (size_t k = 0; k < 3; ++k) {
1971 if (fabs(sb_coords[k]) > SB_COORD_TOL) {
1972 tmp_type_flag++;
1973 tgt_point_data->src_points[tmp_type_flag] = src_triangle_vertices[k];
1974 tgt_point_data->sb_coords[tmp_type_flag] = sb_coords[k];
1975 // check the mask for the source vertices
1976 if (src_mask != NULL) is_masked |= !src_mask[src_triangle_vertices[k]];
1977 }
1978 }
1979 tgt_point_data->type_flag = (enum interp_type_flag)tmp_type_flag;
1980 tgt_point_data->is_masked = is_masked;
1981 } else {
1982 tgt_point_data->is_masked = 1;
1983 }
1984 }
1985 free(num_cells_per_vertex);
1986 free(vertex_to_cell_offsets);
1987 free(vertex_to_cell);
1988
1989 return tgt_point_data - num_tgt_points;
1990}
1991
1993 size_t num_vertices, size_t * vertices,
1994 struct yac_interp_grid * interp_grid) {
1995
1996 yac_const_coordinate_pointer src_point_coords =
1998
1999 struct vertex_interp_data * vertex_data =
2000 xmalloc(num_vertices * sizeof(*vertex_data));
2001
2002 for (size_t i = 0; i < num_vertices; ++i) {
2003 size_t vertex_idx = vertices[i];
2004 memcpy(vertex_data[i].coordinate_xyz,
2005 src_point_coords[vertex_idx], sizeof(*src_point_coords));
2006 vertex_data[i].src_point = vertex_idx;
2007 }
2008
2009 return vertex_data;
2010}
2011
2013 size_t num_edges, size_t (*edge_to_vertex_d_data)[2],
2014 struct vertex_interp_data * vertex_data) {
2015
2016 struct edge_interp_data * edge_data =
2017 xmalloc(num_edges * sizeof(*edge_data));
2018
2019 for (size_t i = 0; i < num_edges; ++i) {
2020
2021 // link edge to derivative data of its two vertices
2022 edge_data[i].vertex_d_data[0] = vertex_data + edge_to_vertex_d_data[i][0];
2023 edge_data[i].vertex_d_data[1] = vertex_data + edge_to_vertex_d_data[i][1];
2024
2025 double * vertex_coordinates_xyz[2] = {
2026 edge_data[i].vertex_d_data[0]->coordinate_xyz,
2027 edge_data[i].vertex_d_data[1]->coordinate_xyz};
2028
2029 // compute middle point of edge
2030 double * middle_point = edge_data[i].middle_d_data.coordinate_xyz;
2031 middle_point[0] = vertex_coordinates_xyz[0][0] +
2032 vertex_coordinates_xyz[1][0];
2033 middle_point[1] = vertex_coordinates_xyz[0][1] +
2034 vertex_coordinates_xyz[1][1];
2035 middle_point[2] = vertex_coordinates_xyz[0][2] +
2036 vertex_coordinates_xyz[1][2];
2037 double sb_coord =
2038 1.0 / sqrt(middle_point[0] * middle_point[0] +
2039 middle_point[1] * middle_point[1] +
2040 middle_point[2] * middle_point[2]);
2041 middle_point[0] *= sb_coord;
2042 middle_point[1] *= sb_coord;
2043 middle_point[2] *= sb_coord;
2044 edge_data[i].sb_coord_middle_point = sb_coord;
2045
2046 // compute vector that is orthogonal to the edge and a tangent on the
2047 // sphere in the middle point
2048 crossproduct_d(vertex_coordinates_xyz[0], vertex_coordinates_xyz[1],
2049 &(edge_data[i].g[0]));
2050 }
2051
2052 return edge_data;
2053}
2054
2056 size_t num_triangles, size_t (*triangle_to_edge_data)[3],
2057 struct edge_interp_data * edge_data) {
2058
2059 struct triangle_interp_data * triangle_data =
2060 xmalloc(num_triangles * sizeof(*triangle_data));
2061
2062 for (size_t i = 0; i < num_triangles; ++i) {
2063
2064 // link triangle to the data of its edges
2065 triangle_data[i].edges[0] = edge_data + triangle_to_edge_data[i][0];
2066 triangle_data[i].edges[1] = edge_data + triangle_to_edge_data[i][1];
2067 triangle_data[i].edges[2] = edge_data + triangle_to_edge_data[i][2];
2068 }
2069
2070 return triangle_data;
2071}
2072
2073static size_t do_search_hcsbb (struct interp_method * method,
2074 struct yac_interp_grid * interp_grid,
2075 size_t * tgt_points, size_t count,
2076 struct yac_interp_weights * weights) {
2077
2078 UNUSED(method);
2079
2080 YAC_ASSERT(
2081 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
2082 "ERROR(do_search_hcsbb): invalid number of source fields")
2083
2084 YAC_ASSERT(
2085 (yac_interp_grid_get_src_field_location(interp_grid, 0) ==
2086 YAC_LOC_CORNER) ||
2087 (yac_interp_grid_get_src_field_location(interp_grid, 0) ==
2088 YAC_LOC_CELL),
2089 "ERROR(do_search_hcsbb): unsupported source field location type")
2090
2091 // get coordinates of target points
2092 yac_coordinate_pointer tgt_coords = xmalloc(count * sizeof(*tgt_coords));
2094 interp_grid, tgt_points, count, tgt_coords);
2095
2096 size_t * size_t_buffer = xmalloc(2 * count * sizeof(*size_t_buffer));
2097 size_t * src_cells = size_t_buffer;
2098 size_t * reorder_idx = size_t_buffer + count;
2099
2100 // get matching source cells for all target points
2101 // (this search assumes that all edges of the grid are great circles)
2103 interp_grid, tgt_coords, count, src_cells);
2104
2105 // sort target points, for which we found a source cell, to the beginning of
2106 // the array
2107 for (size_t i = 0; i < count; ++i) reorder_idx[i] = i;
2108 yac_quicksort_index_size_t_size_t(src_cells, count, reorder_idx);
2109
2110 // count the number of target for which a valid source cell was found
2111 size_t temp_result_count = 0;
2112 for (temp_result_count = 0; temp_result_count < count; ++temp_result_count)
2113 if (src_cells[temp_result_count] == SIZE_MAX) break;
2114
2115 // initialise search data for all target point, for which a matching source
2116 // cell was found
2117 struct tgt_point_search_data * tgt_point_data =
2119 temp_result_count, reorder_idx, tgt_points, tgt_coords, src_cells,
2120 interp_grid);
2121 free(tgt_coords);
2122
2123 // sort target point data (first by whether they are marked or not, second
2124 // by interp_type, and third by vertex indices)
2125 qsort(tgt_point_data, temp_result_count, sizeof(*tgt_point_data),
2127
2128 // move target points which cannot be interpolated to the end of the list
2129 for (size_t i = temp_result_count; i < count; ++i)
2130 reorder_idx[i] = tgt_points[reorder_idx[i]];
2131 memcpy(tgt_points + temp_result_count, reorder_idx + temp_result_count,
2132 (count - temp_result_count) * sizeof(*tgt_points));
2133 // put remaining target point to the front of the list (can also contain
2134 // failed target points)
2135 for (size_t i = 0; i < temp_result_count; ++i)
2136 tgt_points[i] = tgt_point_data[i].tgt_point;
2137 free(size_t_buffer);
2138
2139 // count the number of target points that can be interpolated
2140 size_t result_count = 0;
2141 size_t result_counts[3] = {0,0,0};
2142 for (result_count = 0; result_count < temp_result_count; ++result_count) {
2143 if (tgt_point_data[result_count].is_masked) break;
2144 else result_counts[(int)(tgt_point_data[result_count].type_flag)]++;
2145 }
2146
2147 // get all edges and triangles required for the interpolation
2148 // (each triangle requires three edges, usually each triangle edges is shared
2149 // between two triangles)
2150 size_t * unique_vertices;
2151 size_t (*unique_edge_to_unique_vertex)[2];
2152 size_t (*unique_triangle_to_unique_edge)[3];
2153 size_t num_unique_vertices, num_unique_edges, num_unique_triangles;
2155 tgt_point_data + result_counts[0], result_counts[1], result_counts[2],
2156 &unique_vertices, &num_unique_vertices,
2157 &unique_edge_to_unique_vertex, &num_unique_edges,
2158 &unique_triangle_to_unique_edge, &num_unique_triangles);
2159
2160 // initialise data for all unique vertices, edges and triangles
2161 struct vertex_interp_data * vertex_data =
2163 num_unique_vertices, unique_vertices, interp_grid);
2164 struct edge_interp_data * edge_data =
2166 num_unique_edges, unique_edge_to_unique_vertex, vertex_data);
2167 struct triangle_interp_data * triangle_data =
2169 num_unique_triangles, unique_triangle_to_unique_edge, edge_data);
2170 free(unique_vertices);
2171 free(unique_edge_to_unique_vertex);
2172 free(unique_triangle_to_unique_edge);
2173
2174 // link tgt_point_data to the edge_data and triangle_data
2175 for (size_t i = result_counts[0];
2176 i < result_counts[0] + result_counts[1]; ++i)
2177 tgt_point_data[i].interp_data.edge =
2178 edge_data + tgt_point_data[i].interp_data.idx;
2179 for (size_t i = result_counts[0] + result_counts[1];
2180 i < result_counts[0] + result_counts[1] + result_counts[2]; ++i)
2181 tgt_point_data[i].interp_data.triangle =
2182 triangle_data + tgt_point_data[i].interp_data.idx;
2183
2184 // for the interpolation we will need to estimate derivatives at all required
2185 // vertices and the middle points of all edges
2186 size_t num_gauss_nnn_patches;
2187 struct gauss_nnn_patch * gauss_nnn_patches; // have to free this later
2189 num_unique_vertices, vertex_data, num_unique_edges, edge_data, interp_grid,
2190 &num_gauss_nnn_patches, &gauss_nnn_patches);
2191
2192 // compute the actual weights
2194 size_t * num_weights_per_tgt, total_num_weights;
2195 compute_weights(tgt_point_data, result_count,
2196 edge_data, num_unique_edges,
2197 triangle_data, num_unique_triangles,
2198 &weight_vector, &num_weights_per_tgt, &total_num_weights);
2199 free(vertex_data);
2200 if (num_unique_edges > 0) free(edge_data->c[0].data);
2201 free(edge_data);
2202 for (size_t i = 0; i < num_unique_triangles; ++i)
2203 for (size_t j = 0; j < 3; ++j) free(triangle_data[i].c_111_a[j].data);
2204 free(triangle_data);
2205 free(tgt_point_data);
2206
2207 // free memory associated with the gauss_nnn_patches
2208 free_gauss_nnn_patches(gauss_nnn_patches, num_gauss_nnn_patches);
2209
2210 size_t * src_points = xmalloc(total_num_weights * sizeof(*src_points));
2211 double * weight_data = xmalloc(total_num_weights * sizeof(*weight_data));
2212 for (size_t i = 0; i < total_num_weights; ++i) {
2213 weight_data[i] = weight_vector[i].weight;
2214 src_points[i] = weight_vector[i].point_id;
2215 }
2216 free(weight_vector);
2217
2218 struct remote_point * srcs =
2220 interp_grid, 0, src_points, total_num_weights);
2221
2222 struct remote_points tgts = {
2223 .data =
2225 interp_grid, tgt_points, result_count),
2226 .count = result_count};
2227
2228 // store weights
2230 weights, &tgts, num_weights_per_tgt, srcs, weight_data);
2231
2232 free(tgts.data);
2233 free(srcs);
2234 free(weight_data);
2235 free(src_points);
2236 free(num_weights_per_tgt);
2237
2238 return result_count;
2239}
2240
2242
2243 struct interp_method_hcsbb * method = xmalloc(1 * sizeof(*method));
2244
2246
2247 return (struct interp_method*)method;
2248}
2249
2250static void delete_hcsbb(struct interp_method * method) {
2251 free(method);
2252}
void yac_compute_bnd_triangle(double *vertices, size_t num_vertices, double triangle[][3], size_t num_tests)
#define UNUSED(x)
Definition core.h:71
int const * const_int_pointer
yac_int const * const_yac_int_pointer
#define ENSURE_ARRAY_SIZE(arrayp, curr_array_size, req_size)
void yac_triangulate_cell_indices(size_t const *corner_indices, size_t num_corners, size_t start_corner, size_t(*triangle_indices)[3])
Definition grid_cell.c:144
static void crossproduct_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
#define yac_angle_tol
Definition geometry.h:34
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
size_t yac_interp_grid_get_num_src_fields(struct yac_interp_grid *interp_grid)
const_int_pointer yac_interp_grid_get_src_field_mask(struct yac_interp_grid *interp_grid, size_t src_field_idx)
void yac_interp_grid_get_aux_grid_src(struct yac_interp_grid *interp_grid, size_t *cells, size_t count, size_t **vertex_to_cell, size_t **vertex_to_cell_offsets, int **num_cells_per_vertex)
const_yac_int_pointer yac_interp_grid_get_src_field_global_ids(struct yac_interp_grid *interp_grid, size_t src_field_idx)
struct remote_point * yac_interp_grid_get_tgt_remote_points(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count)
enum yac_location yac_interp_grid_get_src_field_location(struct yac_interp_grid *interp_grid, size_t src_field_idx)
struct remote_point * yac_interp_grid_get_src_remote_points(struct yac_interp_grid *interp_grid, size_t src_field_idx, size_t *src_points, size_t count)
yac_const_coordinate_pointer yac_interp_grid_get_src_field_coords(struct yac_interp_grid *interp_grid, size_t src_field_idx)
void yac_interp_grid_do_nnn_search_src(struct yac_interp_grid *interp_grid, yac_coordinate_pointer search_coords, size_t count, size_t n, size_t *src_points)
void yac_interp_grid_get_tgt_coordinates(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, yac_coordinate_pointer tgt_coordinates)
void yac_interp_grid_do_points_search_gc(struct yac_interp_grid *interp_grid, yac_coordinate_pointer search_coords, size_t count, size_t *src_cells)
struct yac_const_basic_grid_data * yac_interp_grid_get_basic_grid_data_src(struct yac_interp_grid *interp_grid)
static struct triangle_interp_data * init_triangle_interp_data(size_t num_triangles, size_t(*triangle_to_edge_data)[3], struct edge_interp_data *edge_data)
static struct edge_interp_data * init_edge_interp_data(size_t num_edges, size_t(*edge_to_vertex_d_data)[2], struct vertex_interp_data *vertex_data)
static int compare_tgt_point_search_data(const void *a, const void *b)
static void compute_weights(struct tgt_point_search_data *tgt_point_data, size_t num_tgt_points, struct edge_interp_data *edge_data, size_t num_edges, struct triangle_interp_data *triangle_data, size_t num_triangles, struct weight_vector_data **weights, size_t **num_weights_per_tgt, size_t *total_num_weights)
#define FAC_0
static void compute_gauss_nnn_patch(struct gauss_nnn_patch *gauss_nnn_patch, double(*gauss_points)[3], size_t *nnn_search_results, double w_nnn[][HCSBB_NUM_GAUSS_POINTS], size_t *src_point_buffer, yac_int *global_id_buffer, yac_const_coordinate_pointer src_point_coords, const_yac_int_pointer src_point_global_ids)
static struct vertex_interp_data * init_vertex_interp_data(size_t num_vertices, size_t *vertices, struct yac_interp_grid *interp_grid)
static int compare_HCSBB_D_NNN_size_t(const void *a, const void *b)
static int compare_size_t(const void *a, const void *b)
static void generate_d_data_triangle(size_t count, yac_coordinate_pointer coords, struct yac_interp_grid *interp_grid, size_t *num_unique_triangles_, double(**unique_bnd_triangles_)[3][3], size_t **coord_to_triangle)
#define HCSBB_D_NNN
static struct interp_method_vtable interp_method_hcsbb_vtable
#define HCSBB_GAUSS_NNN
#define HCSBB_NUM_SBB_COEFFS
#define HCSBB_NUM_GAUSS_POINTS
@ ON_EDGE
the target point is on an edge of the source point grid
@ ON_TRIANGLE
static void compute_triangle_coefficients(struct triangle_interp_data *triangle_data, size_t num_triangles)
#define COMPACT_WEIGHTS(edge_idx)
#define COPY_D_DATA(edge_idx)
static void generate_derivative_data(size_t num_vertices, struct vertex_interp_data *vertex_data, size_t num_edges, struct edge_interp_data *edge_data, struct yac_interp_grid *interp_grid, size_t *num_gauss_nnn_patches_, struct gauss_nnn_patch **gauss_nnn_patches_)
static void compute_hcsbb_weights_vertex(struct tgt_point_search_data *tgt_point_data, struct weight_vector_data *weights, size_t *n)
#define MULT_COEFF(edge_idx)
static void compute_sb_coords(double *sb_coords, size_t num_vertices, double triangle[3][3])
static void generate_bnd_triangle(size_t *points, yac_const_coordinate_pointer coords, double bnd_triangle[][3])
static void compute_hcsbb_weights_edge(struct tgt_point_search_data *tgt_point_data, struct weight_vector_data *weights, size_t *n)
static void inverse(double A[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_SBB_COEFFS])
static void compute_hcsbb_weights(struct tgt_point_search_data *tgt_point_data, struct weight_vector_data *weights, size_t *n, struct weight_vector_data_pos *compact_buffer)
static void compact_weight_vector_data(struct weight_vector_data *weights, size_t *n, struct weight_vector_data_pos *buffer)
static void get_unique_interp_data(struct tgt_point_search_data *tgt_point_data, size_t num_edges, size_t num_triangles, size_t **unique_vertices, size_t *num_unique_vertices, size_t(**unique_edge_to_unique_vertex)[2], size_t *num_unique_edges, size_t(**unique_triangle_to_unique_edge)[3], size_t *num_unique_triangles)
static void delete_hcsbb(struct interp_method *method)
static void compute_edge_coefficients(struct edge_interp_data *edge_data, size_t num_edges)
static int compare_bnd_triangles(void const *a, void const *b)
static void compute_sbb_polynomials_gauss_3d(double sb_coords[HCSBB_NUM_GAUSS_POINTS][3], double sbb_polynomials[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_GAUSS_POINTS])
#define COPY_COEFF(c_ijk, edge_idx, b_idx, bw_factor)
#define FAC_2
static int compare_2_size_t(const void *a, const void *b)
static void init_gauss_nnn_patch(struct gauss_nnn_patch *gauss_nnn_patch, double(*bnd_triangle)[3], double(*gauss_vertices)[3], double gauss_sbb_polynomials[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_GAUSS_POINTS], double C[HCSBB_NUM_SBB_COEFFS][HCSBB_NUM_SBB_COEFFS])
static int check_polygon(double *check_coord, int num_vertices, size_t const *vertices, size_t vertex_start_idx, const_yac_int_pointer global_ids, yac_const_coordinate_pointer point_coords, size_t *triangle, double *sb_coords)
static int compare_3_size_t(const void *a, const void *b)
struct interp_method * yac_interp_method_hcsbb_new()
static int compare_weight_vector_data_point(void const *a, void const *b)
static void compute_hcsbb_weights_triangle(struct tgt_point_search_data *tgt_point_data, struct weight_vector_data *weights, size_t *n)
static size_t do_search_hcsbb(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static void generate_gauss_legendre_points(double gauss_vertices[HCSBB_NUM_GAUSS_POINTS][3], double gauss_sb_coords[HCSBB_NUM_GAUSS_POINTS][3], double bnd_triangle[3][3])
#define FAC_3
static int compare_weight_vector_data_pos_pos(void const *a, void const *b)
#define POW_2(x)
#define POW_1(x)
static size_t get_max_num_weights(struct tgt_point_search_data *tgt_point_data)
static void evaluate_blending_function(double *sb_coords, double *A)
static int compare_weight_vector_data_pos_weight(void const *a, void const *b)
static int get_matching_aux_cell_triangle(double *tgt_coord, size_t src_cell, struct yac_const_basic_grid_data *src_grid_data, yac_const_coordinate_pointer src_point_coords, size_t *vertex_to_cell, size_t *vertex_to_cell_offsets, int *num_cells_per_vertex, size_t *src_triangle, double *sb_coords)
#define FAC_1
#define POW_0(x)
static void free_gauss_nnn_patches(struct gauss_nnn_patch *gauss_nnn_patches, size_t num_gauss_nnn_patches)
#define SB_COORD_TOL
static int get_matching_vertex_triangle(double *tgt_coord, size_t src_cell, struct yac_const_basic_grid_data *src_grid_data, yac_const_coordinate_pointer src_point_coords, size_t *src_triangle_vertices, double *sb_coords)
#define POW_3(x)
static struct tgt_point_search_data * init_tgt_point_data(size_t num_tgt_points, size_t *selected_tgt, size_t *tgt_points, yac_coordinate_pointer tgt_coords, size_t *src_cells, struct yac_interp_grid *interp_grid)
static void get_derivative_weights(struct gauss_nnn_patch *gauss_nnn_patch, double coordinate_xyz[3], double direction[3], struct weight_vector_data *weights, double mult_factor)
static void compute_d_sbb_polynomials_3d(double bnd_triangle[3][3], double direction[3], double coordinate_xyz[3], double d_sbb_polynomials[HCSBB_NUM_SBB_COEFFS])
void yac_interp_weights_add_wsum(struct yac_interp_weights *weights, struct remote_points *tgts, size_t *num_src_per_tgt, struct remote_point *srcs, double *w)
yac_location
Definition location.h:12
@ YAC_LOC_CORNER
Definition location.h:15
@ YAC_LOC_CELL
Definition location.h:14
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xmalloc(size)
Definition ppm_xfuncs.h:66
double coordinate_xyz[3]
point at which a derivative is to be estimated
struct edge_interp_data::@20 middle_d_data
derivative data of edge middle point
struct weight_vector c[2]
struct vertex_interp_data * vertex_d_data[2]
derivative data of both vertices
struct gauss_nnn_patch * gauss_nnn_patch
gauss integration point patch used to estimate the derivative
double bnd_triangle[3][3]
bounding triangle
struct interp_method_vtable * vtable
size_t(* do_search)(struct interp_method *method, struct yac_interp_grid *grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
struct remote_point * data
enum interp_type_flag type_flag
depending on type flag interp_data contains a point to the required data
struct triangle_interp_data * triangle
union tgt_point_search_data::@21 interp_data
struct edge_interp_data * edge
struct weight_vector c_111_a[3]
weights for alpha values used to compute c_111
struct edge_interp_data * edges[3]
edges of the triangle (edge[0] is the on opposing corner src_points[0])
size_t src_point
index of the associacted source point
double coordinate_xyz[3]
coordinate of the source point
struct gauss_nnn_patch * gauss_nnn_patch
struct weight_vector_data * weight
struct weight_vector_data * data
const_size_t_pointer cell_to_vertex_offsets
const_size_t_pointer cell_to_vertex
const const_int_pointer num_vertices_per_cell
const const_yac_int_pointer ids[3]
#define MIN(a, b)
void yac_quicksort_index_yac_int_size_t(yac_int *a, size_t n, size_t *idx)
static void yac_remove_duplicates_size_t_3(size_t(*array)[3], size_t *n)
Definition utils_core.h:169
static void yac_remove_duplicates_size_t(size_t *array, size_t *n)
Definition utils_core.h:108
void yac_quicksort_index_size_t_size_t(size_t *a, size_t n, size_t *idx)
static void yac_remove_duplicates_size_t_2(size_t(*array)[2], size_t *n)
Definition utils_core.h:135
static struct user_input_data_points ** points
Definition yac.c:120
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16
Xt_int yac_int
Definition yac_types.h:15
double const (* yac_const_coordinate_pointer)[3]
Definition yac_types.h:20
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19