YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_partial_areas.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#include <stdlib.h>
6#include <stdio.h>
7#include <math.h>
8
9#include "grids/basic_grid.h"
10#include "clipping.h"
11#include "area.h"
12#include "tests.h"
13#include "geometry.h"
14#include "test_common.h"
15
21static void utest_check_partial_areas_deg(
22 double * tgt_lon, double * tgt_lat,
23 enum yac_edge_type * tgt_edges, int tgt_count,
24 double ** src_lons, double ** src_lats,
25 enum yac_edge_type ** src_edges, int * src_counts, size_t nSourceCells);
26static void utest_check_partial_areas_rad(
27 double * tgt_lon, double * tgt_lat,
28 enum yac_edge_type * tgt_edges, int tgt_count,
29 double ** src_lons, double ** src_lats,
30 enum yac_edge_type ** src_edges, int * src_counts, size_t nSourceCells);
31
32int main (void) {
33
34 enum yac_edge_type gc_edges[] = {YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
35 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
36 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE,
37 YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE};
38 enum yac_edge_type latlon_edges[] = {YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE,
39 YAC_LAT_CIRCLE_EDGE, YAC_LON_CIRCLE_EDGE};
40
41 /* 4 source elements overlapping a target cell
42
43 Target Cell
44
45 0.0 (lon) [1]
46 0.7 (lat)
47 / \
48 / \
49 -0.5 [2] / \ 0.5 [0]
50 0.0 \ / 0.0
51 \ /
52 \ /
53 0.0 [3]
54 -0.5
55 */
56 /* source cell test data
57
58 0.0 1.0 (lon)
59 1.8 1.8 (lat)
60 +-------+--------+
61 | | |
62 | 1 | 0 |
63 | | |
64 | | | 1.0
65 +-------+--------+ 0.6
66 | 0.0 |
67 | 0.6 |
68 | |
69 | 2 | 3 |
70 | | |
71 +-------+--------+ 1.0
72 -0.6
73 */
74 {
75 enum {nSourceCells = 4};
76 utest_check_partial_areas_deg(
77 (double[]){0.5, 0.0, -0.5, 0.0},
78 (double[]){0.0, 0.7, 0.0, -0.5}, gc_edges, 4,
79 (double*[nSourceCells]){(double[]){1.0, 0.0, 0.0, 1.0},
80 (double[]){0.0, -1.0, -1.0, 0.0},
81 (double[]){0.0, -1.0, -1.0, 0.0},
82 (double[]){1.0, 0.0, 0.0, 1.0}},
83 (double*[nSourceCells]){(double[]){1.8, 1.8, 0.6, 0.6},
84 (double[]){1.8, 1.8, 0.6, 0.6},
85 (double[]){0.6, 0.6, -0.6, -0.6},
86 (double[]){0.6, 0.6, -0.6, -0.6}},
87 (enum yac_edge_type*[nSourceCells])
88 {gc_edges, gc_edges, gc_edges, gc_edges},
89 (int[nSourceCells]){4, 4, 4, 4}, nSourceCells);
90 }
91
92 /* Next test is over the pole
93
94 Target cell located over the pole o with
95 vertices x and source cells 0 to 3 arranged as
96
97
98 135.0 90.0 45.0 (lon)
99 89.0 88.5 (lat)
100 ------------------
101 | | |
102 | | |
103 | 1x |
104 | | 0 |
105-180.0 ---x----o----x---- 0.0 pole @ -180.0; 90.0
106 | 2 | |
107 | | |
108 | 3x |
109 | | |
110 | | |
111-135.0 ------------------ -45.0
112 -90.0 -0.6
113
114 */
115 {
116 enum {nSourceCells = 4};
117 utest_check_partial_areas_deg(
118 (double[]){ 0.0, 90.0, -180.0, -90.0},
119 (double[]){89.5, 89.5, 89.5, 89.5}, gc_edges, 4,
120 (double*[nSourceCells]){(double[]){45.0, 90.0, -180.0, 0.0},
121 (double[]){90.0, 135.0, -180.0, -180.0},
122 (double[]){-180.0, -180.0, -135.0, -90.0},
123 (double[]){ 0.0, -180.0, -90.0, -45.0}},
124 (double*[nSourceCells]){(double[]){88.5, 89.0, 90.0, 89.0},
125 (double[]){89.0, 88.5, 89.0, 90.0},
126 (double[]){90.0, 89.0, 88.5, 89.0},
127 (double[]){89.0, 90.0, 89.0, 88.5}},
128 (enum yac_edge_type*[nSourceCells])
129 {gc_edges, gc_edges, gc_edges, gc_edges},
130 (int[nSourceCells]){4, 4, 4, 4}, nSourceCells);
131 }
132
133 // checking overlap between identical cells
134 {
135 enum {nSourceCells = 1};
136 utest_check_partial_areas_deg(
137 (double[]){-180.0, -150.0, -150.0, -180.0},
138 (double[]){ -90.0, -90.0, -60.0, -60.0}, gc_edges, 4,
139 (double*[nSourceCells]){(double[]){-180.0, -150.0, -150.0, -180.0}},
140 (double*[nSourceCells]){(double[]){ -90.0, -90.0, -60.0, -60.0}},
141 (enum yac_edge_type*[nSourceCells]){gc_edges},
142 (int[nSourceCells]){4}, nSourceCells);
143 }
144 {
145 enum {nSourceCells = 1};
146 utest_check_partial_areas_deg(
147 (double[]){-180.0, -150.0, -150.0, -180.0},
148 (double[]){ 60.0, 60.0, 90.0, 90.0}, gc_edges, 4,
149 (double*[nSourceCells]){(double[]){-180.0, -150.0, -150.0, -180.0}},
150 (double*[nSourceCells]){(double[]){ 60.0, 60.0, 90.0, 90.0}},
151 (enum yac_edge_type*[nSourceCells]){gc_edges},
152 (int[nSourceCells]){4}, nSourceCells);
153 }
154
155 // one out of three source cells completely overlaps with the target cell
156 {
157 enum {nSourceCells = 3};
158 utest_check_partial_areas_deg(
159 (double[]){25.0, 30.0, 30.0, 25.0},
160 (double[]){75.0, 75.0, 80.0, 80.0}, gc_edges, 4,
161 (double*[nSourceCells]){(double[]){30.0, 35.0, 35.0, 30.0},
162 (double[]){25.0, 30.0, 30.0, 25.0},
163 (double[]){20.0, 25.0, 25.0, 20.0}},
164 (double*[nSourceCells]){(double[]){75.0, 75.0, 80.0, 80.0},
165 (double[]){75.0, 75.0, 80.0, 80.0},
166 (double[]){75.0, 75.0, 80.0, 80.0}},
167 (enum yac_edge_type*[nSourceCells])
168 {gc_edges, gc_edges, gc_edges},
169 (int[nSourceCells]){4, 4, 4}, nSourceCells);
170 }
171
172 // two cell are very close to each other
173 {
174 struct yac_grid_cell TargetCell =
175 generate_cell_rad((double[]){0.098707396293595859,
176 0.117826295355522680,
177 0.098174496441371703},
178 (double[]){0.333033372840672190,
179 0.345484775276814930,
180 0.353335609258583650}, gc_edges, 3);
181 struct yac_grid_cell SourceCell =
182 generate_cell_rad((double[]){0.073631077818510776,
183 0.098174770424681035,
184 0.098174770424681035,
185 0.073631077818510776},
186 (double[]){0.343611696486383620,
187 0.343611696486383620,
188 0.368155389092553910,
189 0.368155389092553910}, gc_edges, 4);
190
191 double overlap_area;
192 yac_compute_overlap_areas(
193 1, &SourceCell, TargetCell, &overlap_area);
194
195 if (fabs(overlap_area) > 1e-9)
196 PUT_ERR("ERROR in yac_compute_overlap_areas");
197
198 yac_free_grid_cell(&SourceCell);
199 yac_free_grid_cell(&TargetCell);
200 }
201
202 // overlapping lon-lat cells
203 {
204 enum {nSourceCells = 4};
205 utest_check_partial_areas_deg(
206 (double[]){-0.5, 0.5, 0.5, -0.5},
207 (double[]){-0.5, -0.5, 0.5, 0.5}, latlon_edges, 4,
208 (double*[nSourceCells]){(double[]){-1, 0, 0, -1},
209 (double[]){ 0, 1, 1, 0},
210 (double[]){-1, 0, 0, -1},
211 (double[]){0, 1, 1, 0}},
212 (double*[nSourceCells]){(double[]){-1, -1, 0, 0},
213 (double[]){-1, -1, 0, 0},
214 (double[]){ 0, 0, 1, 1},
215 (double[]){ 0, 0, 1, 1}},
216 (enum yac_edge_type*[nSourceCells])
217 {latlon_edges, latlon_edges, latlon_edges, latlon_edges},
218 (int[nSourceCells]){4, 4, 4, 4}, nSourceCells);
219 }
220
221 // overlapping lon-lat cells
222 {
223 enum {nSourceCells = 1};
224 utest_check_partial_areas_deg(
225 (double[]){-0.5, 0.5, 0.5, -0.5},
226 (double[]){-0.5, -0.5, 0.5, 0.5}, latlon_edges, 4,
227 (double*[nSourceCells]){(double[]){-1, 1, 1, -1}},
228 (double*[nSourceCells]){(double[]){-1, -1, 1, 1}},
229 (enum yac_edge_type*[nSourceCells])
230 {latlon_edges},
231 (int[nSourceCells]){4}, nSourceCells);
232 }
233
234 // lon-lat cells overlapping with gc cell
235 {
236 enum {nSourceCells = 1};
237 utest_check_partial_areas_deg(
238 (double[]){-0.5, 0.5, 0.0},
239 (double[]){-0.5, -0.5, 0.5}, gc_edges, 3,
240 (double*[nSourceCells]){(double[]){-1, 1, 1, -1}},
241 (double*[nSourceCells]){(double[]){-1, -1, 1, 1}},
242 (enum yac_edge_type*[nSourceCells])
243 {latlon_edges},
244 (int[nSourceCells]){4}, nSourceCells);
245 }
246
247 // test case provided by a user
248 {
249 enum {nSourceCells = 4};
250 utest_check_partial_areas_rad(
251 (double[]){2.5034566458293663,
252 2.5280003384355365,
253 2.5280003384355365,
254 2.5034566458293663},
255 (double[]){0.5890486225480862,
256 0.5890486225480862,
257 0.6135923151542565,
258 0.6135923151542565}, latlon_edges, 4,
259 (double*[nSourceCells]){(double[]){2.5132750564642796,
260 2.5132751216401110,
261 2.5530451220055963},
262 (double[]){2.5530451220055963,
263 2.5132751216401110,
264 2.5514531161521616},
265 (double[]){2.4735050741253155,
266 2.5132751216401110,
267 2.5132750564642796},
268 (double[]){2.4750971948235070,
269 2.5132751216401110,
270 2.4735050741253155}},
271 (double*[nSourceCells]){(double[]){0.6227096831849960,
272 0.5888365165503746,
273 0.6044728875052012},
274 (double[]){0.6044728875052012,
275 0.5888365165503746,
276 0.5706926293017953},
277 (double[]){0.6044725603601205,
278 0.5888365165503746,
279 0.6227096831849960},
280 (double[]){0.5706923461529829,
281 0.5888365165503746,
282 0.6044725603601205}},
283 (enum yac_edge_type*[nSourceCells])
284 {gc_edges, gc_edges, gc_edges, gc_edges},
285 (int[nSourceCells]){3, 3, 3, 3}, nSourceCells);
286 }
287
288 // test case provided by a user
289 {
290 enum {nSourceCells = 3};
291 utest_check_partial_areas_rad(
292 (double[]){0.6626797003665970,
293 0.6872233929727672,
294 0.6872233929727672,
295 0.6626797003665970},
296 (double[]){1.2271846303085130,
297 1.2271846303085130,
298 1.2517283229146832,
299 1.2517283229146832}, latlon_edges, 4,
300 (double*[nSourceCells]){(double[]){0.6910166804158336,
301 0.6283648515564006,
302 0.7428780447168533},
303 (double[]){0.5657152198049494,
304 0.6283648515564006,
305 0.6910166804158336},
306 (double[]){0.6283674880262698,
307 0.5657152198049494,
308 0.6910166804158336}},
309 (double*[nSourceCells]){(double[]){1.2515918327919406,
310 1.2207207913058662,
311 1.2190404189581385},
312 (double[]){1.2515932566235379,
313 1.2207207913058662,
314 1.2515918327919406},
315 (double[]){1.2831450892755811,
316 1.2515932566235379,
317 1.2515918327919406}},
318 (enum yac_edge_type*[nSourceCells])
319 {gc_edges, gc_edges, gc_edges},
320 (int[nSourceCells]){3, 3, 3}, nSourceCells);
321 }
322
323 // test case provided by a user
324 {
325 enum {nSourceCells = 7};
326 utest_check_partial_areas_rad(
327 (double[]){-2.1798479726998332,
328 -2.1991148570983681,
329 -2.2190472361409284},
330 (double[]){-0.0311624518327736,
331 0.0000000043633749,
332 -0.0338654325037920}, gc_edges, 3,
333 (double*[nSourceCells]){(double[]){ 4.0987966652304335,
334 4.1233403578366037,
335 4.1233403578366037,
336 4.0987966652304335},
337 (double[]){ 4.0742529726242633,
338 4.0987966652304335,
339 4.0987966652304335,
340 4.0742529726242633},
341 (double[]){ 4.0742529726242633,
342 4.0987966652304335,
343 4.0987966652304335,
344 4.0742529726242633},
345 (double[]){ 4.0987966652304335,
346 4.1233403578366037,
347 4.1233403578366037,
348 4.0987966652304335},
349 (double[]){ 4.0742529726242633,
350 4.0987966652304335,
351 4.0987966652304335,
352 4.0742529726242633},
353 (double[]){ 4.0497092800180932,
354 4.0742529726242633,
355 4.0742529726242633,
356 4.0497092800180932},
357 (double[]){ 4.0497092800180932,
358 4.0742529726242633,
359 4.0742529726242633,
360 4.0497092800180932}},
361 (double*[nSourceCells]){(double[]){-0.0490873852123405,
362 -0.0490873852123405,
363 -0.0245436926061703,
364 -0.0245436926061703},
365 (double[]){-0.0490873852123405,
366 -0.0490873852123405,
367 -0.0245436926061703,
368 -0.0245436926061703},
369 (double[]){ 0.0000000000000000,
370 0.0000000000000000,
371 0.0245436926061703,
372 0.0245436926061703},
373 (double[]){-0.0245436926061703,
374 -0.0245436926061703,
375 0.0000000000000000,
376 0.0000000000000000},
377 (double[]){-0.0245436926061703,
378 -0.0245436926061703,
379 0.0000000000000000,
380 0.0000000000000000},
381 (double[]){-0.0245436926061703,
382 -0.0245436926061703,
383 0.0000000000000000,
384 0.0000000000000000},
385 (double[]){-0.0490873852123405,
386 -0.0490873852123405,
387 -0.0245436926061703,
388 -0.0245436926061703}},
389 (enum yac_edge_type*[nSourceCells])
390 {latlon_edges, latlon_edges, latlon_edges, latlon_edges,
391 latlon_edges, latlon_edges, latlon_edges},
392 (int[nSourceCells]){4, 4, 4, 4, 4, 4, 4}, nSourceCells);
393 }
394
395 // test case provided by a user
396 {
397 enum {nSourceCells = 4};
398 utest_check_partial_areas_rad(
399 (double[]){ 3.4361169648638361,
400 3.4606606574700067,
401 3.4606606574700067,
402 3.4361169648638361},
403 (double[]){-0.0245436926061703,
404 -0.0245436926061703,
405 0.0000000000000000,
406 0.0000000000000000}, latlon_edges, 4,
407 (double*[nSourceCells]){(double[]){-2.8075010104343816,
408 -2.8274333887952685,
409 -2.8467002736213392},
410 (double[]){-2.8081665040027728,
411 -2.8274333887952685,
412 -2.7882305020025413},
413 (double[]){-2.7882305020025413,
414 -2.8274333887952685,
415 -2.8075010104343816},
416 (double[]){-2.8467002736213392,
417 -2.8274333887952685,
418 -2.8666362753717674}},
419 (double*[nSourceCells]){(double[]){-0.0338654341705874,
420 0.0000000011622542,
421 -0.0311624540855508},
422 (double[]){ 0.0311624566805337,
423 0.0000000011622542,
424 -0.0026744729717723},
425 (double[]){-0.0026744729717723,
426 0.0000000011622542,
427 -0.0338654341705874},
428 (double[]){-0.0311624540855508,
429 0.0000000011622542,
430 0.0026744751289656}},
431 (enum yac_edge_type*[nSourceCells])
432 {gc_edges, gc_edges, gc_edges, gc_edges},
433 (int[nSourceCells]){3, 3, 3, 3}, nSourceCells);
434 }
435
436 // test case provided by a user
437 // (the overlap between the first source cell and the target cell
438 // is very small, which is why it is potentially dropped by YAC)
439 {
440 enum {nSourceCells = 3};
441 utest_check_partial_areas_rad(
442 (double[]){ 3.1170489609836229,
443 3.1415926535897931,
444 3.1415926535897931,
445 3.1170489609836229},
446 (double[]){-1.0062913968529805,
447 -1.0062913968529805,
448 -0.9817477042468103,
449 -0.9817477042468103}, latlon_edges, 4,
450 (double*[nSourceCells]){(double[]){ 3.1415926521557882,
451 3.1415926520689506,
452 -3.0774466652769172},
453 (double[]){ 3.0774466624490269,
454 3.1415926520689506,
455 3.1415926521557882},
456 (double[]){ 3.0808931330453300,
457 3.0774466624490269,
458 3.1415926521557882}},
459 (double*[nSourceCells]){(double[]){-0.9805860393425995,
460 -1.0172219678978514,
461 -0.9979427097227050},
462 (double[]){-0.9979427096430752,
463 -1.0172219678978514,
464 -0.9805860393425995},
465 (double[]){-0.9613498550843829,
466 -0.9979427096430752,
467 -0.9805860393425995}},
468 (enum yac_edge_type*[nSourceCells])
469 {gc_edges, gc_edges, gc_edges},
470 (int[nSourceCells]){3, 3, 3}, nSourceCells);
471 }
472
473 // overlap of a triangle with a lot of very narrow lon lat cells
474 {
475 enum {nSourceCells = 198};
476 double ** src_lons = xmalloc(2 * nSourceCells * sizeof(*src_lons));
477 double ** src_lats = src_lons + nSourceCells;
478 enum yac_edge_type ** src_edges =
479 xmalloc(nSourceCells * sizeof(*src_edges));
480 int * src_counts = xmalloc(nSourceCells * sizeof(*src_counts));
481
482 src_lons[0] = xmalloc(2 * 4 * nSourceCells * sizeof(**src_lons));
483 double dx = ((49.2 - 22.8) /((double)nSourceCells));
484 for (size_t i = 0; i < nSourceCells; ++i) {
485 src_lons[i] = src_lons[0] + i * 4;
486 src_lats[i] = src_lons[0] + i * 4 + 4 * nSourceCells;
487 src_lons[i][0] = 22.8000 + dx*((double)(i+0));
488 src_lons[i][1] = 22.8000 + dx*((double)(i+1));
489 src_lons[i][2] = 22.8000 + dx*((double)(i+1));
490 src_lons[i][3] = 22.8000 + dx*((double)(i+0));
491 src_lats[i][0] = 89.7333;
492 src_lats[i][1] = 89.7333;
493 src_lats[i][2] = 89.8667;
494 src_lats[i][3] = 89.8667;
495 src_edges[i] = latlon_edges;
496 src_counts[i] = 4;
497 }
498
499 utest_check_partial_areas_deg(
500 (double[]){22.8313, 49.1687, 36.0000},
501 (double[]){89.7418, 89.7418, 89.8335}, gc_edges, 3,
502 src_lons, src_lats, src_edges, src_counts, nSourceCells);
503
504 free(src_counts);
505 free(src_edges);
506 free(src_lons[0]);
507 free(src_lons);
508 }
509
510 // test case provided by a user
511 {
512 enum {nSourceCells = 4};
513 utest_check_partial_areas_rad(
514 (double[]){ 5.0069132916587327,
515 5.0130492148102759,
516 5.0130492148102759,
517 5.0069132916587327},
518 (double[]){-1.2946797849754812,
519 -1.2946797849754812,
520 -1.2885438618239387,
521 -1.2885438618239387}, latlon_edges, 4,
522 (double*[nSourceCells]){(double[]){ 4.9999997254688875,
523 4.9999997254688875,
524 5.0078539201557488,
525 5.0078539201557488},
526 (double[]){ 5.0078539201557488,
527 5.0078539201557488,
528 5.0157075822103927,
529 5.0157075822103927},
530 (double[]){ 5.0078539201557488,
531 5.0078539201557488,
532 5.0157075822103927,
533 5.0157075822103927},
534 (double[]){ 4.9999997254688875,
535 4.9999997254688875,
536 5.0078539201557488,
537 5.0078539201557488}},
538 (double*[nSourceCells]){(double[]){-1.2878426515089207,
539 -1.2946756570757916,
540 -1.2946797849754812,
541 -1.2878469125666649},
542 (double[]){-1.2946797849754812,
543 -1.3015127905423520,
544 -1.3015170516000960,
545 -1.2946840460332254},
546 (double[]){-1.2878469125666649,
547 -1.2946797849754812,
548 -1.2946840460332254,
549 -1.2878513067824635},
550 (double[]){-1.2946756570757916,
551 -1.3015085294846078,
552 -1.3015127905423520,
553 -1.2946797849754812}},
554 (enum yac_edge_type*[nSourceCells])
555 {gc_edges, gc_edges, gc_edges, gc_edges},
556 (int[nSourceCells]){4, 4, 4, 4}, nSourceCells);
557 }
558
559 // test case provided by a user
560 {
561 enum {nSourceCells = 2};
562 utest_check_partial_areas_rad(
563 (double[]){-2.8460897445189417,
564 -2.8659852121839959,
565 -2.8853055611946927},
566 (double[]){-0.0961602003314924,
567 -0.0622820530745658,
568 -0.0933150814318732}, gc_edges, 3,
569 (double*[nSourceCells]){(double[]){-2.8655528147976832,
570 -2.9038195735214831,
571 -2.8254194450162888},
572 (double[]){-2.9440830021466091,
573 -2.9038195735214831,
574 -2.8655528147976832}},
575 (double*[nSourceCells]){(double[]){-0.0615856127966505,
576 -0.1228534567110898,
577 -0.1292806814915037},
578 (double[]){-0.0558979744839555,
579 -0.1228534567110898,
580 -0.0615856127966505}},
581 (enum yac_edge_type*[nSourceCells]){gc_edges, gc_edges},
582 (int[nSourceCells]){3, 3}, nSourceCells);
583 }
584
585 // test case provided by a user
586 {
587 enum {nSourceCells = 6};
588 utest_check_partial_areas_rad(
589 (double[]){3.0127596408133406,
590 2.9702182441873108,
591 2.9943996418061438},
592 (double[]){0.5498596133493898,
593 0.5469967117356207,
594 0.5140000458840724}, gc_edges, 3,
595 (double*[nSourceCells]){(double[]){2.9198084624304621,
596 2.9720969737232750,
597 3.0078074923914948},
598 (double[]){2.8889398582596519,
599 2.9720969737232750,
600 2.9198084624304621},
601 (double[]){3.0187962726049782,
602 2.9720969737232750,
603 2.9378118571213374},
604 (double[]){2.9378118571213374,
605 2.9720969737232750,
606 2.8889398582596519},
607 (double[]){3.0565195473916882,
608 2.9720969737232750,
609 3.0187962726049782},
610 (double[]){3.0078074923914948,
611 2.9720969737232750,
612 3.0565195473916882}},
613 (double*[nSourceCells]){(double[]){0.6123031304013200,
614 0.5471403770103002,
615 0.6192726969308253},
616 (double[]){0.5392518296263771,
617 0.5471403770103002,
618 0.6123031304013200},
619 (double[]){0.4805223860987507,
620 0.5471403770103002,
621 0.4745836829538209},
622 (double[]){0.4745836829538209,
623 0.5471403770103002,
624 0.5392518296263771},
625 (double[]){0.5519553785110132,
626 0.5471403770103002,
627 0.4805223860987507},
628 (double[]){0.6192726969308253,
629 0.5471403770103002,
630 0.5519553785110132}},
631 (enum yac_edge_type*[nSourceCells])
632 {gc_edges, gc_edges, gc_edges, gc_edges, gc_edges, gc_edges},
633 (int[nSourceCells]){3, 3, 3, 3, 3, 3}, nSourceCells);
634 }
635
636 { // test the barycenter returned by yac_compute_overlap_info for a
637 // target cell with more than 4 edges
638 struct yac_grid_cell TargetCell =
639 generate_cell_deg((double[]){-0.5, 0.5, 0.75, 0.5, -0.5, -0.75},
640 (double[]){-0.5, -0.5, 0.0 , 0.5, 0.5, 0.0},
641 gc_edges, 6);
642 struct yac_grid_cell SourceCells[] = {
643 generate_cell_deg((double[]){0.3, -0.3, 0.0},
644 (double[]){0.3, 0.3, -0.3}, gc_edges, 3),
645 generate_cell_deg((double[]){-0.4, 0.4, 0.4, -0.4},
646 (double[]){-0.4, -0.4, 0.4, 0.4}, gc_edges, 4),
647 generate_cell_deg((double[]){-0.6, 0.6, 0.8, 0.6, -0.6, -0.8},
648 (double[]){-0.6, -0.6, 0.0, 0.6, 0.6, 0.0},
649 gc_edges, 6)};
650 enum {nSourceCells = sizeof(SourceCells)/sizeof(SourceCells[0])};
651
652 double ref_area[nSourceCells];
653 double ref_barycenter[nSourceCells][3];
654 for (int i = 0; i < nSourceCells; ++i)
655 for (int j = 0; j < 3; ++j)
656 ref_barycenter[i][j] = 0.0;
657 {
658 ref_area[0] =
659 yac_huiliers_area_info(SourceCells[0], ref_barycenter[0], 1.0);
660 ref_area[1] = yac_huiliers_area(SourceCells[1]);
661 ref_area[2] = yac_huiliers_area(TargetCell);
662 LLtoXYZ(0, 0, ref_barycenter[1]);
663 LLtoXYZ(0, 0, ref_barycenter[2]);
664 }
665
666 double overlap_areas[nSourceCells];
667 double overlap_barycenters[nSourceCells][3];
668 yac_compute_overlap_info(
669 nSourceCells, SourceCells, TargetCell,
670 overlap_areas, &overlap_barycenters[0]);
671
672 for (int i = 0; i < nSourceCells; ++i) {
673 double area_tolerance = yac_huiliers_area(SourceCells[i]) * 1e-6;
674 if (fabs(ref_area[i] - overlap_areas[i]) > area_tolerance)
675 PUT_ERR("ERROR in yac_compute_overlap_info (area)");
676 if (overlap_areas[i] > area_tolerance) {
677 normalise_vector(ref_barycenter[i]);
678 normalise_vector(overlap_barycenters[i]);
679 if(get_vector_angle(
680 ref_barycenter[i], overlap_barycenters[i]) > 1e-9)
681 PUT_ERR("ERROR in yac_compute_overlap_info (barycenter)");
682 }
683 }
684
685 for (size_t i = 0; i < nSourceCells; ++i)
686 yac_free_grid_cell(&SourceCells[i]);
687 yac_free_grid_cell(&TargetCell);
688 }
689
690 { // test the area and barycenter returned by
691 // yac_compute_overlap_info for a triangle target cell
692
693 struct yac_grid_cell TargetCells[] =
694 // normal triangle
695 {generate_cell_deg((double[]){-0.5, 0.5, -0.5},
696 (double[]){-0.5, 0.5, 0.5}, gc_edges, 3),
697 // triangle with a duplicated vertex to simulate a degenerated cell
698 generate_cell_deg((double[]){-0.5, 0.5, -0.5, -0.5},
699 (double[]){-0.5, 0.5, 0.5, 0.5}, gc_edges, 4)};
700 enum {nTargetCells = sizeof(TargetCells)/sizeof(TargetCells[0])};
701
702 struct yac_grid_cell SourceCells[] = {
703 generate_cell_deg((double[]){-0.6, 0.7, -0.6},
704 (double[]){-0.7, 0.6, 0.6}, gc_edges, 3),
705 generate_cell_deg((double[]){-0.4, 0.3, -0.4},
706 (double[]){-0.3, 0.4, 0.4}, gc_edges, 3),
707 generate_cell_deg((double[]){-0.5, 0.5, -0.5},
708 (double[]){-0.5, 0.5, 0.5}, gc_edges, 3),
709 generate_cell_deg((double[]){0.0, 1.0, 1.0, 0.0},
710 (double[]){0.0, 0.0, 1.0, 1.0}, latlon_edges, 4),
711 generate_cell_deg((double[]){-1.0, 0.0, 0.0, -1.0},
712 (double[]){-1.0, -1.0, 1.0, 1.0}, latlon_edges, 4),
713 generate_cell_deg((double[]){1.0, 2.0, 2.0, 1.0},
714 (double[]){1.0, 1.0, 2.0, 2.0}, latlon_edges, 4)};
715 enum {nSourceCells = sizeof(SourceCells)/sizeof(SourceCells[0])};
716
717 double ref_area[nSourceCells];
718 double ref_barycenter[nSourceCells][3];
719 for (int i = 0; i < nSourceCells; ++i)
720 for (int j = 0; j < 3; ++j)
721 ref_barycenter[i][j] = 0.0;
722 {
723 ref_area[0] = yac_huiliers_area_info(TargetCells[0], ref_barycenter[0], 1.0);
724 ref_area[1] = yac_huiliers_area_info(SourceCells[1], ref_barycenter[1], 1.0);
725 ref_area[2] = yac_huiliers_area_info(TargetCells[0], ref_barycenter[2], 1.0);
726 {
727 double intersection[3], intersection_lon, intersection_lat;
728 if (!intersect(YAC_GREAT_CIRCLE_EDGE, -0.5, 0.5, 0.5, 0.5,
729 YAC_GREAT_CIRCLE_EDGE, 0.0, 0.0, 0.0, 1.0, intersection))
730 return EXIT_FAILURE;
731 XYZtoLL(intersection, &intersection_lon, &intersection_lat);
732 intersection_lon /= YAC_RAD;
733 intersection_lat /= YAC_RAD;
734 struct yac_grid_cell overlap_cell =
735 generate_cell_deg((double[]){0.0, 0.5, intersection_lon},
736 (double[]){0.0, 0.5, intersection_lat},
737 gc_edges, 3);
738 ref_area[3] =
739 yac_huiliers_area_info(overlap_cell, ref_barycenter[3], 1.0);
740 yac_free_grid_cell(&overlap_cell);
741 overlap_cell =
742 generate_cell_deg((double[]){-0.5, -0.5, 0.0, intersection_lon},
743 (double[]){ 0.5, -0.5, 0.0, intersection_lat},
744 gc_edges, 4);
745 ref_area[4] =
746 yac_huiliers_area_info(overlap_cell, ref_barycenter[4], 1.0);
747 yac_free_grid_cell(&overlap_cell);
748 }
749 ref_area[5] = 0.0;
750 ref_barycenter[5][0] = 0.0;
751 ref_barycenter[5][1] = 0.0;
752 ref_barycenter[5][2] = 0.0;
753 }
754
755 for (int i = 0; i < nTargetCells; ++i) {
756
757 double overlap_areas[nSourceCells];
758 double overlap_barycenters[nSourceCells][3];
759 yac_compute_overlap_info(
760 nSourceCells, SourceCells, TargetCells[i],
761 overlap_areas, &overlap_barycenters[0]);
762
763 for (int i = 0; i < nSourceCells; ++i) {
764 double area_tolerance = yac_huiliers_area(SourceCells[i]) * 1e-6;
765 if (fabs(ref_area[i] - overlap_areas[i]) > area_tolerance)
766 PUT_ERR("ERROR in yac_compute_overlap_info (area)");
767 if (overlap_areas[i] > area_tolerance) {
768 normalise_vector(ref_barycenter[i]);
769 normalise_vector(overlap_barycenters[i]);
770 if(get_vector_angle(
771 ref_barycenter[i], overlap_barycenters[i]) > 1e-9)
772 PUT_ERR("ERROR in yac_compute_overlap_info (barycenter)");
773 }
774 }
775 }
776 for (int i = 0; i < nSourceCells; ++i)
777 yac_free_grid_cell(&SourceCells[i]);
778 for (int i = 0; i < nTargetCells; ++i)
779 yac_free_grid_cell(&TargetCells[i]);
780 }
781
782 // cleanup
783 yac_compute_overlap_buf_free();
784
785 return TEST_EXIT_CODE;
786}
787
788static void utest_check_partial_areas(
789 double * tgt_lon, double * tgt_lat,
790 enum yac_edge_type * tgt_edges, int tgt_count,
791 double ** src_lons, double ** src_lats,
792 enum yac_edge_type ** src_edges, int * src_counts, size_t nSourceCells,
793 struct yac_grid_cell (*generate_cell)(
794 double*, double*, enum yac_edge_type*,size_t)) {
795
796 /* For the test at the Equator we cannot get better than 1.0e-11, while
797 for the Pole test we are somewhat better.
798
799 tolerance test for the deviation from 1 of the sum of weights: 1.0e-11
800
801 */
802 double const epsilon = 1.0e-10; // relative precision
803
804 struct yac_grid_cell TargetCell =
805 generate_cell(tgt_lon, tgt_lat, tgt_edges, tgt_count);
806
807 struct yac_grid_cell * SourceCells =
808 xmalloc(nSourceCells * sizeof(*SourceCells));
809 for (size_t i = 0; i < nSourceCells; ++i)
810 SourceCells[i] =
811 generate_cell(
812 src_lons[i], src_lats[i], src_edges[i], src_counts[i]);
813
814 double tgt_area = yac_huiliers_area(TargetCell);
815 double area_tolerance = MAX(tgt_area * 1e-6, 1e-10);
816
817 double partial_areas[nSourceCells];
818 yac_compute_overlap_areas(
819 nSourceCells, SourceCells, TargetCell, partial_areas);
820
821 double partial_areas_sum = 0.0;
822 double weights[nSourceCells];
823 for (size_t i = 0; i < nSourceCells; ++i) {
824 partial_areas_sum += partial_areas[i];
825 weights[i] = partial_areas[i] / tgt_area;
826 }
827
828 // check partial areas
829 if (fabs(partial_areas_sum - tgt_area) > area_tolerance)
830 PUT_ERR("ERROR in sum of partial areas\n");
831
832 // test yac_correct_weights
833 yac_correct_weights(nSourceCells, weights);
834 double weights_sum = 0.0;
835 for (size_t i = 0; i < nSourceCells; ++i) weights_sum += weights[i];
836 if ( fabs(weights_sum-1.0) > epsilon )
837 PUT_ERR("ERROR in sum of corrected weights\n");
838
839 for (size_t i = 0; i < nSourceCells; ++i)
840 yac_free_grid_cell(&SourceCells[i]);
841 yac_free_grid_cell(&TargetCell);
842 free(SourceCells);
843}
844
845static void utest_check_partial_areas_deg(
846 double * tgt_lon, double * tgt_lat,
847 enum yac_edge_type * tgt_edges, int tgt_count,
848 double ** src_lons, double ** src_lats,
849 enum yac_edge_type ** src_edges, int * src_counts, size_t nSourceCells) {
850
851 utest_check_partial_areas(
852 tgt_lon, tgt_lat, tgt_edges, tgt_count,
853 src_lons, src_lats, src_edges, src_counts, nSourceCells,
854 generate_cell_deg);
855}
856
857static void utest_check_partial_areas_rad(
858 double * tgt_lon, double * tgt_lat,
859 enum yac_edge_type * tgt_edges, int tgt_count,
860 double ** src_lons, double ** src_lats,
861 enum yac_edge_type ** src_edges, int * src_counts, size_t nSourceCells) {
862
863 utest_check_partial_areas(
864 tgt_lon, tgt_lat, tgt_edges, tgt_count,
865 src_lons, src_lats, src_edges, src_counts, nSourceCells,
866 generate_cell_rad);
867}
yac_huiliers_area
double yac_huiliers_area(struct yac_grid_cell cell)
Area calculation on a unit sphere taken from ESMF based on L'Huilier's Theorem.
Definition area.c:449
yac_huiliers_area_info
double yac_huiliers_area_info(struct yac_grid_cell cell, double *barycenter, double sign)
Definition area.c:539
area.h
Structs and interfaces for area calculations.
basic_grid.h
yac_correct_weights
void yac_correct_weights(size_t nSourceCells, double *weight)
correct interpolation weights
Definition clipping.c:1445
yac_compute_overlap_info
void yac_compute_overlap_info(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, double *overlap_areas, double(*overlap_barycenters)[3])
calculates partial areas for all overlapping parts of the source cells with arbitrary target cells,...
Definition clipping.c:119
yac_compute_overlap_areas
void yac_compute_overlap_areas(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, double *partial_areas)
calculates partial areas for all overlapping parts of the source cells with arbitrary target cells,...
Definition clipping.c:278
generate_cell
static void generate_cell(struct point_list *list, struct yac_grid_cell *cell)
Definition clipping.c:1758
yac_compute_overlap_buf_free
void yac_compute_overlap_buf_free()
Definition clipping.c:88
YAC_RAD
#define YAC_RAD
Definition dynamic_config.c:12
normalise_vector
static void normalise_vector(double v[])
Definition geometry.h:635
get_vector_angle
static double get_vector_angle(double const a[3], double const b[3])
Definition geometry.h:340
yac_free_grid_cell
void yac_free_grid_cell(struct yac_grid_cell *cell)
Definition grid_cell.c:44
yac_edge_type
yac_edge_type
Definition grid_cell.h:12
YAC_GREAT_CIRCLE_EDGE
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
YAC_LAT_CIRCLE_EDGE
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
YAC_LON_CIRCLE_EDGE
@ YAC_LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:15
test_core.weights
weights
Definition test_core.py:63
test_raw_exchange.i
int i
Definition test_raw_exchange.py:89
XYZtoLL
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition perf_toy_icon.c:355
xmalloc
#define xmalloc(size)
Definition ppm_xfuncs.h:66
yac_grid_cell
Definition grid_cell.h:18
generate_cell_rad
struct yac_grid_cell generate_cell_rad(double *lon, double *lat, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:42
generate_cell_deg
struct yac_grid_cell generate_cell_deg(double *lon, double *lat, enum yac_edge_type *edge_type, size_t num_corners)
Definition test_common.c:36
intersect
int intersect(enum yac_edge_type edge_type_a, double lon_a, double lat_a, double lon_b, double lat_b, enum yac_edge_type edge_type_b, double lon_c, double lat_c, double lon_d, double lat_d, double *intersection)
Definition test_common.c:64
test_common.h
gc_edges
static enum yac_edge_type gc_edges[]
Definition test_compute_overlap_area.c:16
latlon_edges
static enum yac_edge_type latlon_edges[4]
Definition test_compute_overlap_area.c:26
TEST_EXIT_CODE
#define TEST_EXIT_CODE
Definition tests.h:14
PUT_ERR
#define PUT_ERR(string)
Definition tests.h:10
LLtoXYZ
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition toy_scrip.c:587
MAX
#define MAX(a, b)
Definition utils_common.h:66