YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
interp_method_conserv.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
14#include "clipping.h"
15#include "area.h"
16#include "float.h"
17
18#define AREA_TOL_FACTOR (1e-6)
19
20static size_t do_search_conserv_1st_order(struct interp_method * method,
21 struct yac_interp_grid * interp_grid,
22 size_t * tgt_points, size_t count,
23 struct yac_interp_weights * weights);
24static size_t do_search_conserv_2nd_order(struct interp_method * method,
25 struct yac_interp_grid * interp_grid,
26 size_t * tgt_points, size_t count,
27 struct yac_interp_weights * weights);
28static void delete_conserv(struct interp_method * method);
29
30static struct interp_method_vtable
34
35static struct interp_method_vtable
39
47
53
59
62 size_t n;
63};
64
66 struct {
67 size_t local_id;
69 } src, tgt;
70 double norm_area;
71 double area;
72 double barycenter[3];
74};
75
77 struct yac_const_basic_grid_data * basic_grid_data) {
78
79 int max_num_vertices_per_cell = 0;
80 for (size_t i = 0; i < basic_grid_data->count[YAC_LOC_CELL]; ++i)
81 if (basic_grid_data->num_vertices_per_cell[i] >
82 max_num_vertices_per_cell)
83 max_num_vertices_per_cell =
84 basic_grid_data->num_vertices_per_cell[i];
85 return max_num_vertices_per_cell;
86}
87
88static void get_cell_buffers(
89 struct yac_interp_grid * interp_grid, size_t max_num_src_per_tgt,
90 struct yac_grid_cell * tgt_grid_cell, struct yac_grid_cell ** src_grid_cells) {
91
92 struct yac_const_basic_grid_data * src_basic_grid_data =
94 struct yac_const_basic_grid_data * tgt_basic_grid_data =
96
97 *src_grid_cells = xmalloc(max_num_src_per_tgt * sizeof(**src_grid_cells));
98 enum yac_edge_type * edge_type_buffer;
99 double (*coordinates_xyz_buffer)[3];
100
101 // prepare source grid cell buffer
102 {
103 int max_num_vertices_per_cell =
104 MAX(get_max_num_vertices_per_cell(src_basic_grid_data),
105 get_max_num_vertices_per_cell(tgt_basic_grid_data));
106
107 edge_type_buffer =
108 xmalloc((max_num_src_per_tgt + 1) * (size_t)max_num_vertices_per_cell *
109 sizeof(*edge_type_buffer));
110 coordinates_xyz_buffer =
111 xmalloc((max_num_src_per_tgt + 1) * (size_t)max_num_vertices_per_cell *
112 sizeof(*coordinates_xyz_buffer));
113
114 tgt_grid_cell->coordinates_xyz = coordinates_xyz_buffer;
115 tgt_grid_cell->edge_type = edge_type_buffer;
116 tgt_grid_cell->array_size = max_num_vertices_per_cell;
117 for (size_t i = 0; i < max_num_src_per_tgt; ++i) {
118 (*src_grid_cells)[i].coordinates_xyz =
119 coordinates_xyz_buffer + (i + 1) * max_num_vertices_per_cell;
120 (*src_grid_cells)[i].edge_type =
121 edge_type_buffer + (i + 1) * max_num_vertices_per_cell;
122 (*src_grid_cells)[i].array_size = max_num_vertices_per_cell;
123 }
124 }
125}
126
128 struct yac_interp_grid * interp_grid,
129 struct yac_grid_cell * tgt_grid_cell, struct yac_grid_cell * src_grid_cell) {
130
131 struct yac_const_basic_grid_data * src_basic_grid_data =
133 struct yac_const_basic_grid_data * tgt_basic_grid_data =
135
136 int max_num_vertices_per_cell =
137 MAX(get_max_num_vertices_per_cell(src_basic_grid_data),
138 get_max_num_vertices_per_cell(tgt_basic_grid_data));
139
140 enum yac_edge_type * edge_type_buffer =
141 xmalloc(2 * (size_t)max_num_vertices_per_cell *
142 sizeof(*edge_type_buffer));
143 yac_coordinate_pointer coordinates_xyz_buffer =
144 xmalloc(2 * (size_t)max_num_vertices_per_cell *
145 sizeof(*coordinates_xyz_buffer));
146
147 tgt_grid_cell->coordinates_xyz = coordinates_xyz_buffer;
148 tgt_grid_cell->edge_type = edge_type_buffer;
149 tgt_grid_cell->array_size = max_num_vertices_per_cell;
150
151 src_grid_cell->coordinates_xyz =
152 coordinates_xyz_buffer + max_num_vertices_per_cell;
153 src_grid_cell->edge_type =
154 edge_type_buffer + max_num_vertices_per_cell;
155 src_grid_cell->array_size = max_num_vertices_per_cell;
156}
157
159 struct yac_const_basic_grid_data * tgt_basic_grid_data, size_t tgt_cell,
160 struct yac_const_basic_grid_data * src_basic_grid_data, size_t src_count,
161 size_t * src_cells, struct yac_grid_cell tgt_grid_cell_buffer,
162 struct yac_grid_cell * src_grid_cell_buffer,
163 double * weights, size_t * num_weights, int partial_coverage,
165 int enforced_conserv) {
166
168 tgt_basic_grid_data, tgt_cell, &tgt_grid_cell_buffer);
169 for (size_t i = 0; i < src_count; ++i)
171 src_basic_grid_data, src_cells[i], src_grid_cell_buffer + i);
172
173 double * area = weights;
175 src_count, src_grid_cell_buffer, tgt_grid_cell_buffer, area);
176
177 size_t num_valid_weights = 0;
178 for (size_t i = 0; i < src_count; ++i) {
179
180 if (area[i] > 0.0) {
181 if (i != num_valid_weights) {
182 area[num_valid_weights] = area[i];
183 src_cells[num_valid_weights] = src_cells[i];
184 }
185 ++num_valid_weights;
186 }
187 }
188 *num_weights = num_valid_weights;
189 if (num_valid_weights == 0) return 0;
190
191 double tgt_cell_area = yac_huiliers_area(tgt_grid_cell_buffer);
192 double norm_factor;
193
195 (normalisation == YAC_INTERP_CONSERV_DESTAREA) ||
196 (normalisation == YAC_INTERP_CONSERV_FRACAREA),
197 "ERROR(compute_weights_order_first_conserv_no_partial): "
198 "invalid normalisation option in conservative remapping")
199 switch(normalisation) {
201 norm_factor = 1.0 / tgt_cell_area;
202 break;
203 default:
205 double fracarea = 0.0;
206 for (size_t i = 0; i < num_valid_weights; ++i) fracarea += area[i];
207 norm_factor = 1.0 / fracarea;
208 break;
209 }
210 };
211
212 if (partial_coverage) {
213 for (size_t i = 0; i < num_valid_weights; ++i) weights[i] *= norm_factor;
214 return 1;
215 } else {
216 double tgt_cell_area_diff = tgt_cell_area;
217 double area_tol = tgt_cell_area * AREA_TOL_FACTOR;
218 for (size_t i = 0; i < num_valid_weights; ++i) {
219 double curr_area = area[i];
220 tgt_cell_area_diff -= curr_area;
221 weights[i] = curr_area * norm_factor;
222 }
223 int successful = fabs(tgt_cell_area_diff) <= area_tol;
224 if (successful && enforced_conserv)
225 yac_correct_weights(num_valid_weights, weights);
226 return successful;
227 }
228}
229
230static size_t do_search_conserv_1st_order (struct interp_method * method,
231 struct yac_interp_grid * interp_grid,
232 size_t * tgt_points, size_t count,
233 struct yac_interp_weights * weights) {
234
235 struct interp_method_conserv * method_conserv =
236 (struct interp_method_conserv *)method;
237
239 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
240 "ERROR(do_search_conserv): invalid number of source fields")
241
244 "ERROR(do_search_conserv): unsupported source field location type")
245
248 "ERROR(do_search_conserv): unsupported target field location type")
249
250
251 size_t * src_cells = NULL;
252 size_t * num_src_per_tgt = xmalloc(count * sizeof(*num_src_per_tgt));
253
254 // search matching cells
256 interp_grid, tgt_points, count, &src_cells, num_src_per_tgt);
257
258 // we did a search on the interp_grid, therefore we have to re-get the basic
259 // grid data
260 struct yac_const_basic_grid_data * tgt_basic_grid_data =
262 struct yac_const_basic_grid_data * src_basic_grid_data =
264
265 size_t total_num_weights = 0;
266 size_t max_num_src_per_tgt = 0;
267 for (size_t i = 0; i < count; ++i) {
268 size_t curr_num_src_per_tgt = num_src_per_tgt[i];
269 if (curr_num_src_per_tgt > max_num_src_per_tgt)
270 max_num_src_per_tgt = curr_num_src_per_tgt;
271 total_num_weights += num_src_per_tgt[i];
272 }
273
274 // to ensure that the interpolation always procduces the same result, we
275 // sort the source cells for each target point by their global ids
276 {
277 yac_int * temp_src_global_ids =
278 xmalloc(max_num_src_per_tgt * sizeof(*temp_src_global_ids));
279
280 for (size_t i = 0, offset = 0; i < count; ++i) {
281
282 size_t curr_num_src_per_tgt = num_src_per_tgt[i];
283 size_t * curr_src_cells = src_cells + offset;
284 offset += curr_num_src_per_tgt;
285
286 for (size_t j = 0; j < curr_num_src_per_tgt; ++j)
287 temp_src_global_ids[j] =
288 src_basic_grid_data->ids[YAC_LOC_CELL][curr_src_cells[j]];
289
291 temp_src_global_ids, curr_num_src_per_tgt, curr_src_cells);
292 }
293
294 free(temp_src_global_ids);
295 }
296
297 double * w = xmalloc(total_num_weights * sizeof(*w));
298 size_t result_count = 0;
299 size_t * failed_tgt = xmalloc(count * sizeof(*failed_tgt));
300 total_num_weights = 0;
301
302 int partial_coverage = method_conserv->partial_coverage;
303 enum yac_interp_method_conserv_normalisation normalisation =
304 method_conserv->normalisation;
305 int enforced_conserv = method_conserv->enforced_conserv;
306
307 struct yac_grid_cell tgt_grid_cell;
308 struct yac_grid_cell * src_grid_cells;
310 interp_grid, max_num_src_per_tgt, &tgt_grid_cell, &src_grid_cells);
311
312 // compute overlaps
313 for (size_t i = 0, offset = 0, result_offset = 0; i < count; ++i) {
314
315 size_t curr_src_count = num_src_per_tgt[i];
316 size_t curr_tgt_point = tgt_points[i];
317 size_t num_weights;
318
319 // if weight computation was successful
321 tgt_basic_grid_data, curr_tgt_point,
322 src_basic_grid_data, curr_src_count, src_cells + offset,
323 tgt_grid_cell, src_grid_cells, w + result_offset, &num_weights,
324 partial_coverage, normalisation, enforced_conserv)) {
325
326 if (offset != result_offset) {
327
328 memmove(
329 src_cells + result_offset, src_cells + offset,
330 num_weights * sizeof(*src_cells));
331 }
332 tgt_points[result_count] = curr_tgt_point;
333 num_src_per_tgt[result_count] = num_weights;
334 result_count++;
335 result_offset += num_weights;
336 total_num_weights += num_weights;
337 } else {
338 failed_tgt[i - result_count] = curr_tgt_point;
339 }
340
341 offset += curr_src_count;
342 }
343
344 free(tgt_grid_cell.edge_type);
345 free(tgt_grid_cell.coordinates_xyz);
346 free(src_grid_cells);
347
348 if (result_count != count)
349 memcpy(tgt_points + result_count, failed_tgt,
350 (count - result_count) * sizeof(*tgt_points));
351 free(failed_tgt);
352
353 struct remote_points tgts = {
354 .data =
356 interp_grid, tgt_points, result_count),
357 .count = result_count};
358 struct remote_point * srcs =
360 interp_grid, 0, src_cells, total_num_weights);
361
362 // store weights
364 weights, &tgts, num_src_per_tgt, srcs, w);
365
366 free(tgts.data);
367 free(srcs);
368 free(src_cells);
369 free(num_src_per_tgt);
370 free(w);
371
372 return result_count;
373}
374
375static int
376compare_supermesh_cell_src_local_ids(const void * a, const void * b) {
377
378 struct supermesh_cell * a_ = (struct supermesh_cell *)a;
379 struct supermesh_cell * b_ = (struct supermesh_cell *)b;
380
381 int ret = (a_->src.local_id > b_->src.local_id) -
382 (a_->src.local_id < b_->src.local_id);
383 if (ret) return ret;
384 return (a_->tgt.global_id > b_->tgt.global_id) -
385 (a_->tgt.global_id < b_->tgt.global_id);
386}
387
388static int
389compare_supermesh_cell_tgt_local_ids(const void * a, const void * b) {
390
391 struct supermesh_cell * a_ = (struct supermesh_cell *)a;
392 struct supermesh_cell * b_ = (struct supermesh_cell *)b;
393
394 int ret = (a_->tgt.local_id > b_->tgt.local_id) -
395 (a_->tgt.local_id < b_->tgt.local_id);
396 if (ret) return ret;
397 return (a_->src.global_id > b_->src.global_id) -
398 (a_->src.global_id < b_->src.global_id);
399}
400
402 double * src_cell_centroid, struct weight_vector_3d * G_i,
403 struct weight_vector_data_3d * buffer) {
404
405 // This routine computes: (I_3 - C_i * C_i^-1) * G_i
406 //
407 // O(g_i) = g_i - C_i * (C_i^-1 * g_i)
408 // where: C_i is the centeroid of a source cell
409 // g_i is the gradient of the source field in C_i
410 // O(g_i) is the projection of g_i into the plane perpendicular to C_i
411 // g_i = G_i * f
412 // where: G_i is the weight matrix to compute g_i
413 // f is the source field vector
414 // => O(g_i) = (I_3 - C_i * C_i^-1) * G_i * f
415 // where: I_3 is the identity matrix of size 3 x 3
416 //
417 // M = I_3 - C_i * C_i^-1
418
419 double M[3][3];
420 for (size_t k = 0; k < 3; ++k)
421 for (size_t l = 0; l < 3; ++l)
422 M[k][l] = - src_cell_centroid[k] * src_cell_centroid[l];
423 for (size_t k = 0; k < 3; ++k)
424 M[k][k] += 1.0;
425
426 struct weight_vector_data_3d * G_i_data = G_i->data;
427
428 size_t N = G_i->n;
429 for (size_t i = 0; i < N; ++i) {
430 buffer[i].local_id = G_i_data[i].local_id;
431 buffer[i].global_id = G_i_data[i].global_id;
432 for (size_t j = 0; j < 3; ++j) buffer[i].weight[j] = 0.0;
433 }
434
435 for (size_t n = 0; n < N; ++n)
436 for (size_t i = 0; i < 3; ++i)
437 for (size_t j = 0; j < 3; ++j)
438 buffer[n].weight[i] += G_i_data[n].weight[j] * M[i][j];
439
440 memcpy(G_i->data, buffer, N * sizeof(*buffer));
441}
442
444 void const * a, void const * b) {
445
446 struct weight_vector_data const * weight_a =
447 (struct weight_vector_data const *)a;
448 struct weight_vector_data const * weight_b =
449 (struct weight_vector_data const *)b;
450
451 int ret = weight_a->global_id - weight_b->global_id;
452 if (ret) return ret;
453 double abs_weight_a = fabs(weight_a->weight);
454 double abs_weight_b = fabs(weight_b->weight);
455 ret = (abs_weight_a > abs_weight_b) - (abs_weight_a < abs_weight_b);
456 if (ret) return ret;
457 return (weight_a->weight > weight_b->weight) -
458 (weight_a->weight < weight_b->weight);
459}
460
462 void const * a, void const * b) {
463
464 struct weight_vector_data const * weight_a =
465 (struct weight_vector_data const *)a;
466 struct weight_vector_data const * weight_b =
467 (struct weight_vector_data const *)b;
468
469 return weight_a->global_id - weight_b->global_id;
470}
471
473 struct weight_vector_data * weights, size_t * n) {
474
475 size_t n_ = *n;
476
477 if (n_ <= 1) return;
478
479 // sort weights by global_id then by weight
480 qsort(weights, n_, sizeof(*weights), compare_weight_vector_data_weight);
481
482 size_t new_n = 1;
483 struct weight_vector_data * prev_weight_data = weights;
484 struct weight_vector_data * curr_weight_data = weights + 1;
485 for (size_t i = 1; i < n_; ++i, ++curr_weight_data) {
486
487 // if both weights refer to the same source point (by global_id)
488 if (!compare_weight_vector_data(prev_weight_data, curr_weight_data)) {
489 prev_weight_data->weight += curr_weight_data->weight;
490 } else {
491 ++new_n;
492 ++prev_weight_data;
493 *prev_weight_data = *curr_weight_data;
494 }
495 }
496
497 n_ = new_n;
498 new_n = 0;
499
500 // check for zero-weights
501 for (size_t i = 0; i < n_; ++i) {
502
503 if (weights[i].weight == 0.0) continue;
504 if (i != new_n) weights[new_n] = weights[i];
505 ++new_n;
506 }
507
508 *n = new_n;
509}
510
512 struct supermesh_cell * super_cell, struct weight_vector_data * weights) {
513
514 weights[0].weight = super_cell->norm_area;
515 weights[0].global_id = super_cell->src.global_id;
516 weights[0].local_id = super_cell->src.local_id;
517
518 struct weight_vector_3d * src_cell_gradient = super_cell->src_cell_gradient;
519
520 size_t N = src_cell_gradient->n;
521 // in case we have no gradient for the current supermesh cell,
522 // we assume a constant field across the whole associated source cell
523 if (N > 0) {
524 struct weight_vector_data_3d * gradient_weights = src_cell_gradient->data - 1;
525 double * overlap_barycenter = super_cell->barycenter;
526
527 for (size_t n = 1; n <= N; ++n) {
528 weights[n].weight =
529 (gradient_weights[n].weight[0] * overlap_barycenter[0] +
530 gradient_weights[n].weight[1] * overlap_barycenter[1] +
531 gradient_weights[n].weight[2] * overlap_barycenter[2]) *
532 super_cell->norm_area;
533 weights[n].global_id = gradient_weights[n].global_id;
534 weights[n].local_id = gradient_weights[n].local_id;
535 }
536 }
537
538 return 1 + N;
539}
540
542 struct yac_const_basic_grid_data * grid_data, size_t cell_idx,
543 double barycenter[3]) {
544
545 size_t num_vertices = grid_data->num_vertices_per_cell[cell_idx];
546 size_t const * vertices =
547 grid_data->cell_to_vertex + grid_data->cell_to_vertex_offsets[cell_idx];
548
549 barycenter[0] = 0.0;
550 barycenter[1] = 0.0;
551 barycenter[2] = 0.0;
552
553 for (size_t i = 0; i < num_vertices; ++i) {
554 double const * curr_vertex_coordinate =
555 grid_data->vertex_coordinates[vertices[i]];
556 barycenter[0] += curr_vertex_coordinate[0];
557 barycenter[1] += curr_vertex_coordinate[1];
558 barycenter[2] += curr_vertex_coordinate[2];
559 }
560 normalise_vector(barycenter);
561}
562
564 struct yac_interp_grid * interp_grid, size_t * tgt_points, size_t count,
565 struct supermesh_cell ** super_cells_, size_t * num_super_cells,
566 int * interp_fail_flag, size_t ** src_cells, size_t * num_src_cells,
568 int partial_coverage) {
569
571 (normalisation == YAC_INTERP_CONSERV_DESTAREA) ||
572 (normalisation == YAC_INTERP_CONSERV_FRACAREA),
573 "ERROR(compute_super_cells): "
574 "invalid normalisation option in conservative remapping")
575
576 size_t * num_src_per_tgt = xmalloc(count * sizeof(*num_src_per_tgt));
577
578 // search for all source cell overlapping with the the target cells
580 interp_grid, tgt_points, count, src_cells, num_src_per_tgt);
581
582 // determine the number of unique matching source cells
583 size_t total_num_overlaps = 0;
584 for (size_t i = 0; i < count; ++i) total_num_overlaps += num_src_per_tgt[i];
585 yac_quicksort_index_size_t_int(*src_cells, total_num_overlaps, NULL);
586 *num_src_cells = total_num_overlaps;
587 yac_remove_duplicates_size_t(*src_cells, num_src_cells);
588
589 size_t * num_tgt_per_src =
590 xrealloc(num_src_per_tgt, *num_src_cells * sizeof(*num_tgt_per_src));
591
592 // for some required source cells we may have not all required supermesh cells
593 // in that case we have to find the respective target cells and compute the
594 // missing supermesh cells from them
595 size_t * tgt_cells = NULL;
597 interp_grid, *src_cells, *num_src_cells, &tgt_cells, num_tgt_per_src);
598
599 total_num_overlaps = 0;
600 for (size_t i = 0; i < *num_src_cells; ++i)
601 total_num_overlaps += num_tgt_per_src[i];
602
603 struct supermesh_cell * super_cells =
604 xmalloc(total_num_overlaps * sizeof(*super_cells));
605
606 struct yac_const_basic_grid_data * src_basic_grid_data =
608 struct yac_const_basic_grid_data * tgt_basic_grid_data =
610
611 for (size_t i = 0, j = 0; i < *num_src_cells; ++i) {
612
613 size_t curr_num_overlaps = num_tgt_per_src[i];
614 size_t curr_src_cell = (*src_cells)[i];
615 yac_int curr_src_global_id =
616 src_basic_grid_data->ids[YAC_LOC_CELL][curr_src_cell];
617
618 for (size_t k = 0; k < curr_num_overlaps; ++k, ++j) {
619
620 size_t curr_tgt_cell = tgt_cells[j];
621 struct supermesh_cell * curr_super_cell = super_cells + j;
622 curr_super_cell->src.local_id = curr_src_cell;
623 curr_super_cell->src.global_id = curr_src_global_id;
624 curr_super_cell->tgt.local_id = curr_tgt_cell;
625 curr_super_cell->tgt.global_id =
626 tgt_basic_grid_data->ids[YAC_LOC_CELL][curr_tgt_cell];
627 }
628 }
629 free(tgt_cells);
630 free(num_tgt_per_src);
631
632 // sort supermesh_cell first by local ids of the target cells and
633 // second by global cell id
634 qsort(super_cells, total_num_overlaps, sizeof(*super_cells),
636
637 struct yac_grid_cell tgt_grid_cell;
638 struct yac_grid_cell src_grid_cell;
639 get_cell_buffers_(interp_grid, &tgt_grid_cell, &src_grid_cell);
640
641 src_basic_grid_data = yac_interp_grid_get_basic_grid_data_src(interp_grid);
642 tgt_basic_grid_data = yac_interp_grid_get_basic_grid_data_tgt(interp_grid);
643
644 // For all supermesh cells compute the area and normalised area.
645 // Additionally, remove all empty supermesh cells.
646 size_t new_num_super_cells = 0;
647 size_t tgt_idx = 0;
648 for (size_t i = 0, j = 0; i < total_num_overlaps;) {
649
650 // get information about the current target cell
651 size_t curr_tgt_cell = super_cells[i].tgt.local_id;
653 tgt_basic_grid_data, curr_tgt_cell, &tgt_grid_cell);
654 double curr_tgt_cell_coverage = 0.0;
655
656 // for all supermesh cells overlapping with the current target cell
657 for (;(i < total_num_overlaps) &&
658 (super_cells[i].tgt.local_id == curr_tgt_cell); ++i) {
659
660 // get the current source cell
662 src_basic_grid_data, super_cells[i].src.local_id, &src_grid_cell);
663
664 // compute area of the current supermesh cell
665 double super_cell_area;
666 double barycenter[3];
668 1, &src_grid_cell, tgt_grid_cell, &super_cell_area, &barycenter);
669
670 // if there is an overlap between the current source and target cell
671 if (super_cell_area > 0.0) {
672
673 super_cells[new_num_super_cells].src = super_cells[i].src;
674 super_cells[new_num_super_cells].tgt = super_cells[i].tgt;
675 super_cells[new_num_super_cells].area = super_cell_area;
676 memcpy(super_cells[new_num_super_cells].barycenter, barycenter,
677 3 * sizeof(double));
678 super_cells[new_num_super_cells].src_cell_gradient = NULL;
679 ++new_num_super_cells;
680
681 curr_tgt_cell_coverage += super_cell_area;
682 }
683 }
684
685 double curr_tgt_cell_area = yac_huiliers_area(tgt_grid_cell);
686
687 // if there was an overlap
688 if (new_num_super_cells != j) {
689
690 double norm_factor;
692 (normalisation == YAC_INTERP_CONSERV_DESTAREA) ||
693 (normalisation == YAC_INTERP_CONSERV_FRACAREA),
694 "ERROR(compute_super_cells): invalid normalisation")
695 switch (normalisation) {
696 default:
698 norm_factor = 1.0 / curr_tgt_cell_area;
699 break;
701 norm_factor = 1.0 / curr_tgt_cell_coverage;
702 break;
703 }
704 // compute normalised area
705 for (; j < new_num_super_cells; ++j)
706 super_cells[j].norm_area = super_cells[j].area * norm_factor;
707 }
708
709
710 // for target cell that do not overlap with any source cell
711 while ((tgt_idx < count) && (tgt_points[tgt_idx] < curr_tgt_cell))
712 interp_fail_flag[tgt_idx++] = 1;
713
714 if ((tgt_idx < count) && (tgt_points[tgt_idx] == curr_tgt_cell)) {
715
716 double area_tol = curr_tgt_cell_area * AREA_TOL_FACTOR;
717
718 if (partial_coverage) {
719 interp_fail_flag[tgt_idx] = curr_tgt_cell_coverage < area_tol;
720 } else {
721 interp_fail_flag[tgt_idx] =
722 fabs(curr_tgt_cell_area - curr_tgt_cell_coverage) > area_tol;
723 }
724 ++tgt_idx;
725 }
726 }
727 // for all remaining target cell that do not overlap with any source cell
728 for (; tgt_idx < count; ++tgt_idx) interp_fail_flag[tgt_idx] = 1;
729 *num_super_cells = new_num_super_cells;
730 *super_cells_ =
731 xrealloc(super_cells, new_num_super_cells * sizeof(*super_cells));
732 free(tgt_grid_cell.coordinates_xyz);
733 free(tgt_grid_cell.edge_type);
734}
735
737 struct yac_interp_grid * interp_grid,
738 size_t * src_cells, int * skip_src_cell, size_t num_src_cells,
739 struct supermesh_cell * super_cells, size_t num_super_cells) {
740
741 yac_coordinate_pointer src_cell_centroids =
742 xmalloc(num_src_cells * sizeof(*src_cell_centroids));
743
744 // sort supermesh cells first by source local id and second by
745 // target global id
746 qsort(super_cells, num_super_cells, sizeof(*super_cells),
748
749 struct yac_grid_cell src_grid_cell, dummy;
750 get_cell_buffers_(interp_grid, &src_grid_cell, &dummy);
751
752 struct yac_const_basic_grid_data * src_basic_grid_data =
754
755 // compute centroids of source cells
756 // C_i = N(S_(U_k) (A_k*C_k))
757 // where: N(C) = C*(C*C)^-0.5 // normalisation
758 // U_k all supermesh cells overlaping with the respective source cell
759 // A_k area of supermesh cell
760 // C_k barycenter of supermesh cell (normalisation of the sum of all
761 // vertices of the cell)
762 for (size_t i = 0, offset = 0; i < num_src_cells; ++i) {
763
764 if (skip_src_cell[i]) continue;
765
766 size_t curr_src_cell = src_cells[i];
768 src_basic_grid_data, curr_src_cell, &src_grid_cell);
769 double src_cell_area = yac_huiliers_area(src_grid_cell);
770
771 struct supermesh_cell * curr_super_cells = super_cells + offset;
772 size_t curr_num_super_cells = offset;
773 while ((offset < num_super_cells) &&
774 (super_cells[offset].src.local_id == curr_src_cell)) ++offset;
775 curr_num_super_cells = offset - curr_num_super_cells;
776
777 double src_cell_area_diff = src_cell_area;
778 double src_cell_centroid[3] = {0.0, 0.0, 0.0};
779
780 for (size_t j = 0; j < curr_num_super_cells; ++j) {
781
782 double super_cell_area = curr_super_cells[j].area;
783 double * super_cell_barycenter = curr_super_cells[j].barycenter;
784 src_cell_centroid[0] += super_cell_area * super_cell_barycenter[0];
785 src_cell_centroid[1] += super_cell_area * super_cell_barycenter[1];
786 src_cell_centroid[2] += super_cell_area * super_cell_barycenter[2];
787 src_cell_area_diff -= super_cell_area;
788 }
789
790 normalise_vector(src_cell_centroid);
791 src_cell_centroids[i][0] = src_cell_centroid[0];
792 src_cell_centroids[i][1] = src_cell_centroid[1];
793 src_cell_centroids[i][2] = src_cell_centroid[2];
794
795 }
796 free(src_grid_cell.coordinates_xyz);
797 free(src_grid_cell.edge_type);
798
799 return src_cell_centroids;
800}
801
803 struct yac_interp_grid * interp_grid, size_t * src_cells,
804 yac_coordinate_pointer src_cell_centroids, int * skip_src_cell,
805 size_t num_src_cells, size_t * src_cell_neighbours) {
806
807 struct weight_vector_3d * src_cell_gradients =
808 xmalloc(num_src_cells * sizeof(*src_cell_gradients));
809
810 struct yac_const_basic_grid_data * src_basic_grid_data =
812
813 size_t total_num_gradient_weights = 0;
814 size_t max_num_neigh_per_src = 0;
815 for (size_t i = 0; i < num_src_cells; ++i) {
816 src_cell_gradients[i].data = NULL;
817 src_cell_gradients[i].n = 0;
818 size_t curr_num_neigh =
819 src_basic_grid_data->num_vertices_per_cell[src_cells[i]];
820 total_num_gradient_weights += curr_num_neigh + 1;
821 if (max_num_neigh_per_src < curr_num_neigh)
822 max_num_neigh_per_src = curr_num_neigh;
823 }
824 struct weight_vector_data_3d * weight_vector_data_buffer =
825 (total_num_gradient_weights > 0)?
826 (xmalloc(total_num_gradient_weights *
827 sizeof(*weight_vector_data_buffer))):NULL;
828 for (size_t i = 0; i < total_num_gradient_weights; ++i)
829 for (size_t j = 0; j < 3; ++j)
830 weight_vector_data_buffer[i].weight[j] = 0.0;
831
832 struct weight_vector_data_3d * orth_buffer =
833 xmalloc((max_num_neigh_per_src + 1) * sizeof(*orth_buffer));
834
835 // compute gradient in the centroid for each source cell
836 // g_i = O_i(g'_i)
837 // where: O_i(g) = g - C_i*(C_i*g) // makes g orthogonal to C_i
838 // g'_i = (A_(C_k))^-1 * S_ijk(((f_j + f_k) * 0.5 - f_i) *
839 // (C_j x C_k) * |C_j x C_k|^-1 *
840 // asin(|C_j x C_k|))
841 // where: g'_i is the estimated centroid gradient
842 // A_(C_k) is the area of the cell generated by connecting the
843 // barycenters of all neighbour cells
844 // S_ijk is the sum of all edges of the previously described cell in
845 // counterclockwise order
846 // f_i, f_j, and f_k are the mean values of the field over the area of
847 // the respective source cell area
848 // asin(|C_j x C_k|) / |C_j x C_k| ~ 1.0
849 // => g'_i = S_ijk(((f_j + f_k) * 0.5 - f_i) * (C_j x C_k)) / A_(C_k)
850 //
851 // g_i = G_i * f
852 // G_i = S_ijk(((I_j + I_k) * 0.5 - I_i) * (C_j x C_k)) / A_(C_k)
853 // where: G_i is the gradient weight matrix
854 // I_x is a vector of the same size as f, which is all zero except at
855 // position x, where it is one
856 // f is the source field vector
857 for (size_t i = 0, offset = 0, weight_vector_data_buffer_offset = 0;
858 i < num_src_cells; ++i) {
859
860 size_t curr_num_neigh =
861 src_basic_grid_data->num_vertices_per_cell[src_cells[i]];
862 size_t * curr_neighs = src_cell_neighbours + offset;
863 offset += curr_num_neigh;
864 struct weight_vector_3d * G_i = src_cell_gradients + i;
865 G_i->data = weight_vector_data_buffer + weight_vector_data_buffer_offset;
866
867 if (skip_src_cell[i]) continue;
868
869 weight_vector_data_buffer_offset += curr_num_neigh + 1;
870 G_i->n = curr_num_neigh + 1;
871 G_i->data[0].local_id = src_cells[i];
872 G_i->data[0].global_id =
873 src_basic_grid_data->ids[YAC_LOC_CELL][src_cells[i]];
874 for (size_t j = 0; j < curr_num_neigh; ++j) {
875 size_t curr_neigh = curr_neighs[j];
876 // if the current edge has a neighbour
877 if (curr_neigh != SIZE_MAX) {
878 G_i->data[j+1].local_id = curr_neigh;
879 G_i->data[j+1].global_id =
880 src_basic_grid_data->ids[YAC_LOC_CELL][curr_neigh];
881 } else {
882 // if the current edge has no neighbour, use current cell instead
883 G_i->data[j+1].local_id = G_i->data[0].local_id;
884 G_i->data[j+1].global_id = G_i->data[0].global_id;
885 }
886 }
887
888 // area of the polygon that is formed by connecting the barycenters of
889 // the neighbouring cells
890 double A_C_k = 0.0;
891
892 struct yac_grid_cell centroid_triangle = {
893 .coordinates_xyz = (double[3][3]){{0}},
894 .edge_type =
896 .num_corners = 3, .array_size = 0};
897 centroid_triangle.coordinates_xyz[0][0] = src_cell_centroids[i][0];
898 centroid_triangle.coordinates_xyz[0][1] = src_cell_centroids[i][1];
899 centroid_triangle.coordinates_xyz[0][2] = src_cell_centroids[i][2];
900
901 double edge_direction = 0.0;
902
903 // We split the cell that is comprised of the barycenters of the edge
904 // neigbours into triangles, which has the centroid of the current cell
905 // as one corner. The sum of the areas of these triangles is A_C_K.
906 for (size_t j = 0; j < curr_num_neigh; ++j) {
907
908 size_t neigh_idx[2] = {j + 1, (j+1)%curr_num_neigh+1};
909 size_t neigh_local_ids[2] =
910 {G_i->data[neigh_idx[0]].local_id,
911 G_i->data[neigh_idx[1]].local_id};
912
913 if (neigh_local_ids[0] == neigh_local_ids[1]) continue;
914
915 double neigh_cell_barycenters[2][3];
917 src_basic_grid_data, neigh_local_ids[0], neigh_cell_barycenters[0]);
919 src_basic_grid_data, neigh_local_ids[1], neigh_cell_barycenters[1]);
920
921 centroid_triangle.coordinates_xyz[1][0] = neigh_cell_barycenters[0][0];
922 centroid_triangle.coordinates_xyz[1][1] = neigh_cell_barycenters[0][1];
923 centroid_triangle.coordinates_xyz[1][2] = neigh_cell_barycenters[0][2];
924 centroid_triangle.coordinates_xyz[2][0] = neigh_cell_barycenters[1][0];
925 centroid_triangle.coordinates_xyz[2][1] = neigh_cell_barycenters[1][1];
926 centroid_triangle.coordinates_xyz[2][2] = neigh_cell_barycenters[1][2];
927
928 A_C_k += yac_huiliers_area(centroid_triangle);
929
930 // C_j x C_k
931 double C_j_x_C_k[3];
933 neigh_cell_barycenters[0], neigh_cell_barycenters[1], C_j_x_C_k);
934
935 double curr_edge_direction =
936 C_j_x_C_k[0] * centroid_triangle.coordinates_xyz[0][0] +
937 C_j_x_C_k[1] * centroid_triangle.coordinates_xyz[0][1] +
938 C_j_x_C_k[2] * centroid_triangle.coordinates_xyz[0][2];
939
940 if (fabs(curr_edge_direction) > fabs(edge_direction))
941 edge_direction = curr_edge_direction;
942
943 // -I_i * (C_j x C_k)
944 G_i->data[0].weight[0] -= C_j_x_C_k[0];
945 G_i->data[0].weight[1] -= C_j_x_C_k[1];
946 G_i->data[0].weight[2] -= C_j_x_C_k[2];
947
948 for (size_t l = 0; l < 3; ++l) C_j_x_C_k[l] *= 0.5;
949
950 // 0.5 * I_j * (C_j x C_k)
951 G_i->data[neigh_idx[0]].weight[0] += C_j_x_C_k[0];
952 G_i->data[neigh_idx[0]].weight[1] += C_j_x_C_k[1];
953 G_i->data[neigh_idx[0]].weight[2] += C_j_x_C_k[2];
954
955 // 0.5 * I_k * (C_j x C_k)
956 G_i->data[neigh_idx[1]].weight[0] += C_j_x_C_k[0];
957 G_i->data[neigh_idx[1]].weight[1] += C_j_x_C_k[1];
958 G_i->data[neigh_idx[1]].weight[2] += C_j_x_C_k[2];
959 } // curr_num_neigh
960
961 // if the neighbours were ordered in the wrong direction
962 if (edge_direction > 0.0)
963 for (size_t j = 0; j <= curr_num_neigh; ++j)
964 for (size_t k = 0; k < 3; ++k)
965 G_i->data[j].weight[k] *= -1.0;
966
967 double inv_A_C_K = (A_C_k > YAC_AREA_TOL)?(1.0/A_C_k):0.0;
968 // ((I_j + I_k) * 0.5 - I_i) * (C_j x C_k) / A_(C_k)
969 for (size_t k = 0; k <= curr_num_neigh; ++k)
970 for (size_t l = 0; l < 3; ++l)
971 G_i->data[k].weight[l] *= inv_A_C_K;
972
974 src_cell_centroids[i], G_i, orth_buffer);
975 } // num_src_cells
976 free(orth_buffer);
977
978 return src_cell_gradients;
979}
980
982 size_t * tgt_cells, int * interp_fail_flag, size_t num_tgt_cells,
983 struct supermesh_cell * super_cells, size_t num_super_cells,
984 size_t ** src_per_tgt, double ** weights, size_t * num_src_per_tgt) {
985
986 size_t num_interpolated_tgt = 0;
987
988 // sort supermesh cells first by target local id and second by
989 // source global id
990 qsort(super_cells, num_super_cells, sizeof(*super_cells),
992
993 size_t max_num_weights_per_tgt = 0;
994 size_t max_num_total_weights = 0;
995
996 // count maximum number of weights
997 for (size_t i = 0, j = 0; i < num_tgt_cells; ++i) {
998
999 if (interp_fail_flag[i]) continue;
1000
1001 size_t curr_tgt_cell = tgt_cells[i];
1002 size_t curr_num_weights = 0;
1003
1004 // skip supermesh cells not overlapping with current target cell
1005 while ((j < num_super_cells) &&
1006 (super_cells[j].tgt.local_id < curr_tgt_cell)) ++j;
1007
1008 // for all supermesh cells overlapping with the current target cell
1009 while ((j < num_super_cells) &&
1010 (super_cells[j].tgt.local_id == curr_tgt_cell))
1011 curr_num_weights += 1 + super_cells[j++].src_cell_gradient->n;
1012
1013 max_num_total_weights += curr_num_weights;
1014 if (max_num_weights_per_tgt < curr_num_weights)
1015 max_num_weights_per_tgt = curr_num_weights;
1016 num_interpolated_tgt++;
1017 }
1018
1019 struct weight_vector_data * weight_buffer =
1020 xmalloc(max_num_weights_per_tgt * sizeof(*weight_buffer));
1021 *weights = xmalloc(max_num_total_weights * sizeof(**weights));
1022 *src_per_tgt = xmalloc(max_num_total_weights * sizeof(**src_per_tgt));
1023
1024 // sort all target points that can be interpolated to the beginning of the
1025 // tgt_cells array
1026 yac_quicksort_index_int_size_t(interp_fail_flag, num_tgt_cells, tgt_cells);
1027 yac_quicksort_index_size_t_int(tgt_cells, num_interpolated_tgt, NULL);
1028
1029 // compute 2nd order weights
1030 // f_j = SUM(f_k') / A_j
1031 // f_k' = A_k*(f_i + g_i * C_k)
1032 //
1033 // f_k' = A_k * (I_i + C_k * G_i) * f
1034 // where: I_i is a vector of the same length as f, that is all 0.0 except at
1035 // position i, where it is 1.0
1036 // f is the source field vector
1037 // => f_k' / A_j = M_k * f
1038 // where: M_k = A_k / A_j * (I_i + C_k * G_i)
1039 // => f_j = SUM(M_k) * f
1040 // f_j = M_j * f
1041 // where: M_j = SUM(M_k)
1042 size_t w_idx = 0;
1043 for (size_t i = 0, j = 0; i < num_interpolated_tgt; ++i) {
1044
1045 size_t curr_tgt_cell = tgt_cells[i];
1046
1047 // skip supermesh cells not overlapping with current target cell
1048 while ((j < num_super_cells) &&
1049 (super_cells[j].tgt.local_id != curr_tgt_cell)) ++j;
1050
1051 size_t weight_buffer_offset = 0;
1052
1053 // for all supermesh cells overlapping with the current target cell
1054 while ((j < num_super_cells) &&
1055 (super_cells[j].tgt.local_id == curr_tgt_cell)) {
1056
1057 size_t num_weights =
1059 super_cells + j, weight_buffer + weight_buffer_offset);
1060 weight_buffer_offset += num_weights;
1061 ++j;
1062 }
1063
1064 // compact weights and sort them by global_id
1065 compact_weight_vector_data(weight_buffer, &weight_buffer_offset);
1066
1067 num_src_per_tgt[i] = weight_buffer_offset;
1068
1069 for (size_t k = 0; k < weight_buffer_offset; ++k, ++w_idx) {
1070 (*weights)[w_idx] = weight_buffer[k].weight;
1071 (*src_per_tgt)[w_idx] = weight_buffer[k].local_id;
1072 }
1073 }
1074 free(weight_buffer);
1075
1076 return num_interpolated_tgt;
1077}
1078
1079static size_t do_search_conserv_2nd_order (struct interp_method * method,
1080 struct yac_interp_grid * interp_grid,
1081 size_t * tgt_points, size_t count,
1082 struct yac_interp_weights * weights) {
1083
1084 struct interp_method_conserv * method_conserv =
1085 (struct interp_method_conserv *)method;
1086
1087 YAC_ASSERT(
1088 yac_interp_grid_get_num_src_fields(interp_grid) == 1,
1089 "ERROR(do_search_conserv): invalid number of source fields")
1090
1091 YAC_ASSERT(
1093 "ERROR(do_search_conserv): unsupported source field location type")
1094
1095 YAC_ASSERT(
1097 "ERROR(do_search_conserv): unsupported target field location type")
1098
1099 // sort target points
1100 yac_quicksort_index_size_t_int(tgt_points, count, NULL);
1101
1102 int * interp_fail_flag = xmalloc(count * sizeof(*interp_fail_flag));
1103 struct supermesh_cell * super_cells = NULL;
1104 size_t num_super_cells = 0;
1105 size_t * src_cells = NULL;
1106 size_t num_src_cells = 0;
1107
1108 // compute the overlaps between the target cells and the source grid
1110 interp_grid, tgt_points, count, &super_cells, &num_super_cells,
1111 interp_fail_flag, &src_cells, &num_src_cells,
1112 method_conserv->normalisation, method_conserv->partial_coverage);
1113
1114 size_t total_num_src_cell_neighbours = 0;
1115 struct yac_const_basic_grid_data * src_basic_grid_data =
1117 for (size_t i = 0; i < num_src_cells; ++i)
1118 total_num_src_cell_neighbours +=
1119 src_basic_grid_data->num_vertices_per_cell[src_cells[i]];
1120 size_t * src_cell_neighbours =
1121 xmalloc(total_num_src_cell_neighbours * sizeof(*src_cell_neighbours));
1122
1123 // get the neighbours for all source cells overlapping with a target cell
1125 interp_grid, src_cells, num_src_cells, src_cell_neighbours);
1126
1127 const_int_pointer src_cell_mask =
1128 yac_interp_grid_get_src_field_mask(interp_grid, 0);
1129 int * skip_src_cell = xmalloc(num_src_cells * sizeof(*skip_src_cell));
1130
1131 // check whether any required source cell or any of its neighbours is masked
1132 // out in the field mask (deactivate respective cells)
1133 if (src_cell_mask != NULL) {
1134 for (size_t i = 0, offset = 0; i < num_src_cells; ++i) {
1135 size_t curr_src_cell = src_cells[i];
1136 if ((skip_src_cell[i] = !src_cell_mask[curr_src_cell])) continue;
1137 size_t * curr_neighbours = src_cell_neighbours + offset;
1138 size_t curr_num_neigh =
1139 src_basic_grid_data->num_vertices_per_cell[curr_src_cell];
1140 offset += curr_num_neigh;
1141 for (size_t j = 0; j < curr_num_neigh; ++j)
1142 if ((curr_neighbours[j] != SIZE_MAX) &&
1143 (!src_cell_mask[curr_neighbours[j]]))
1144 curr_neighbours[j] = SIZE_MAX;
1145 }
1146 } else {
1147 memset(skip_src_cell, 0, num_src_cells * sizeof(*skip_src_cell));
1148 }
1149
1150 // compute the centroids of all required source cells
1151 yac_coordinate_pointer src_cell_centroids =
1152 compute_src_cell_centroids(interp_grid, src_cells, skip_src_cell,
1153 num_src_cells, super_cells, num_super_cells);
1154
1155 // compute the gradients weights for the source cells
1156 struct weight_vector_3d * src_cell_gradients =
1158 interp_grid, src_cells, src_cell_centroids,
1159 skip_src_cell, num_src_cells, src_cell_neighbours);
1160 free(src_cell_neighbours);
1161 free(src_cell_centroids);
1162
1163 // sort supermesh cells by local src id
1164 qsort(super_cells, num_super_cells, sizeof(*super_cells),
1166 for (size_t i = 0, j = 0; i < num_src_cells; ++i) {
1167 size_t curr_src_cell = src_cells[i];
1168 struct weight_vector_3d * curr_src_cell_gradient = src_cell_gradients + i;
1169 // there should not be any supercell whose source cell is not in the list of
1170 // required source cells
1171 YAC_ASSERT(
1172 (j >= num_super_cells) ||
1173 (super_cells[j].src.local_id >= curr_src_cell),
1174 "ERROR(do_search_conserv_2nd_order): internal error");
1175 while ((j < num_super_cells) &&
1176 (super_cells[j].src.local_id == curr_src_cell)) {
1177 super_cells[j++].src_cell_gradient = curr_src_cell_gradient;
1178 }
1179 }
1180 free(skip_src_cell);
1181 free(src_cells);
1182
1183 size_t * src_per_tgt = NULL;
1184 double * w = NULL;
1185 size_t * num_src_per_tgt = xmalloc(count * sizeof(*num_src_per_tgt));
1186
1187 // compute the weights for the target points
1188 size_t num_interpolated_tgt =
1190 tgt_points, interp_fail_flag, count, super_cells, num_super_cells,
1191 &src_per_tgt, &w, num_src_per_tgt);
1192 if (num_src_cells > 0) free(src_cell_gradients->data);
1193 free(src_cell_gradients);
1194 free(super_cells);
1195 free(interp_fail_flag);
1196
1197 size_t total_num_weights = 0;
1198 for (size_t i = 0; i < num_interpolated_tgt; ++i)
1199 total_num_weights += num_src_per_tgt[i];
1200
1201 struct remote_points tgts = {
1202 .data =
1204 interp_grid, tgt_points, num_interpolated_tgt),
1205 .count = num_interpolated_tgt};
1206 struct remote_point * srcs =
1208 interp_grid, 0, src_per_tgt, total_num_weights);
1209
1210 // store weights
1212 weights, &tgts, num_src_per_tgt, srcs, w);
1213
1214 free(tgts.data);
1215 free(srcs);
1216 free(w);
1217 free(num_src_per_tgt);
1218 free(src_per_tgt);
1219
1220 return num_interpolated_tgt;
1221}
1222
1224 int order, int enforced_conserv, int partial_coverage,
1225 enum yac_interp_method_conserv_normalisation normalisation) {
1226
1227 struct interp_method_conserv * method = xmalloc(1 * sizeof(*method));
1228
1229 YAC_ASSERT(
1230 (enforced_conserv != 1) || (order == 1),
1231 "ERROR(yac_interp_method_conserv_new): interp_method_conserv only "
1232 "supports enforced_conserv with first order conservative remapping")
1233 YAC_ASSERT(
1234 (order == 1) || (order == 2),
1235 "ERROR(yac_interp_method_conserv_new): invalid order")
1236
1237 method->vtable =
1238 (order == 1)?
1242 method->normalisation = normalisation;
1244
1245 return (struct interp_method*)method;
1246}
1247
1248static void delete_conserv(struct interp_method * method) {
1249 free(method);
1250}
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:396
Structs and interfaces for area calculations.
#define YAC_AREA_TOL
Definition area.h:23
static int edge_direction(double *a, double *b)
void yac_correct_weights(size_t nSourceCells, double *weight)
correct interpolation weights
Definition clipping.c:1437
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:109
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:266
void yac_const_basic_grid_data_get_grid_cell(struct yac_const_basic_grid_data *grid_data, size_t cell_idx, struct yac_grid_cell *buffer_cell)
Definition dist_grid.c:2998
int const * const_int_pointer
static void crossproduct_kahan(double const a[], double const b[], double cross[])
Definition geometry.h:332
static void normalise_vector(double v[])
Definition geometry.h:689
yac_edge_type
Definition grid_cell.h:12
@ YAC_GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:13
void yac_interp_grid_do_cell_search_tgt(struct yac_interp_grid *interp_grid, size_t *src_cells, size_t count, size_t **tgt_cells, size_t *num_tgt_per_src)
enum yac_location yac_interp_grid_get_tgt_field_location(struct yac_interp_grid *interp_grid)
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_src_cell_neighbours(struct yac_interp_grid *interp_grid, size_t *src_cells, size_t count, size_t *neighbours)
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)
struct yac_const_basic_grid_data * yac_interp_grid_get_basic_grid_data_tgt(struct yac_interp_grid *interp_grid)
void yac_interp_grid_do_cell_search_src(struct yac_interp_grid *interp_grid, size_t *tgt_cells, size_t count, size_t **src_cells, size_t *num_src_per_tgt)
struct yac_const_basic_grid_data * yac_interp_grid_get_basic_grid_data_src(struct yac_interp_grid *interp_grid)
static void delete_conserv(struct interp_method *method)
static void get_cell_buffers_(struct yac_interp_grid *interp_grid, struct yac_grid_cell *tgt_grid_cell, struct yac_grid_cell *src_grid_cell)
static size_t compute_2nd_order_weights(size_t *tgt_cells, int *interp_fail_flag, size_t num_tgt_cells, struct supermesh_cell *super_cells, size_t num_super_cells, size_t **src_per_tgt, double **weights, size_t *num_src_per_tgt)
static yac_coordinate_pointer compute_src_cell_centroids(struct yac_interp_grid *interp_grid, size_t *src_cells, int *skip_src_cell, size_t num_src_cells, struct supermesh_cell *super_cells, size_t num_super_cells)
static void get_cell_buffers(struct yac_interp_grid *interp_grid, size_t max_num_src_per_tgt, struct yac_grid_cell *tgt_grid_cell, struct yac_grid_cell **src_grid_cells)
#define AREA_TOL_FACTOR
static void compute_cell_barycenter(struct yac_const_basic_grid_data *grid_data, size_t cell_idx, double barycenter[3])
static size_t do_search_conserv_1st_order(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static void compact_weight_vector_data(struct weight_vector_data *weights, size_t *n)
static struct interp_method_vtable interp_method_conserv_2nd_order_vtable
struct interp_method * yac_interp_method_conserv_new(int order, int enforced_conserv, int partial_coverage, enum yac_interp_method_conserv_normalisation normalisation)
static int get_max_num_vertices_per_cell(struct yac_const_basic_grid_data *basic_grid_data)
static void compute_super_cells(struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct supermesh_cell **super_cells_, size_t *num_super_cells, int *interp_fail_flag, size_t **src_cells, size_t *num_src_cells, enum yac_interp_method_conserv_normalisation normalisation, int partial_coverage)
static int compare_supermesh_cell_tgt_local_ids(const void *a, const void *b)
static void orthogonalise_weight_vector(double *src_cell_centroid, struct weight_vector_3d *G_i, struct weight_vector_data_3d *buffer)
static struct interp_method_vtable interp_method_conserv_1st_order_vtable
static int compare_supermesh_cell_src_local_ids(const void *a, const void *b)
static size_t compute_2nd_order_tgt_cell_weights(struct supermesh_cell *super_cell, struct weight_vector_data *weights)
static int compare_weight_vector_data_weight(void const *a, void const *b)
static int compare_weight_vector_data(void const *a, void const *b)
static size_t do_search_conserv_2nd_order(struct interp_method *method, struct yac_interp_grid *interp_grid, size_t *tgt_points, size_t count, struct yac_interp_weights *weights)
static struct weight_vector_3d * compute_src_cell_gradients(struct yac_interp_grid *interp_grid, size_t *src_cells, yac_coordinate_pointer src_cell_centroids, int *skip_src_cell, size_t num_src_cells, size_t *src_cell_neighbours)
static int compute_1st_order_weights(struct yac_const_basic_grid_data *tgt_basic_grid_data, size_t tgt_cell, struct yac_const_basic_grid_data *src_basic_grid_data, size_t src_count, size_t *src_cells, struct yac_grid_cell tgt_grid_cell_buffer, struct yac_grid_cell *src_grid_cell_buffer, double *weights, size_t *num_weights, int partial_coverage, enum yac_interp_method_conserv_normalisation normalisation, int enforced_conserv)
yac_interp_method_conserv_normalisation
@ YAC_INTERP_CONSERV_DESTAREA
@ YAC_INTERP_CONSERV_FRACAREA
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_LOC_CELL
Definition location.h:14
Definition __init__.py:1
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xmalloc(size)
Definition ppm_xfuncs.h:66
enum yac_interp_method_conserv_normalisation normalisation
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
struct weight_vector_3d * src_cell_gradient
struct supermesh_cell::@10 tgt
struct supermesh_cell::@10 src
struct weight_vector_data_3d * data
const yac_const_coordinate_pointer vertex_coordinates
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]
enum yac_edge_type * edge_type
Definition grid_cell.h:20
size_t array_size
Definition grid_cell.h:22
double(* coordinates_xyz)[3]
Definition grid_cell.h:19
#define MAX(a, b)
void yac_quicksort_index_yac_int_size_t(yac_int *a, size_t n, size_t *idx)
void yac_quicksort_index_int_size_t(int *a, size_t n, size_t *idx)
static void yac_remove_duplicates_size_t(size_t *array, size_t *n)
Definition utils_core.h:108
void yac_quicksort_index_size_t_int(size_t *a, size_t n, int *idx)
#define YAC_ASSERT(exp, msg)
Definition yac_assert.h:16
Xt_int yac_int
Definition yac_types.h:15
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19