YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
bnd_circle.c
Go to the documentation of this file.
1// Copyright (c) 2024 The YAC Authors
2//
3// SPDX-License-Identifier: BSD-3-Clause
4
5#ifdef HAVE_CONFIG_H
6// Get the definition of the 'restrict' keyword.
7#include "config.h"
8#endif
9
10#include <math.h>
11#include <stdlib.h>
12#include <string.h>
13#include <float.h>
14
15#include "geometry.h"
16#include "utils_core.h"
17#include "basic_grid.h"
18#include "ensure_array_size.h"
19
25// computes the circumscribe circle of a quad on the sphere
26// it is assumed that the edges are circles of longitude and latitude
28 double a[3], double b[3], double c[3],
29 struct bounding_circle * bnd_circle) {
30
32}
33
34// computes the bounding circle of a quad on the sphere
35// it is assumed that the edges are circles of longitude and latitude
37 double a[3], double b[3], double c[3],
38 struct bounding_circle * bnd_circle) {
39
41}
42
43// computes the circumscribe circle of a triangle on the sphere
44// it is assumed that all edges are great circles
46 double a[3], double b[3], double c[3],
47 struct bounding_circle * bnd_circle) {
48
49 double ab[3] = {a[0]-b[0], a[1]-b[1], a[2]-b[2]},
50 ac[3] = {b[0]-c[0], b[1]-c[1], b[2]-c[2]};
51
52 // it is assumed that the angles of a triangle do not get too small...
53 // crossproduct_kahan(ab, ac, bnd_circle->base_vector);
54 crossproduct_d(ab, ac, bnd_circle->base_vector);
55 normalise_vector(bnd_circle->base_vector);
56
57 // find biggest component of base_vector
58 double fabs_base_vector[3] = {fabs(bnd_circle->base_vector[0]),
59 fabs(bnd_circle->base_vector[1]),
60 fabs(bnd_circle->base_vector[2])};
61 int biggest_component_index =
62 ((fabs_base_vector[0] < fabs_base_vector[1]) |
63 (fabs_base_vector[0] < fabs_base_vector[2])) <<
64 (fabs_base_vector[1] < fabs_base_vector[2]);
65
66 if ((bnd_circle->base_vector[biggest_component_index] > 0) ^
67 (a[biggest_component_index] > 0)) {
68 bnd_circle->base_vector[0] = -bnd_circle->base_vector[0];
69 bnd_circle->base_vector[1] = -bnd_circle->base_vector[1];
70 bnd_circle->base_vector[2] = -bnd_circle->base_vector[2];
71 }
72
73 bnd_circle->inc_angle =
76 bnd_circle->base_vector, a), SIN_COS_TOL);
77 bnd_circle->sq_crd = DBL_MAX;
78}
79
80static inline double get_sin_vector_angle(
81 double a[3], double b[3]) {
82
83 double cross_ab[3];
84 crossproduct_kahan(a, b, cross_ab);
85
86 return sqrt(cross_ab[0]*cross_ab[0] +
87 cross_ab[1]*cross_ab[1] +
88 cross_ab[2]*cross_ab[2]);
89}
90
91// computes the bounding circle of a triangle on the sphere
92// it is assumed that all edges are great circles
94 double a[3], double b[3], double c[3],
95 struct bounding_circle * bnd_circle) {
96
97 double middle_point[3];
98
99 middle_point[0] = a[0] + b[0] + c[0];
100 middle_point[1] = a[1] + b[1] + c[1];
101 middle_point[2] = a[2] + b[2] + c[2];
102
103 normalise_vector(middle_point);
104
105 double cos_angles[3] = {middle_point[0] * a[0] +
106 middle_point[1] * a[1] +
107 middle_point[2] * a[2],
108 middle_point[0] * b[0] +
109 middle_point[1] * b[1] +
110 middle_point[2] * b[2],
111 middle_point[0] * c[0] +
112 middle_point[1] * c[1] +
113 middle_point[2] * c[2]};
114
115 struct sin_cos_angle inc_angle;
116
117 // find the biggest angle
118
119 if (cos_angles[0] < cos_angles[1]) {
120 if (cos_angles[0] < cos_angles[2]) {
121 inc_angle =
122 sin_cos_angle_new(get_sin_vector_angle(middle_point, a), cos_angles[0]);
123 } else {
124 inc_angle =
125 sin_cos_angle_new(get_sin_vector_angle(middle_point, c), cos_angles[2]);
126 }
127 } else {
128 if (cos_angles[1] < cos_angles[2]) {
129 inc_angle =
130 sin_cos_angle_new(get_sin_vector_angle(middle_point, b), cos_angles[1]);
131 } else {
132 inc_angle =
133 sin_cos_angle_new(get_sin_vector_angle(middle_point, c), cos_angles[2]);
134 }
135 }
136
137 bnd_circle->base_vector[0] = middle_point[0];
138 bnd_circle->base_vector[1] = middle_point[1];
139 bnd_circle->base_vector[2] = middle_point[2];
140 bnd_circle->inc_angle = sum_angles_no_check(inc_angle, SIN_COS_TOL);
141 bnd_circle->sq_crd = DBL_MAX;
142}
143
144// computes a) the angle between the middle point of the edge with the middle
145// point of the polygon
146// b) half of the angle between between the two vertices of the edge
147// returns the sum of both angles
149 double * restrict a, double * restrict b, double * restrict middle_point) {
150
151 double edge_middle_point[3] = {a[0] + b[0], a[1] + b[1], a[2] + b[2]};
152 normalise_vector(edge_middle_point);
153
154 struct sin_cos_angle t1 = get_vector_angle_2(edge_middle_point, a);
155 struct sin_cos_angle t2 = get_vector_angle_2(edge_middle_point, middle_point);
156
157 return sum_angles_no_check(t1, t2);
158}
159
161 struct bounding_circle * bnd_circle) {
162
163 double middle_point[3];
164
165 middle_point[0] = 0.0;
166 middle_point[1] = 0.0;
167 middle_point[2] = 0.0;
168
169 size_t num_corners = cell.num_corners;
170
171 // compute the coordinates in rad and 3d
172 for (size_t i = 0; i < num_corners; ++i) {
173
174 middle_point[0] += cell.coordinates_xyz[i][0];
175 middle_point[1] += cell.coordinates_xyz[i][1];
176 middle_point[2] += cell.coordinates_xyz[i][2];
177 }
178
179 normalise_vector(middle_point);
180
181 // compute the angle required for the bounding circle
182 struct sin_cos_angle max_angle =
184 cell.coordinates_xyz[num_corners-1],
185 cell.coordinates_xyz[0], middle_point);
186
187 for (size_t i = 0; i < num_corners-1; ++i) {
188
189 struct sin_cos_angle curr_angle =
191 cell.coordinates_xyz[i], cell.coordinates_xyz[i+1], middle_point);
192
193 if (compare_angles(max_angle, curr_angle) < 0) max_angle = curr_angle;
194 }
195
196 bnd_circle->base_vector[0] = middle_point[0];
197 bnd_circle->base_vector[1] = middle_point[1];
198 bnd_circle->base_vector[2] = middle_point[2];
199
200 bnd_circle->inc_angle = sum_angles_no_check(max_angle, SIN_COS_TOL);
201 bnd_circle->sq_crd = DBL_MAX;
202}
203
205 struct bounding_circle * extent_b) {
206
207 struct sin_cos_angle base_vector_angle =
208 get_vector_angle_2(extent_a->base_vector, extent_b->base_vector);
209
210 struct sin_cos_angle tmp_angle, inc_angle_sum;
211 int big_sum =
212 sum_angles(extent_a->inc_angle, extent_b->inc_angle, &tmp_angle);
213
214 if (big_sum ||
215 sum_angles(tmp_angle, SIN_COS_TOL, &inc_angle_sum))
216 return 1;
217
218 return compare_angles(base_vector_angle, inc_angle_sum) <= 0;
219}
void yac_get_cell_circumscribe_circle_reg_quad(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:27
void yac_get_cell_circumscribe_circle_unstruct_triangle(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:45
void yac_get_cell_bounding_circle_unstruct_triangle(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:93
void yac_get_cell_bounding_circle_reg_quad(double a[3], double b[3], double c[3], struct bounding_circle *bnd_circle)
Definition bnd_circle.c:36
void yac_get_cell_bounding_circle(struct yac_grid_cell cell, struct bounding_circle *bnd_circle)
Definition bnd_circle.c:160
int yac_extents_overlap(struct bounding_circle *extent_a, struct bounding_circle *extent_b)
Definition bnd_circle.c:204
static struct sin_cos_angle compute_edge_inc_angle(double *restrict a, double *restrict b, double *restrict middle_point)
Definition bnd_circle.c:148
static double get_sin_vector_angle(double a[3], double b[3])
Definition bnd_circle.c:80
static struct sin_cos_angle get_vector_angle_2(double const a[3], double const b[3])
Definition geometry.h:415
static void crossproduct_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
static struct sin_cos_angle sum_angles_no_check(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:524
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:332
static const struct sin_cos_angle SIN_COS_TOL
Definition geometry.h:45
static int compare_angles(struct sin_cos_angle a, struct sin_cos_angle b)
Definition geometry.h:428
static void normalise_vector(double v[])
Definition geometry.h:689
static int sum_angles(struct sin_cos_angle a, struct sin_cos_angle b, struct sin_cos_angle *restrict sum)
Definition geometry.h:534
static struct sin_cos_angle sin_cos_angle_new(double sin, double cos)
Definition geometry.h:405
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
Definition geometry.h:61
double base_vector[3]
Definition geometry.h:59
double sq_crd
Definition geometry.h:64
size_t num_corners
Definition grid_cell.h:21
double(* coordinates_xyz)[3]
Definition grid_cell.h:19