YetAnotherCoupler 3.1.1
Loading...
Searching...
No Matches
grid_unstruct.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 <string.h>
6
7#include "basic_grid_data.h"
8#include "geometry.h"
9#include "utils_common.h"
10
15
16static int compare_temp_edges(void const * a, void const * b) {
17
18 struct temp_edge const * edge_a = (struct temp_edge const *)a;
19 struct temp_edge const * edge_b = (struct temp_edge const *)b;
20
21 if (edge_a->vertex[0] != edge_b->vertex[0])
22 return (edge_a->vertex[0] > edge_b->vertex[0])?1:-1;
23 return (edge_a->vertex[1] > edge_b->vertex[1]) -
24 (edge_a->vertex[1] < edge_b->vertex[1]);
25}
26
28 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
29 double *x_vertices, double *y_vertices, int *cell_to_vertex_,
30 void (*LLtoXYZ_ptr)(double, double, double[])) {
31
33 xmalloc(nbr_vertices * sizeof(*vertex_coordinates));
34 for (size_t i = 0; i < nbr_vertices; ++i)
35 LLtoXYZ_ptr(
36 x_vertices[i], y_vertices[i], &(vertex_coordinates[i][0]));
37
39 xmalloc(nbr_cells * sizeof(*num_vertices_per_cell));
40 memcpy(num_vertices_per_cell, num_vertices_per_cell_,
41 nbr_cells * sizeof(*num_vertices_per_cell));
42
43 size_t total_num_cell_corners = 0;
44 size_t * cell_to_vertex_offsets =
45 xmalloc(nbr_cells * sizeof(*cell_to_vertex_offsets));
46 for (size_t i = 0; i < nbr_cells; ++i) {
47 cell_to_vertex_offsets[i] = total_num_cell_corners;
48 total_num_cell_corners += (size_t)(num_vertices_per_cell[i]);
49 }
50
51 for (size_t i = 0; i < total_num_cell_corners; ++i){
52 YAC_ASSERT_F(cell_to_vertex_[i] >= 0,
53 "ERROR(yac_generate_basic_grid_data_unstruct_): "
54 "Index %zu in cell_to_vertex (%d) is negative.",
55 i, cell_to_vertex_[i])
56 YAC_ASSERT_F(cell_to_vertex_[i] < nbr_vertices,
57 "ERROR(yac_generate_basic_grid_data_unstruct_): "
58 "Index %zu in cell_to_vertex (%d) is not smaller than nbr_vertices (%zu).",
59 i, cell_to_vertex_[i], nbr_vertices)
60 }
61
62 size_t * cell_to_vertex =
63 xmalloc(total_num_cell_corners * sizeof(*cell_to_vertex));
64 for (size_t i = 0; i < total_num_cell_corners; ++i)
65 cell_to_vertex[i] = (size_t)(cell_to_vertex_[i]);
66
68 xcalloc(nbr_vertices, sizeof(*num_cells_per_vertex));
69 for (size_t i = 0; i < total_num_cell_corners; ++i)
71
72 size_t * vertex_to_cell =
73 xmalloc(total_num_cell_corners * sizeof(*vertex_to_cell));
74 size_t * vertex_to_cell_offsets =
75 xmalloc((nbr_vertices + 1) * sizeof(*vertex_to_cell_offsets));
77 for (size_t i = 0, accu = 0; i < nbr_vertices; ++i) {
78 vertex_to_cell_offsets[i+1] = accu;
79 accu += num_cells_per_vertex[i];
80 }
81 for (size_t i = 0, k = 0; i < nbr_cells; ++i) {
82 size_t curr_num_vertices = num_vertices_per_cell[i];
83 for (size_t j = 0; j < curr_num_vertices; ++j, ++k)
85 }
86
87 struct temp_edge * temp_edges =
88 xmalloc(total_num_cell_corners * sizeof(*temp_edges));
89 for (size_t i = 0, offset = 0, k = 0; i < nbr_cells; ++i) {
90 size_t * curr_cell_to_vertex = cell_to_vertex + offset;
91 size_t curr_num_edges = num_vertices_per_cell[i];
92 offset += curr_num_edges;
93 for (size_t j = 0; j < curr_num_edges; ++j, ++k) {
94 int order =
95 curr_cell_to_vertex[j] > curr_cell_to_vertex[(j+1)%curr_num_edges];
96 temp_edges[k].vertex[order] = curr_cell_to_vertex[j];
97 temp_edges[k].vertex[order^1] = curr_cell_to_vertex[(j+1)%curr_num_edges];
98 temp_edges[k].cell_to_edge_idx = k;
99 }
100 }
101 qsort(temp_edges, total_num_cell_corners,
102 sizeof(*temp_edges), compare_temp_edges);
103
104 size_t * cell_to_edge =
105 xmalloc(total_num_cell_corners * sizeof(*cell_to_edge));
106 size_t nbr_edges = 0;
107 yac_size_t_2_pointer edge_to_vertex = (yac_size_t_2_pointer)temp_edges;
108 for (size_t i = 0, prev_indices[2] = {SIZE_MAX, SIZE_MAX};
109 i < total_num_cell_corners; ++i) {
110
111 size_t curr_cell_to_edge_idx = temp_edges[i].cell_to_edge_idx;
112 if ((prev_indices[0] != temp_edges[i].vertex[0]) ||
113 (prev_indices[1] != temp_edges[i].vertex[1])) {
114
115 prev_indices[0] = temp_edges[i].vertex[0];
116 prev_indices[1] = temp_edges[i].vertex[1];
117 edge_to_vertex[nbr_edges][0] = prev_indices[0];
118 edge_to_vertex[nbr_edges][1] = prev_indices[1];
119 ++nbr_edges;
120 }
121
122 cell_to_edge[curr_cell_to_edge_idx] = nbr_edges - 1;
123 }
124 edge_to_vertex =
125 xrealloc(edge_to_vertex, nbr_edges * sizeof(*edge_to_vertex));
126
127 enum yac_edge_type * edge_type = xmalloc(nbr_edges * sizeof(*edge_type));
128 for (size_t i = 0; i < nbr_edges; ++i) edge_type[i] = GREAT_CIRCLE_EDGE;
129
130 struct yac_basic_grid_data grid;
132 grid.cell_ids = NULL;
133 grid.vertex_ids = NULL;
134 grid.edge_ids = NULL;
135 grid.num_cells = nbr_cells;
136 grid.num_vertices = nbr_vertices;
137 grid.num_edges = nbr_edges;
138 grid.core_cell_mask = NULL;
139 grid.core_vertex_mask = NULL;
140 grid.core_edge_mask = NULL;
150 grid.edge_type = edge_type;
151 grid.num_total_cells = nbr_cells;
152 grid.num_total_vertices = nbr_vertices;
153 grid.num_total_edges = nbr_edges;
154 return grid;
155}
156
158 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
159 double *x_vertices, double *y_vertices, int *cell_to_vertex_) {
160
161 return
163 nbr_vertices, nbr_cells, num_vertices_per_cell_,
164 x_vertices, y_vertices, cell_to_vertex_, LLtoXYZ);
165}
166
168 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
169 double *x_vertices, double *y_vertices, int *cell_to_vertex_) {
170
171 return
173 nbr_vertices, nbr_cells, num_vertices_per_cell_,
174 x_vertices, y_vertices, cell_to_vertex_, LLtoXYZ_deg);
175}
176
178 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
179 double *x_vertices, double *y_vertices, int *cell_to_vertex_,
180 void (*LLtoXYZ_ptr)(double, double, double[]),
181 double (*get_angle_ptr)(double, double),
182 double angle_tol, double pole_y_vertex) {
183
184 struct yac_basic_grid_data grid_data =
186 nbr_vertices, nbr_cells, num_vertices_per_cell_,
187 x_vertices, y_vertices, cell_to_vertex_, LLtoXYZ_ptr);
188
189 for (size_t i = 0; i < grid_data.num_edges; ++i) {
190 int is_lon_edge =
191 (fabs(
192 get_angle_ptr(
193 x_vertices[grid_data.edge_to_vertex[i][0]],
194 x_vertices[grid_data.edge_to_vertex[i][1]])) < angle_tol) ||
195 ((fabs(fabs(y_vertices[grid_data.edge_to_vertex[i][0]]) -
196 pole_y_vertex) < angle_tol) ^
197 (fabs(fabs(y_vertices[grid_data.edge_to_vertex[i][1]]) -
198 pole_y_vertex) < angle_tol));
199 int is_lat_edge =
200 fabs(y_vertices[grid_data.edge_to_vertex[i][0]] -
201 y_vertices[grid_data.edge_to_vertex[i][1]]) < angle_tol;
203 is_lon_edge || is_lat_edge,
204 "ERROR(yac_generate_basic_grid_data_unstruct_ll): "
205 "edge is neither lon nor lat ((%lf,%lf),(%lf,%lf))",
206 x_vertices[grid_data.edge_to_vertex[i][0]],
207 y_vertices[grid_data.edge_to_vertex[i][0]],
208 x_vertices[grid_data.edge_to_vertex[i][1]],
209 y_vertices[grid_data.edge_to_vertex[i][1]]);
210 grid_data.edge_type[i] = (is_lon_edge)?LON_CIRCLE_EDGE:LAT_CIRCLE_EDGE;
211 }
212
213 return grid_data;
214}
215
217 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
218 double *x_vertices, double *y_vertices, int *cell_to_vertex_) {
219
220 return
222 nbr_vertices, nbr_cells, num_vertices_per_cell_,
223 x_vertices, y_vertices, cell_to_vertex_, LLtoXYZ,
224 get_angle, yac_angle_tol, M_PI_2);
225}
226
228 size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_,
229 double *x_vertices, double *y_vertices, int *cell_to_vertex_) {
230
231 return
233 nbr_vertices, nbr_cells, num_vertices_per_cell_,
234 x_vertices, y_vertices, cell_to_vertex_, LLtoXYZ_deg,
235 get_angle_deg, yac_angle_tol / M_PI * 180.0, 90.0);
236}
static void LLtoXYZ_deg(double lon, double lat, double p_out[])
Definition geometry.h:304
#define yac_angle_tol
Definition geometry.h:34
static double get_angle(double a_lon, double b_lon)
Definition geometry.h:127
static void LLtoXYZ(double lon, double lat, double p_out[])
Definition geometry.h:287
static double get_angle_deg(double a_lon, double b_lon)
Definition geometry.h:143
yac_edge_type
Definition grid_cell.h:10
@ LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:12
@ GREAT_CIRCLE_EDGE
great circle
Definition grid_cell.h:11
@ LON_CIRCLE_EDGE
longitude circle
Definition grid_cell.h:13
static struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_ll_(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_, void(*LLtoXYZ_ptr)(double, double, double[]), double(*get_angle_ptr)(double, double), double angle_tol, double pole_y_vertex)
static struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_, void(*LLtoXYZ_ptr)(double, double, double[]))
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_deg(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_)
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_ll_deg(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_)
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_)
static int compare_temp_edges(void const *a, void const *b)
struct yac_basic_grid_data yac_generate_basic_grid_data_unstruct_ll(size_t nbr_vertices, size_t nbr_cells, int *num_vertices_per_cell_, double *x_vertices, double *y_vertices, int *cell_to_vertex_)
#define xrealloc(ptr, size)
Definition ppm_xfuncs.h:67
#define xcalloc(nmemb, size)
Definition ppm_xfuncs.h:64
#define xmalloc(size)
Definition ppm_xfuncs.h:66
size_t cell_to_edge_idx
yac_int vertex[2]
yac_coordinate_pointer vertex_coordinates
size_t * vertex_to_cell_offsets
yac_size_t_2_pointer edge_to_vertex
enum yac_edge_type * edge_type
size_t * cell_to_vertex_offsets
#define YAC_ASSERT_F(exp, format,...)
Xt_int yac_int
Definition yac_types.h:11
size_t(* yac_size_t_2_pointer)[2]
Definition yac_types.h:19
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:15