YetAnotherCoupler 3.4.0
Loading...
Searching...
No Matches
debug_grid2vtk.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 <unistd.h>
8#include <string.h>
9#include <math.h>
10
11#include <netcdf.h>
12
13#include "yac_utils.h"
14
15// redefine YAC assert macros
16#undef YAC_ASSERT
17#undef YAC_ASSERT_F
18
19#define YAC_RAD (0.01745329251994329576923690768489) // M_PI / 180
20
21static char const * cmd;
22#define STR_USAGE "Usage: %s -n grid_name -f grid_filename -o vtk_filename\n"
23
24#define YAC_ASSERT(exp, msg) \
25 { \
26 if(!((exp))) { \
27 fprintf(stderr, "ERROR: %s\n" STR_USAGE, msg, cmd); \
28 exit(EXIT_FAILURE); \
29 } \
30 }
31
32#define YAC_ASSERT_F(exp, format, ...) \
33 { \
34 if(!((exp))) { \
35 fprintf(stderr, "ERROR: " format "\n" STR_USAGE, __VA_ARGS__, cmd); \
36 exit(EXIT_FAILURE); \
37 } \
38 }
39
40static void parse_arguments(int argc, char ** argv,
41 char const ** grid_filename,
42 char const ** grid_name,
43 char const ** vtk_filename);
44static void read_grid_file(
45 char const * grid_filename, char const * grid_name,
46 size_t * nv, size_t * nc,
47 double ** clon, double ** clat, double ** vlon, double ** vlat,
48 int ** gid, int ** cmk, int ** rnk);
49static void write_vtk_file(
50 char const * vtk_filename, char const * grid_name,
51 size_t nv, size_t nc,
52 double * clon, double * clat, double * vlon, double * vlat,
53 int * gid, int * cmk, int * rnk);
54
55int main(int argc, char ** argv) {
56
57 cmd = argv[0];
58
59 char const * grid_name, * grid_filename, * vtk_filename;
60 parse_arguments(argc, argv, &grid_filename, &grid_name, &vtk_filename);
61
62 size_t nv, nc;
63 double * clon, * clat, * lon, * lat;
64 int * gid, * cmk, * rnk;
65
67 grid_filename, grid_name,
68 &nv, &nc, &clon, &clat, &lon, &lat, &gid, &cmk, &rnk);
69
71 vtk_filename, grid_name,
72 nv, nc, clon, clat, lon, lat, gid, cmk, rnk);
73
74 free(clon);
75 free(clat);
76 free(lon);
77 free(lat);
78 free(gid);
79 free(cmk);
80 free(rnk);
81}
82
83static double * read_coords(int ncid, int varid, size_t varlen) {
84
85 int is_degree = yac_check_coord_units(ncid, varid);
86
87 double * coord = malloc(varlen * sizeof(*coord));
88 YAC_HANDLE_ERROR(nc_get_var_double (ncid, varid, coord));
89
90 // convert to radiant if necessary
91 if (is_degree) for (size_t i = 0; i < varlen; ++i) coord[i] *= YAC_RAD;
92
93 return coord;
94}
95
96static int * read_int(int ncid, int varid, size_t varlen) {
97
98 int * array = malloc(varlen * sizeof(*array));
99 YAC_HANDLE_ERROR(nc_get_var_int(ncid, varid, array));
100 return array;
101}
102
103static void read_grid_file(
104 char const * grid_filename, char const * grid_name,
105 size_t * nv, size_t * nc,
106 double ** clon, double ** clat, double ** lon, double ** lat,
107 int ** gid, int ** cmk, int ** rnk) {
108
109 int status;
110
111 size_t grid_name_len = strlen(grid_name) + 1;
112 char nv_dim_name[3 + grid_name_len];
113 char nc_dim_name[3 + grid_name_len];
114 char cla_var_name[4 + grid_name_len];
115 char clo_var_name[4 + grid_name_len];
116 char lat_var_name[4 + grid_name_len];
117 char lon_var_name[4 + grid_name_len];
118 char gid_var_name[4 + grid_name_len];
119 char cmk_var_name[4 + grid_name_len];
120 char rnk_var_name[4 + grid_name_len];
121
122 snprintf(nv_dim_name, 3 + grid_name_len, "nv_%s", grid_name);
123 snprintf(nc_dim_name, 3 + grid_name_len, "nc_%s", grid_name);
124 snprintf(cla_var_name, 4 + grid_name_len, "%s.cla", grid_name);
125 snprintf(clo_var_name, 4 + grid_name_len, "%s.clo", grid_name);
126 snprintf(lat_var_name, 4 + grid_name_len, "%s.lat", grid_name);
127 snprintf(lon_var_name, 4 + grid_name_len, "%s.lon", grid_name);
128 snprintf(gid_var_name, 4 + grid_name_len, "%s.gid", grid_name);
129 snprintf(cmk_var_name, 4 + grid_name_len, "%s.cmk", grid_name);
130 snprintf(rnk_var_name, 4 + grid_name_len, "%s.rnk", grid_name);
131
133 yac_file_exists(grid_filename), "File %s does not exist.", grid_filename)
134
135 int ncid, dimid, var_id;
136 yac_nc_open(grid_filename, NC_NOWRITE, &ncid);
137
138 yac_nc_inq_dimid(ncid, nv_dim_name, &dimid);
139 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimid, nv));
140 yac_nc_inq_dimid(ncid, nc_dim_name, &dimid);
141 YAC_HANDLE_ERROR(nc_inq_dimlen(ncid, dimid, nc));
142
143 yac_nc_inq_varid(ncid, cla_var_name, &var_id);
144 *clat = read_coords(ncid, var_id, *nv * *nc);
145 yac_nc_inq_varid(ncid, clo_var_name, &var_id);
146 *clon = read_coords(ncid, var_id, *nv * *nc);
147
148 status = nc_inq_varid(ncid, lat_var_name, &var_id);
149 if (status == NC_NOERR) {
150 *lat = read_coords(ncid, var_id, *nc);
151 } else if (status == NC_ENOTVAR) {
152 *lat = NULL;
153 status = NC_NOERR;
154 }
155 YAC_HANDLE_ERROR(status);
156 status = nc_inq_varid(ncid, lon_var_name, &var_id);
157 if (status == NC_NOERR) {
158 *lon = read_coords(ncid, var_id, *nc);
159 } else if (status == NC_ENOTVAR) {
160 *lon = NULL;
161 status = NC_NOERR;
162 }
163 YAC_HANDLE_ERROR(status);
164
165 status = nc_inq_varid(ncid, gid_var_name, &var_id);
166 if (status == NC_NOERR) {
167 *gid = read_int(ncid, var_id, *nc);
168 } else if (status == NC_ENOTVAR) {
169 *gid = NULL;
170 status = NC_NOERR;
171 }
172 YAC_HANDLE_ERROR(status);
173
174 status = nc_inq_varid(ncid, cmk_var_name, &var_id);
175 if (status == NC_NOERR) {
176 *cmk = read_int(ncid, var_id, *nc);
177 } else if (status == NC_ENOTVAR) {
178 *cmk = NULL;
179 status = NC_NOERR;
180 }
181 YAC_HANDLE_ERROR(status);
182
183 YAC_HANDLE_ERROR(nc_inq_varid(ncid, rnk_var_name, &var_id));
184 *rnk = read_int(ncid, var_id, *nc);
185
186 YAC_HANDLE_ERROR(nc_close(ncid));
187}
188
189static void LLtoXYZ(double lon, double lat, double p_out[]) {
190
191 if ((fabs(lon) > 8.0 * M_PI) || (fabs(lat) > 2.0 * M_PI)) {
192 p_out[0] = 0.0;
193 p_out[1] = 0.0;
194 p_out[2] = 0.0;
195 return;
196 }
197
198 while (lon < -M_PI) lon += 2.0 * M_PI;
199 while (lon >= M_PI) lon -= 2.0 * M_PI;
200
201 double cos_lat = cos(lat);
202 p_out[0] = cos_lat * cos(lon);
203 p_out[1] = cos_lat * sin(lon);
204 p_out[2] = sin(lat);
205}
206
207static void write_vtk_file(
208 char const * vtk_filename, char const * grid_name,
209 size_t nv, size_t nc,
210 double * clon, double * clat, double * lon, double * lat,
211 int * gid, int * cmk, int * rnk) {
212
214 malloc((nv * nc + ((lon && lat)?nc:0)) * sizeof(*points));
215
216 for (size_t i = 0; i < nv * nc; ++i) LLtoXYZ(clon[i], clat[i], points[i]);
217 if (lon && lat)
218 for (size_t i = 0; i < nc; ++i)
219 LLtoXYZ(lon[i], lat[i], points[nv * nc + i]);
220
221 unsigned * num_points_per_polygon =
222 malloc(nc * sizeof(*num_points_per_polygon));
223 for (size_t i = 0; i < nc; ++i) num_points_per_polygon[i] = (unsigned)nv;
224
225 unsigned * polygon_data = malloc(nc * nv * sizeof(*polygon_data));
226 for (size_t i = 0; i < nv * nc; ++i) polygon_data[i] = (unsigned)i;
227
228 YAC_VTK_FILE * file = yac_vtk_open(vtk_filename, grid_name);
229 yac_vtk_write_point_data(file, &(points[0][0]), nv * nc + ((lon && lat)?nc:0));
231 file, polygon_data, num_points_per_polygon, nc);
232 if (gid) yac_vtk_write_cell_scalars_int(file, gid, nc, "gid");
233 if (cmk) yac_vtk_write_cell_scalars_int(file, cmk, nc, "cmk");
234 if (rnk) yac_vtk_write_cell_scalars_int(file, rnk, nc, "rnk");
235 yac_vtk_close(file);
236
237 free(polygon_data);
238 free(num_points_per_polygon);
239 free(points);
240}
241
242static void parse_arguments(int argc, char ** argv,
243 char const ** grid_filename,
244 char const ** grid_name,
245 char const ** vtk_filename) {
246
247 *grid_name = NULL;
248 *grid_filename = NULL;
249 *vtk_filename = NULL;;
250
251 int opt;
252 while ((opt = getopt(argc, argv, "n:f:o:")) != -1) {
254 (opt == 'n') ||
255 (opt == 'f') ||
256 (opt == 'o'), "invalid command argument")
257 switch (opt) {
258 default:
259 case 'n':
260 *grid_name = optarg;
261 break;
262 case 'f':
263 *grid_filename = optarg;
264 break;
265 case 'o':
266 *vtk_filename = optarg;
267 break;
268 }
269 }
270 YAC_ASSERT_F(optind >= argc, "non-option ARGV-element: \"%s\"", argv[optind])
271 YAC_ASSERT(argc != 1, "too few arguments")
272 YAC_ASSERT(*grid_name != NULL, "grid_name argument is missing")
273 YAC_ASSERT(*grid_filename != NULL, "grid_filename argument is missing")
274 YAC_ASSERT(*vtk_filename != NULL, "vtk_filename argument is missing")
275}
static void parse_arguments(int argc, char **argv, char const **grid_filename, char const **grid_name, char const **vtk_filename)
static char const * cmd
int main(int argc, char **argv)
static void write_vtk_file(char const *vtk_filename, char const *grid_name, size_t nv, size_t nc, double *clon, double *clat, double *vlon, double *vlat, int *gid, int *cmk, int *rnk)
#define YAC_ASSERT_F(exp, format,...)
static double * read_coords(int ncid, int varid, size_t varlen)
static int * read_int(int ncid, int varid, size_t varlen)
#define YAC_RAD
static void LLtoXYZ(double lon, double lat, double p_out[])
static void read_grid_file(char const *grid_filename, char const *grid_name, size_t *nv, size_t *nc, double **clon, double **clat, double **vlon, double **vlat, int **gid, int **cmk, int **rnk)
#define YAC_ASSERT(exp, msg)
void yac_nc_inq_varid(int ncid, char const *name, int *varidp)
Definition io_utils.c:405
void yac_nc_open(const char *path, int omode, int *ncidp)
Definition io_utils.c:344
void yac_nc_inq_dimid(int ncid, char const *name, int *dimidp)
Definition io_utils.c:379
int yac_file_exists(const char *filename)
Definition utils_core.c:12
#define YAC_HANDLE_ERROR(exp)
Definition io_utils.h:30
int yac_check_coord_units(int ncid, int varid)
Definition read_grid.c:45
void yac_vtk_write_cell_scalars_int(YAC_VTK_FILE *vtk_file, int *scalars, unsigned num_cells, char *name)
Definition vtk_output.c:250
void yac_vtk_write_point_data(YAC_VTK_FILE *vtk_file, double *point_data, unsigned num_points)
Definition vtk_output.c:69
void yac_vtk_close(YAC_VTK_FILE *vtk_file)
Definition vtk_output.c:403
void yac_vtk_write_cell_data(YAC_VTK_FILE *vtk_file, unsigned *cell_corners, unsigned *num_points_per_cell, unsigned num_cells)
Definition vtk_output.c:88
YAC_VTK_FILE * yac_vtk_open(const char *filename, const char *title)
Definition vtk_output.c:45
static struct user_input_data_points ** points
Definition yac.c:136
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19