YAC 3.12.0
Yet Another Coupler
Loading...
Searching...
No Matches
test_compute_overlap_area.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 "tests.h"
6#include "clipping.h"
7#include "geometry.h"
8#include "area.h"
9#include "test_common.h"
10
31
32static void utest_check_overlap(
33 double * lon_a, double * lat_a, enum yac_edge_type * edges_a, int count_a,
34 double * lon_b, double * lat_b, enum yac_edge_type * edges_b, int count_b,
35 double ref_area, double const * ref_barycenter);
36static void utest_check_overlap_polygons(
37 struct yac_grid_cell * Polygons, double ref_area,
38 double const * ref_barycenter);
39
40int main (void) {
41
42 {
43 /* intersection between concave pentagon with a triangle
44
45 10.0 (lon) 11.0 (lon)
46 32.0 (lat) 32.0 (lat)
47 x--------x
48 / \
49 / 10.5 \
50 / 31.0 \
51 9.5 (lon) / x \ 11.5
52 30.0 (lat) x x 30.0
53
54
55 with
56
57
58 10.5 (lon)
59 31.5 (lat)
60 x
61 / \
62 / \
63 10.0 / \ 11.0
64 29.5 x-------x 29.5
65 */
66
67 // the reference area any barycenter
68 double ref_area;
69 double ref_barycenter[3] = {0, 0, 0};
70 {
71 double intersection[2][3], intersection_lon[2], intersection_lat[2];
72
73 if (!intersect(YAC_GREAT_CIRCLE_EDGE, 11.5, 30.0, 10.5, 31.0,
74 YAC_GREAT_CIRCLE_EDGE, 11.0, 29.5, 10.5, 31.5,
75 intersection[0]) ||
76 !intersect(YAC_GREAT_CIRCLE_EDGE, 9.5, 30.0, 10.5, 31.0,
77 YAC_GREAT_CIRCLE_EDGE, 10.0, 29.5, 10.5, 31.5,
78 intersection[1]))
79 return EXIT_FAILURE;
80 for (int i = 0; i < 2; ++i) {
81 XYZtoLL(intersection[i], &intersection_lon[i], &intersection_lat[i]);
82 intersection_lon[i] /= YAC_RAD;
83 intersection_lat[i] /= YAC_RAD;
84 }
85
86 struct yac_grid_cell Polygon =
87 generate_cell_deg(
88 (double[]){10.5, intersection_lon[0], 10.5, intersection_lon[1]},
89 (double[]){31.0, intersection_lat[0], 31.5, intersection_lat[1]},
90 gc_edges, 4);
91 ref_area = yac_huiliers_area_info(Polygon, ref_barycenter, 1.0);
92 normalise_vector(ref_barycenter);
93 yac_free_grid_cell(&Polygon);
94 }
95
96 utest_check_overlap(
97 (double[]){10.0, 11.0, 11.5, 10.5, 9.5},
98 (double[]){32.0, 32.0, 30.0, 31.0, 30.0}, gc_edges, 5,
99 (double[]){10.5, 11.0, 10.0},
100 (double[]){31.5, 29.5, 29.5}, gc_edges, 3,
101 ref_area, ref_barycenter);
102 }
103
104 {
105 /* intersection between square and a triangle
106
107 -5.0 5.0
108 x--------x 5.0
109 | |
110 | |
111 | |
112 x--------x -5.0
113
114 with
115
116 -4.0 0.0 4.0
117 x 6.0
118 / \
119 / \
120 / \
121 x-------x -4.0
122 */
123 // the reference area any barycenter
124 double ref_area;
125 double ref_barycenter[3] = {0, 0, 0};
126 {
127 double intersection[2][3], intersection_lon[2], intersection_lat[2];
128
129 if (!intersect(YAC_GREAT_CIRCLE_EDGE, 0, 6, 4, -4,
130 YAC_GREAT_CIRCLE_EDGE, -5, 5, 5, 5,
131 intersection[0]) ||
132 !intersect(YAC_GREAT_CIRCLE_EDGE, 0, 6, -4, -4,
133 YAC_GREAT_CIRCLE_EDGE, -5, 5, 5, 5,
134 intersection[1]))
135 return EXIT_FAILURE;
136 for (int i = 0; i < 2; ++i) {
137 XYZtoLL(intersection[i], &intersection_lon[i], &intersection_lat[i]);
138 intersection_lon[i] /= YAC_RAD;
139 intersection_lat[i] /= YAC_RAD;
140 }
141
142 struct yac_grid_cell Polygon =
143 generate_cell_deg(
144 (double[]){-4, 4, intersection_lon[0], intersection_lon[1]},
145 (double[]){-4, -4, intersection_lat[0], intersection_lat[1]}, gc_edges, 4);
146 ref_area = yac_huiliers_area_info(Polygon, ref_barycenter, 1.0);
147 normalise_vector(ref_barycenter);
148 yac_free_grid_cell(&Polygon);
149 }
150
151 utest_check_overlap(
152 (double[]){-5, 5, 5, -5}, (double[]){-5, -5, 5, 5}, gc_edges, 4,
153 (double[]){-4, 4, 0}, (double[]){-4, -4, 6}, gc_edges, 3,
154 ref_area, ref_barycenter);
155 }
156
157 {
158 /* set the test cells (the target node is one of the corners of the
159 concave polygon; the center of both shapes
160 is the north pole)
161
162
163 +------+
164 | |
165 | /\ |
166 | / \ |
167 +/ \+
168
169 with
170
171 +------+
172 | |
173 | |
174 | |
175 +------+
176
177 */
178
179 // the reference area and barycenter
180 double ref_area = 0.0;
181 double ref_barycenter[] = {0.0, 0.0, 0.0};
182 for (int i = 0; i < 3; ++i) {
183 double tri_lon[3][3] = {{0,0,90},{90,90,180},{180,180,270}};
184 double tri_lat[3] = {90,89,89};
185 struct yac_grid_cell Tri =
186 generate_cell_deg(tri_lon[i], tri_lat, gc_edges, 3);
187 ref_area += yac_huiliers_area_info(Tri, ref_barycenter, 1.0);
188 yac_free_grid_cell(&Tri);
189 }
190 normalise_vector(ref_barycenter);
191
192 utest_check_overlap(
193 (double[]){0,0,90,180,270}, (double[]){90,88,88,88,88}, gc_edges, 5,
194 (double[]){0,90,180,270}, (double[]){89,89,89,89}, gc_edges, 4,
195 ref_area, ref_barycenter);
196 }
197
198 {
199 /*
200
201 +---+ +---+
202 | | | |
203 | | +--+ | |
204 | | | | | |
205 | | +--+ | |
206 | | | |
207 | +--------+ |
208 | |
209 +----------------+
210
211 */
212
213 // the reference area and barycenter
214 double ref_area = 0.0;
215 double * ref_barycenter = NULL;
216
217 utest_check_overlap(
218 (double[]){2,3,3,-3,-3,-2,-2,2},
219 (double[]){3,3,-3,-3,3,3,-2,-2}, gc_edges, 8,
220 (double[]){1,1,-1,-1},
221 (double[]){1,-1,-1,1}, gc_edges, 4,
222 ref_area, ref_barycenter);
223 }
224
225 {
226 /*
227
228 +------+------+
229 |......|......|
230 |..+---+---+..|
231 |..| |..|
232 |..| +---+ |..|
233 |..| |...| |..|
234 |..| +---+ |..|
235 |..| |..|
236 |..+-------+..|
237 |.............|
238 +-------------+
239
240 */
241
242 // the reference area and barycenter
243 double ref_area = 0.0;
244 double * ref_barycenter = NULL;
245
246 utest_check_overlap(
247 (double[]){0,3,3,-3,-3,0,0,-2,-2,2,2,0},
248 (double[]){3,3,-3,-3,3,3,2,2,-2,-2,2,2}, gc_edges, 12,
249 (double[]){1,1,-1,-1}, (double[]){1,-1,-1,1}, gc_edges, 4,
250 ref_area, ref_barycenter);
251 }
252
253 {
254 /*
255
256 + +
257 |\ /|
258 |..\ /..|
259 |....\ /....|
260 |.....+ +.....|
261 |.....| |.....|
262 |.....| + |.....|
263 |.....| / \ |.....|
264 |...../-------\.....|
265 |.../...........\...|
266 |..+-------------+..|
267 |...................|
268 +-------------------+
269
270 */
271
272 // the reference area and barycenter
273 double ref_area;
274 double ref_barycenter[3] = {0.0, 0.0, 0.0};
275 {
276 struct yac_grid_cell Tri_a =
277 generate_cell_deg(
278 (double[]){90,180,90},
279 (double[]){90,87.5,87.5}, gc_edges, 3);
280 struct yac_grid_cell Tri_b =
281 generate_cell_deg(
282 (double[]){90,180,90},
283 (double[]){90,88,88}, gc_edges, 3);
284 ref_area = yac_huiliers_area_info(Tri_a, ref_barycenter, 1.0) +
285 yac_huiliers_area_info(Tri_b, ref_barycenter, -1.0);
286 if (ref_area < 0.0) {
287 ref_area = -ref_area;
288 ref_barycenter[0] = -ref_barycenter[0];
289 ref_barycenter[1] = -ref_barycenter[1];
290 ref_barycenter[2] = -ref_barycenter[2];
291 }
292 normalise_vector(ref_barycenter);
293 yac_free_grid_cell(&Tri_b);
294 yac_free_grid_cell(&Tri_a);
295 }
296
297 utest_check_overlap(
298 (double[]){0,90,180,270,270,180,90,0},
299 (double[]){88,88,88,88,87,87,87,87}, gc_edges, 8,
300 (double[]){90,180,90},
301 (double[]){90,87.5,87.5}, gc_edges, 3,
302 ref_area, ref_barycenter);
303 }
304
305 {
306 /*
307
308 +---------+---------+
309 |.........|.........|
310 |.........|.........|
311 |.........|.........|
312 |.....+---+---+.....|
313 |.....| |.....|
314 |.....| + |.....|
315 |.....| / \ |.....|
316 |...../-------\.....|
317 |.../...........\...|
318 |..+-------------+..|
319 |...................|
320 +-------------------+
321
322 */
323
324 // the reference area and barycenter
325 double ref_area = 0.0;
326 double ref_barycenter[3] = {0.0, 0.0, 0.0};
327 {
328 struct yac_grid_cell Quad =
329 generate_cell_deg(
330 (double[]){90,180,180,90},
331 (double[]){87.5,87.5,89,89}, gc_edges, 4);
332 ref_area = yac_huiliers_area_info(Quad, ref_barycenter, 1.0);
333 normalise_vector(ref_barycenter);
334 yac_free_grid_cell(&Quad);
335 }
336
337 utest_check_overlap(
338 (double[]){315,270,180,90,0,315,
339 315,0,90,180,270,315},
340 (double[]){89,89,89,89,89,89,
341 86,86,86,86,86,86}, gc_edges, 12,
342 (double[]){90,90,180},
343 (double[]){90,87.5,87.5}, gc_edges, 3,
344 ref_area, ref_barycenter);
345 }
346
347 {
348 /*
349
350 +---------+---------+
351 |.........|.........|
352 |.........|......+..|
353 |.........|...../|..|
354 |.....+---+---/..|..|
355 |.....| / |..|..|
356 +.....+ / +..|..+
357 |.....| / |..|..|
358 |...../---+---+..|..|
359 |.../............|..|
360 |..+-------------+..|
361 |...................|
362 +---------+---------+
363
364 */
365
366 // the reference area and barycenter
367 for (int i = 0; i < 4; ++i) {
368
369 double tri_rotation = (double)i * 90.0;
370
371 double ref_area = 0.0;
372 double ref_barycenter[3] = {0.0, 0.0, 0.0};
373 {
374 struct yac_grid_cell Poly =
375 generate_cell_deg(
376 (double[]){0+tri_rotation,90+tri_rotation,180+tri_rotation,
377 180+tri_rotation,135+tri_rotation,90+tri_rotation,
378 45+tri_rotation,0+tri_rotation},
379 (double[]){87.5,87.5,87.5,89,89,89,89,89}, gc_edges, 8);
380 ref_area = yac_huiliers_area_info(Poly, ref_barycenter, 1.0);
381 normalise_vector(ref_barycenter);
382 yac_free_grid_cell(&Poly);
383 }
384
385 utest_check_overlap(
386 (double[]){315,270,225,180,135,90,45,0,315,
387 315,0,45,90,135,180,225,270,315},
388 (double[]){89,89,89,89,89,89,89,89,89,
389 86,86,86,86,86,86,86,86,86}, gc_edges, 18,
390 (double[]){0+tri_rotation,90+tri_rotation,180+tri_rotation},
391 (double[]){87.5,87.5,87.5}, gc_edges, 3,
392 ref_area, ref_barycenter);
393 }
394 }
395
396 {
397 /*
398
399 +---------+---------+
400 |.........|.........|
401 |..+------|------+..|
402 |..|......|......|..|
403 |..|..+---+---+..|..|
404 |..|..| |..|..|
405 +..|..+ +..|..+
406 |..|..| |..|..|
407 |..|..+---+---+..|..|
408 |..|.............|..|
409 |..+-------------+..|
410 |...................|
411 +---------+---------+
412
413 */
414
415 // the reference area and barycenter
416 double ref_area = 0.0;
417 double ref_barycenter[3];
418 {
419 struct yac_grid_cell Square =
420 generate_cell_deg(
421 (double[]){0,90,180,270},
422 (double[]){87.5,87.5,87.5,87.5}, gc_edges, 4);
423 struct yac_grid_cell Polygon =
424 generate_cell_deg(
425 (double[]){0,45,90,135,180,225,270,315},
426 (double[]){89,89,89,89,89,89,89,89}, gc_edges, 8);
427 ref_area = yac_huiliers_area(Square) - yac_huiliers_area(Polygon);
428 yac_free_grid_cell(&Polygon);
429 yac_free_grid_cell(&Square);
430 LLtoXYZ_deg(90, 90, ref_barycenter);
431 }
432
433 utest_check_overlap(
434 (double[]){315,270,225,180,135,90,45,0,315,
435 315,0,45,90,135,180,225,270,315},
436 (double[]){89,89,89,89,89,89,89,89,89,
437 86,86,86,86,86,86,86,86,86}, gc_edges, 18,
438 (double[]){0,90,180,270},
439 (double[]){87.5,87.5,87.5,87.5}, gc_edges, 4,
440 ref_area, ref_barycenter);
441 }
442
443 {
444 /*
445
446 +---------+---------+
447 |.........|.........|
448 |.........|.........|
449 |.........|.........|
450 |.....+---+---+.....|
451 |.....| |.....|
452 +.....+ +.....+
453 |.....| |.....|
454 |.....+---+---+.....|
455 |...................|
456 |...................|
457 |...................|
458 +---------+---------+
459
460 +------+------+
461 |.............|
462 |.............|
463 |.....+++.....|
464 +.....+ +.....+
465 |.....+++.....|
466 |......|......|
467 |......|......|
468 +------+------+
469
470 */
471
472 // the reference area and barycenter
473 double ref_area;
474 double ref_barycenter[3];
475 {
476 struct yac_grid_cell Polygon =
477 generate_cell_deg(
478 (double[]){315,270,225,180,135,90,45,0,315,
479 315,0,45,90,135,180,225,270,315},
480 (double[]){88,88,88,88,88,88,88,88,88,
481 87,87,87,87,87,87,87,87,87}, gc_edges, 18);
482 ref_area = yac_huiliers_area(Polygon);
483 yac_free_grid_cell(&Polygon);
484 LLtoXYZ_deg(90, 90, ref_barycenter);
485 }
486
487 utest_check_overlap(
488 (double[]){135,90,45,0,315,270,225,180,135,
489 135,180,225,270,315,0,45,90,135},
490 (double[]){88,88,88,88,88,88,88,88,88,
491 86,86,86,86,86,86,86,86,86}, gc_edges, 18,
492 (double[]){315,270,225,180,135,90,45,0,315,
493 315,0,45,90,135,180,225,270,315},
494 (double[]){89,89,89,89,89,89,89,89,89,
495 87,87,87,87,87,87,87,87,87}, gc_edges, 18,
496 ref_area, ref_barycenter);
497 }
498
499
500
501 {
502 /*
503
504
505 + 5
506 / \
507 / \
508 / \
509 +-------------+ 0
510 | |
511 | |
512 | |
513 +-------------+ -5
514 -5 5
515
516 */
517
518 // the reference area and barycenter
519 double ref_area = 0.0;
520 double ref_barycenter[3] = {0.0, 0.0, 0.0};
521
522 utest_check_overlap(
523 (double[]){-5.0,5.0,5.0,-5.0},
524 (double[]){-5.0,-5.0,0.0,0.0}, latlon_edges, 4,
525 (double[]){-5.0,5.0,0.0},
526 (double[]){0.0,0.0,5.0}, gc_edges, 3,
527 ref_area, ref_barycenter);
528 }
529
530 { // test reproducing an issue found in toy_multi example
531 struct yac_grid_cell cells[2] =
532 {{.coordinates_xyz =
533 (double[3][3])
534 {{-0.11609962857526578,
535 0.1297445660881493,
536 -0.98472697932740894},
537 {-0.1873171704582404,
538 0.13116089479097318,
539 -0.97350351685504954},
540 {-0.16162842294547278,
541 0.064722709158966607,
542 -0.98472697932740894}},
543 .edge_type =
544 (enum yac_edge_type[3])
545 {YAC_GREAT_CIRCLE_EDGE,
546 YAC_GREAT_CIRCLE_EDGE,
547 YAC_GREAT_CIRCLE_EDGE},
548 .num_corners = 3,
549 .array_size = 3},
550 {.coordinates_xyz =
551 (double[4][3])
552 {{-0.10626911810438046,
553 0.10117661128547845,
554 -0.98917650996478101},
555 {-0.10872011017156734,
556 0.098538164069449888,
557 -0.98917650996478101},
558 {-0.12667440386975459,
559 0.11481098733454023,
560 -0.98527764238894122},
561 {-0.12381864923052141,
562 0.11788515390505606,
563 -0.98527764238894122}},
564 .edge_type = latlon_edges,
565 .num_corners = 4,
566 .array_size = 4}};
567 struct yac_grid_cell overlap_buffer =
568 {.coordinates_xyz = NULL,
569 .edge_type = NULL,
570 .num_corners = 0,
571 .array_size = 0};
572 yac_cell_clipping(1, &cells[0], cells[1], &overlap_buffer);
573 double ref_overlap_areas = yac_huiliers_area(overlap_buffer);
574 double ref_overlap_barycenters[3] = {0.0, 0.0, 0.0};
575 utest_check_overlap_polygons(cells, ref_overlap_areas, ref_overlap_barycenters);
576 yac_free_grid_cell(&overlap_buffer);
577 }
578
579 { // test reproducing an issue found in toy_multi example
580 struct yac_grid_cell cells[2] =
581 {{.coordinates_xyz =
582 (double[4][3])
583 {{0.12226322617998441,
584 0.0060064071448744901,
585 0.99247953459870997},
586 {0.12207899817177321,
587 0.0090050879009709785,
588 0.99247953459870997},
589 {0.097751558641606923,
590 0.0072105881536315003,
591 0.99518472667219682},
592 {0.097899074391390048,
593 0.0048094731201957178,
594 0.99518472667219682}},
595 .edge_type = latlon_edges,
596 .num_corners = 4,
597 .array_size = 4},
598 {.coordinates_xyz =
599 (double[4][3])
600 {{ 0.13173419827179747,
601 -0.0024316641355669943,
602 0.99128209305687476},
603 { 0.12602575200455376,
604 0.0081957056315718202,
605 0.99199311501687715},
606 { 0.11521731465958367,
607 0.0019341931455107773,
608 0.99333842636812875},
609 { 0.12147807203778695,
610 -0.0083897691028181811,
611 0.99255865810962707}},
612 .edge_type =
613 (enum yac_edge_type[4])
614 {YAC_GREAT_CIRCLE_EDGE,
615 YAC_GREAT_CIRCLE_EDGE,
616 YAC_GREAT_CIRCLE_EDGE,
617 YAC_GREAT_CIRCLE_EDGE},
618 .num_corners = 4,
619 .array_size = 4}};
620 struct yac_grid_cell overlap_buffer =
621 {.coordinates_xyz = NULL,
622 .edge_type = NULL,
623 .num_corners = 0,
624 .array_size = 0};
625 yac_cell_clipping(1, &cells[0], cells[1], &overlap_buffer);
626 double ref_overlap_areas = yac_huiliers_area(overlap_buffer);
627 double ref_overlap_barycenters[3] = {0.0, 0.0, 0.0};
628 yac_free_grid_cell(&overlap_buffer);
629 utest_check_overlap_polygons(cells, ref_overlap_areas, ref_overlap_barycenters);
630 }
631
632 { // test reproducing an issue found in toy_multi example
633 struct yac_grid_cell cells[2] =
634 {{.coordinates_xyz =
635 (double[4][3])
636 {{-0.67635734277981696,
637 -0.34044967630401168,
638 0.65317284295377664},
639 {-0.66779858328675834,
640 -0.35694577933893479,
641 0.65317284295377664},
642 {-0.65346055329338182,
643 -0.34928194263987855,
644 0.67155895484701855},
645 {-0.66183555116521386,
646 -0.33314002067992043,
647 0.67155895484701855}},
648 .edge_type = latlon_edges,
649 .num_corners = 4,
650 .array_size = 4},
651 {.coordinates_xyz =
652 (double[4][3])
653 {{-0.659061024045835,
654 -0.36232186404812883,
655 0.65906102404583489},
656 {-0.66524659712730871,
657 -0.33896007142593126,
658 0.66524659712730871},
659 {-0.64171186251713952,
660 -0.34818611452313608,
661 0.68335372623412638},
662 {-0.63542071129861122,
663 -0.37199386619312252,
664 0.67665433063526625}},
665 .edge_type =
666 (enum yac_edge_type[4])
667 {YAC_GREAT_CIRCLE_EDGE,
668 YAC_GREAT_CIRCLE_EDGE,
669 YAC_GREAT_CIRCLE_EDGE,
670 YAC_GREAT_CIRCLE_EDGE},
671 .num_corners = 4,
672 .array_size = 4}};
673 struct yac_grid_cell overlap_buffer =
674 {.coordinates_xyz = NULL,
675 .edge_type = NULL,
676 .num_corners = 0,
677 .array_size = 0};
678 yac_cell_clipping(1, &cells[0], cells[1], &overlap_buffer);
679 double ref_overlap_areas = yac_huiliers_area(overlap_buffer);
680 double ref_overlap_barycenters[3] = {0.0, 0.0, 0.0};
681 yac_free_grid_cell(&overlap_buffer);
682 utest_check_overlap_polygons(cells, ref_overlap_areas, ref_overlap_barycenters);
683 }
684
685 { // testing issue found with FESOM grid
686
687 /*
688
689 +
690 / \
691 / \
692 + +
693 - -
694 + +
695 - -
696 + +
697 \ /
698 \ /
699 +
700 */
701
702 // the reference area and barycenter
703 double ref_area = 0.0;
704 double ref_barycenter[3];
705 {
706 struct yac_grid_cell Polygon =
707 generate_cell_deg(
708 (double[]){0.0,0.5,2.0,0.5,0.0,-0.5,-2.0,-0.5},
709 (double[]){2.0,0.5,0.0,-0.5,-2.0,-0.5,0.0,0.5}, gc_edges, 8);
710 ref_area = yac_huiliers_area(Polygon);
711 yac_free_grid_cell(&Polygon);
712 LLtoXYZ_deg(0, 0, ref_barycenter);
713 }
714 utest_check_overlap(
715 (double[]){0.0,0.5,2.0,0.5,0.0,-0.5,-2.0,-0.5},
716 (double[]){2.0,0.5,0.0,-0.5,-2.0,-0.5,0.0,0.5}, gc_edges, 8,
717 (double[]){-3.0,3.0,3.0,-3.0},
718 (double[]){3.0,3.0,-3.0,-3.0}, gc_edges, 4,
719 ref_area, ref_barycenter);
720 }
721
722 { // testing issue found with FESOM grid
723
724 /*
725
726 +
727 / \
728 / \
729 + +
730 - -
731 + +
732 - -
733 + +
734 \ /
735 \ /
736 +
737 */
738
739 // the reference area and barycenter
740 double ref_area = 0.0;
741 double ref_barycenter[3];
742 {
743 struct yac_grid_cell Polygon =
744 generate_cell_deg(
745 (double[]){-0.5,0.5,0.5,-0.5},
746 (double[]){0.5,0.5,-0.5,-0.5}, gc_edges, 4);
747 ref_area = yac_huiliers_area(Polygon);
748 yac_free_grid_cell(&Polygon);
749 LLtoXYZ_deg(0, 0, ref_barycenter);
750 }
751 utest_check_overlap(
752 (double[]){0.0,0.5,2.0,0.5,0.0,-0.5,-2.0,-0.5},
753 (double[]){2.0,0.5,0.0,-0.5,-2.0,-0.5,0.0,0.5}, gc_edges, 8,
754 (double[]){-0.5,0.5,0.5,-0.5},
755 (double[]){0.5,0.5,-0.5,-0.5}, gc_edges, 4,
756 ref_area, ref_barycenter);
757 }
758
759 {
760 struct yac_grid_cell src_cell =
761 {.coordinates_xyz =
762 (double[4][3])
763 {{5.06731853971386536628e-01,
764 5.06731853971386536628e-01,
765 6.97456562333055085645e-01},
766 {4.56966311668627278575e-01,
767 6.28960169645094047119e-01,
768 6.28960169645094047119e-01},
769 {5.77350269189625842081e-01,
770 5.77350269189625731059e-01,
771 5.77350269189625842081e-01},
772 {6.28960169645094047119e-01,
773 4.56966311668627389597e-01,
774 6.28960169645094047119e-01}},
775 .edge_type =
776 (enum yac_edge_type[4])
777 {YAC_GREAT_CIRCLE_EDGE,
778 YAC_GREAT_CIRCLE_EDGE,
779 YAC_GREAT_CIRCLE_EDGE,
780 YAC_GREAT_CIRCLE_EDGE},
781 .num_corners = 4,
782 .array_size = 4};
783 struct yac_grid_cell tgt_cell =
784 {.coordinates_xyz =
785 (double[3][3])
786 {{5.26518818746262273756e-01,
787 4.92120108892697694092e-01,
788 6.93250122199397744716e-01},
789 {5.28392100138220688343e-01,
790 4.99031313570251600087e-01,
791 6.86854814781020395209e-01},
792 {5.20491669253990818511e-01,
793 4.99252122636297701597e-01,
794 6.92701768642426274347e-01}},
795 .edge_type =
796 (enum yac_edge_type[3])
797 {YAC_GREAT_CIRCLE_EDGE,
798 YAC_GREAT_CIRCLE_EDGE,
799 YAC_GREAT_CIRCLE_EDGE},
800 .num_corners = 3,
801 .array_size = 3};
802
803 double areas[2];
804 double barycenter[3];
805 yac_compute_overlap_info(
806 1, &src_cell, tgt_cell, &areas[0], NULL);
807 yac_compute_overlap_info(
808 1, &src_cell, tgt_cell, &areas[1], &barycenter);
809
810 if (fabs(areas[0] - areas[1]) > (areas[0] * 1e-6))
811 PUT_ERR("error in yac_compute_overlap_info");
812 }
813
814 { // test an issue found in OYAC
815 // (when the intersection between two cells is a polygon with two corners,
816 // where one edge is gc while the other is lat, an assert was raised)
817
818 /* pole triangle with lat edge
819 x 90
820 / \
821 / \
822 / \
823 x x 85
824 \_____/
825
826 square with gc edges very close to the triangle
827
828 x--------x 84.99
829 | |
830 | |
831 | |
832 x--------x
833 */
834
835 double lat_edge_middle[3];
836 double gc_edge_middle[3];
837 double overlap_intersections[2][3];
838
839 LLtoXYZ_deg(0.0, 85.0, lat_edge_middle);
840 if (!intersect(YAC_GREAT_CIRCLE_EDGE, -5.0, 84.99, 5.0, 84.99,
841 YAC_LON_CIRCLE_EDGE, 0.0, 90.0, 0.0, 80.0,
842 gc_edge_middle) ||
843 !intersect(YAC_GREAT_CIRCLE_EDGE, -5.0, 84.99, 5.0, 84.99,
844 YAC_LAT_CIRCLE_EDGE, -10.0, 85.0, 0.0, 85.0,
845 overlap_intersections[0]) ||
846 !intersect(YAC_GREAT_CIRCLE_EDGE, -5.0, 84.99, 5.0, 84.99,
847 YAC_LAT_CIRCLE_EDGE, 0.0, 85.0, 10.0, 85.0,
848 overlap_intersections[1]))
849 return EXIT_FAILURE;
850
851 struct yac_grid_cell ref_cell_overlap =
852 {.coordinates_xyz = overlap_intersections,
853 .edge_type =
854 (enum yac_edge_type[2]){YAC_GREAT_CIRCLE_EDGE, YAC_LAT_CIRCLE_EDGE},
855 .num_corners = 2,
856 .array_size = 2};
857 double ref_area = yac_huiliers_area(ref_cell_overlap);
858
859 // rough estimate of barycenter of the intersection between the two cells
860 double ref_barycenter[3] = {lat_edge_middle[0] + gc_edge_middle[0],
861 lat_edge_middle[1] + gc_edge_middle[1],
862 lat_edge_middle[2] + gc_edge_middle[2]};
863 normalise_vector(ref_barycenter);
864
865 utest_check_overlap(
866 (double[]){-5.0,5.0,5.0,-5.0},
867 (double[]){90.0,90.0,85.0,85.0}, latlon_edges, 4,
868 (double[]){-5.0,5.0,5.0,-5.0},
869 (double[]){84.99,84.99,80.0,80.0}, gc_edges, 4,
870 ref_area, ref_barycenter);
871 }
872
873 // cleanup
874 yac_compute_overlap_buf_free();
875
876 return TEST_EXIT_CODE;
877}
878
879static void utest_check_overlap_polygons(
880 struct yac_grid_cell * Polygons, double ref_area,
881 double const * ref_barycenter) {
882
883 double area_tolerance = MAX(ref_area * 1e-6, 1e-10);
884
885 for (int test_order = 0; test_order < 2; ++test_order) {
886 double area;
887 if (ref_barycenter) {
888 double barycenter[3] = {0.0,0.0,0.0};
889 yac_compute_overlap_info(
890 1, &Polygons[test_order], Polygons[test_order^1],
891 &area, &barycenter);
892 if (fabs(ref_area-area) > area_tolerance)
893 PUT_ERR("error in yac_compute_overlap_info "
894 "area difference is too large.\n");
895 if ((barycenter[0] != barycenter[0]) ||
896 (barycenter[1] != barycenter[1]) ||
897 (barycenter[2] != barycenter[2]))
898 PUT_ERR("error in yac_compute_overlap_info "
899 "barycenter contains nan's.\n");
900 if ((ref_barycenter[0] != 0.0) ||
901 (ref_barycenter[1] != 0.0) ||
902 (ref_barycenter[2] != 0.0)) {
903 double angle_diff = get_vector_angle(ref_barycenter, barycenter);
904 if (fabs(angle_diff) > 1e-9)
905 PUT_ERR("error in yac_compute_overlap_info "
906 "barycenter difference is too large.\n");
907 }
908 }
909 yac_compute_overlap_areas(
910 1, &Polygons[test_order], Polygons[test_order^1], &area);
911 if (fabs(ref_area-area) > area_tolerance)
912 PUT_ERR("error in yac_compute_overlap_areas "
913 "area difference is too large.\n");
914 }
915}
916
917static void utest_check_overlap(
918 double * lon_a, double * lat_a, enum yac_edge_type * edges_a, int count_a,
919 double * lon_b, double * lat_b, enum yac_edge_type * edges_b, int count_b,
920 double ref_area, double const * ref_barycenter) {
921
922 for (int order_a = -1; order_a <=1; order_a += 2) {
923 for (int order_b = -1; order_b <=1; order_b += 2) {
924 for (int start_idx_a = 0; start_idx_a < count_a; ++start_idx_a) {
925 for (int start_idx_b = 0; start_idx_b < count_b; ++start_idx_b) {
926
927 double curr_lon_a[count_a], curr_lat_a[count_a];
928 double curr_lon_b[count_b], curr_lat_b[count_b];
929 enum yac_edge_type curr_edges_a[count_a];
930 enum yac_edge_type curr_edges_b[count_b];
931
932 for (int i = 0; i < count_a; ++i) {
933 int idx = (start_idx_a + order_a * i + count_a) % count_a;
934 curr_lon_a[i] = lon_a[idx];
935 curr_lat_a[i] = lat_a[idx];
936 idx = (idx - (order_a < 0) + count_a) % count_a;
937 curr_edges_a[i] = edges_a[idx];
938 }
939 for (int i = 0; i < count_b; ++i) {
940 int idx = (start_idx_b + order_b * i + count_b) % count_b;
941 curr_lon_b[i] = lon_b[idx];
942 curr_lat_b[i] = lat_b[idx];
943 idx = (idx - (order_b < 0) + count_b) % count_b;
944 curr_edges_b[i] = edges_b[idx];
945 }
946
947 struct yac_grid_cell Polygons[2] = {
948 generate_cell_deg(curr_lon_a, curr_lat_a, curr_edges_a, count_a),
949 generate_cell_deg(curr_lon_b, curr_lat_b, curr_edges_b, count_b)};
950
951 utest_check_overlap_polygons(Polygons, ref_area, ref_barycenter);
952
953 for (int i = 0; i < 2; ++i) yac_free_grid_cell(&Polygons[i]);
954 }
955 }
956 }
957 }
958}
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.
yac_cell_clipping
void yac_cell_clipping(size_t N, struct yac_grid_cell *source_cell, struct yac_grid_cell target_cell, struct yac_grid_cell *overlap_buffer)
cell clipping to get the cells describing the intersections
Definition clipping.c:1088
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
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
LLtoXYZ_deg
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:269
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_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
yac_grid_cell
Definition grid_cell.h:18
yac_grid_cell::coordinates_xyz
double(* coordinates_xyz)[3]
Definition grid_cell.h:19
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
MAX
#define MAX(a, b)
Definition utils_common.h:66