YetAnotherCoupler 3.2.0
Loading...
Searching...
No Matches
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 <stdio.h>
6#include <string.h>
7#include <float.h>
8
9#include "grid2vtk.h"
10#include "vtk_output.h"
11#include "geometry.h"
12#include "ensure_array_size.h"
13#include "utils_common.h"
14
15static void
17 size_t num_points, double * points) {
18
19 double n_x[3] = {1, 0, 0};
20 double n_y[3] = {0, 1, 0};
21 double angle_x = get_vector_angle(n_x, bnd_circle.base_vector);
22 double angle_y = get_vector_angle(n_y, bnd_circle.base_vector);
23 double * n = (angle_x > angle_y)?n_x:n_y;
24
25 double temp_vector[3];
26 double start_vector[3];
27
28 crossproduct_d(bnd_circle.base_vector, n, temp_vector);
30 temp_vector, bnd_circle.inc_angle, bnd_circle.base_vector, start_vector);
31
32 double angle = 2.0 * M_PI / (num_points + 1);
33
34 for (size_t i = 0; i < num_points; ++i)
36 bnd_circle.base_vector, angle * (double)i, start_vector, points + 3 * i);
37}
38
39static void
41
42 // the equator is approximated by a polygon with num_points corners
43
44 struct bounding_circle equator;
45 equator.base_vector[0] = 0;
46 equator.base_vector[1] = 0;
47 equator.base_vector[2] = 1;
48 equator.inc_angle = SIN_COS_M_PI_2;
49 equator.sq_crd = DBL_MAX;
50
52}
53
55 struct yac_basic_grid_data * grid, char const * name) {
56
57 size_t num_grid_corners = grid->num_vertices;
58 size_t num_grid_cells = grid->num_cells;
59
60 size_t const num_points_per_great_circle = 71;
61
62 //-------------------------------------
63 // generate point data for the vtk file
64 // it will contain data for all corners of the grid and the points for the equator
65 //-------------------------------------
67 xmalloc((num_grid_corners + num_points_per_great_circle) *sizeof(*points));
68
69 // get the point data of the grid
70 memcpy(points, grid->vertex_coordinates, num_grid_corners * sizeof(*points));
71 // get the point data of the equator
73 num_points_per_great_circle, &(points[num_grid_corners][0]));
74
75 //------------------------------------
76 // generate cell data for the vtk file
77 //------------------------------------
78
79 // has entries for all cells of the grid and one for the equator
80 unsigned * num_points_per_cell =
81 xmalloc((grid->num_cells + 1) * sizeof(*num_points_per_cell));
82
83 // get the number of points per grid cell
84 for (size_t i = 0; i < grid->num_cells; ++i)
85 num_points_per_cell[i] = (unsigned)(grid->num_vertices_per_cell[i]);
86
87 // get the number of corners for the equator
88 num_points_per_cell[grid->num_cells] = (unsigned)num_points_per_great_circle;
89
90 // count the total number of cell points
91 size_t total_num_cell_points = 0;
92 for (size_t i = 0; i < grid->num_cells; ++i)
93 total_num_cell_points += num_points_per_cell[i];
94
95 unsigned * cell_data =
96 xmalloc((total_num_cell_points + num_points_per_great_circle) *
97 sizeof(*cell_data));
98
99 // get the cell data of the grid
100 for (size_t i = 0; i < total_num_cell_points; ++i)
101 cell_data[i] = (unsigned)(grid->cell_to_vertex[i]);
102
103 // write the cell data for the great circle
104 for (size_t i = 0; i < num_points_per_great_circle; ++i)
105 cell_data[i + total_num_cell_points] = (unsigned)(num_grid_corners + i);
106
107 //------------------------------------
108 // generate scalar data
109 // we write a cell type that differentiates grid cells and the equator
110 //------------------------------------
111 unsigned * cell_type = xmalloc((num_grid_cells + 1) * sizeof(*cell_type));
112 int * core_cell_mask = xmalloc((num_grid_cells + 1) * sizeof(*core_cell_mask));
113
114 for (size_t i = 0; i < num_grid_cells; ++i) cell_type[i] = 0;
115 cell_type[num_grid_cells] = 1;
116
117 if (grid->core_cell_mask != NULL) {
118 memcpy(core_cell_mask, grid->core_cell_mask,
119 num_grid_cells * sizeof(*core_cell_mask));
120 } else {
121 for (size_t i = 0; i < num_grid_cells; ++i) core_cell_mask[i] = 1;
122 }
123 core_cell_mask[num_grid_cells] = 0;
124
125 //------------------------------------
126 // generate the file name
127 //------------------------------------
128
129 char * filename = xcalloc((strlen(name) + 5), sizeof(*filename));
130 strcpy(filename, name);
131 strcat(filename, ".vtk");
132
133 //------------------------------------
134 // generate the actual vtk file
135 //------------------------------------
136 YAC_VTK_FILE * file;
137
138 file = yac_vtk_open(filename, "grid data");
140 file, &(points[0][0]), num_grid_corners + num_points_per_great_circle);
141 yac_vtk_write_cell_data(file, cell_data, num_points_per_cell, num_grid_cells + 1);
142 yac_vtk_write_cell_scalars_uint(file, cell_type, num_grid_cells + 1, "cell_type");
143 yac_vtk_write_cell_scalars_int(file, core_cell_mask, num_grid_cells + 1, "core_cell_mask");
144
145 yac_vtk_close(file);
146
147 //------------------------------------
148 // some final cleanup
149 //------------------------------------
150
151 free(filename);
152 free(core_cell_mask);
153 free(cell_type);
154 free(cell_data);
155 free(num_points_per_cell);
156 free(points);
157}
158
159static void get_edge_points(
160 double * a, double * b, enum yac_edge_type edge_type, double (**points)[3],
161 size_t * points_array_size, size_t * num_points, size_t num_points_per_edge) {
162
163 double edge_length = get_vector_angle(a, b);
164
165 if (edge_length <= yac_angle_tol) {
166
168 *points, *points_array_size, *num_points + 1);
169
170 (*points)[*num_points][0] = a[0];
171 (*points)[*num_points][1] = a[1];
172 (*points)[*num_points][2] = a[2];
173
174 *num_points += 1;
175
176 return;
177 }
178
180 *points, *points_array_size, *num_points + num_points_per_edge);
181
182 if (edge_type == YAC_LAT_CIRCLE_EDGE) {
183
184 double a_lon, b_lon, lat;
185 XYZtoLL(a, &a_lon, &lat);
186 XYZtoLL(b, &b_lon, &lat);
187
188 double d_lon_angle = get_angle(b_lon, a_lon) / (double)num_points_per_edge;
189
190 for (size_t i = 0, offset = *num_points; i < num_points_per_edge;
191 ++i, ++offset)
192 LLtoXYZ(a_lon + d_lon_angle * (double)i, lat, (*points)[offset]);
193
194 } else {
195
196 double rotation_axis[3];
197 crossproduct_d(a, b, rotation_axis);
198 normalise_vector(rotation_axis);
199
200 double d_angle = edge_length / (double)num_points_per_edge;
201
202 for (size_t i = 0, offset = *num_points; i < num_points_per_edge;
203 ++i, ++offset)
204 rotate_vector(rotation_axis, d_angle * (double)i, a, (*points)[offset]);
205 }
206
207 *num_points += num_points_per_edge;
208}
209
211 struct yac_grid_cell * cells, size_t num_cells, char * name,
212 size_t num_points_per_edge) {
213
214 if (num_cells == 0) return;
215
216 num_points_per_edge = MAX(num_points_per_edge, 3);
217
218 unsigned * num_points_per_cell =
219 xmalloc(num_cells * sizeof(*num_points_per_cell));
220
221 //-------------------------------------
222 // get the points for all cells
223 //-------------------------------------
224
225 double (*points)[3] = NULL;
226 size_t points_array_size = 0;
227 size_t num_points = 0;
228
229 for (size_t i = 0; i < num_cells; ++i) {
230
231 struct yac_grid_cell * curr_cell = cells + i;
232 size_t num_edges = curr_cell->num_corners;
233
234 size_t prev_num_points = num_points;
235
236 // for all edges of the current cell
237 for (size_t j = 0; j < num_edges; ++j)
238 get_edge_points(curr_cell->coordinates_xyz[j],
239 curr_cell->coordinates_xyz[(j + 1) % num_edges],
240 curr_cell->edge_type[j], &points, &points_array_size,
241 &num_points, num_points_per_edge);
242
243 // get the number of points for the current cell
244 num_points_per_cell[i] = num_points - prev_num_points;
245 }
246
247 //------------------------------------
248 // generate cell data for the vtk file
249 //------------------------------------
250
251 unsigned * cell_data = xmalloc(num_points * sizeof(*cell_data));
252
253 // write the cell data for the great circle
254 for (size_t i = 0; i < num_points; ++i) cell_data[i] = (unsigned)i;
255
256 //------------------------------------
257 // generate the file name
258 //------------------------------------
259
260 char * filename = xcalloc((strlen(name) + 5), sizeof(*filename));
261 strcpy(filename, name);
262 strcat(filename, ".vtk");
263
264 //------------------------------------
265 // generate the actual vtk file
266 //------------------------------------
267
268 YAC_VTK_FILE * file = yac_vtk_open(filename, "grid data");
269 yac_vtk_write_point_data(file, &(points[0][0]), (unsigned)num_points);
270 yac_vtk_write_cell_data(file, cell_data, num_points_per_cell, (unsigned)num_cells);
271 yac_vtk_close(file);
272
273 //------------------------------------
274 // some final cleanup
275 //------------------------------------
276
277 free(filename);
278 free(cell_data);
279 free(num_points_per_cell);
280 free(points);
281}
282
284 struct yac_basic_grid * grid, char const * name) {
285
287 yac_basic_grid_get_data(grid), name);
288}
struct yac_basic_grid_data * yac_basic_grid_get_data(struct yac_basic_grid *grid)
Definition basic_grid.c:132
#define ENSURE_ARRAY_SIZE(arrayp, curr_array_size, req_size)
static void crossproduct_d(const double a[], const double b[], double cross[])
Definition geometry.h:347
static void rotate_vector(double axis[], double angle, double v_in[], double v_out[])
Definition geometry.h:732
static const struct sin_cos_angle SIN_COS_M_PI_2
Definition geometry.h:48
static void rotate_vector2(double axis[], struct sin_cos_angle angle, double v_in[], double v_out[])
Definition geometry.h:704
#define yac_angle_tol
Definition geometry.h:34
static void normalise_vector(double v[])
Definition geometry.h:689
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_vector_angle(double const a[3], double const b[3])
Definition geometry.h:364
static void XYZtoLL(double const p_in[], double *lon, double *lat)
Definition geometry.h:318
static void get_equator_point_data(size_t num_points, double *points)
Definition grid2vtk.c:40
static void get_bounding_circle_point_data(struct bounding_circle bnd_circle, size_t num_points, double *points)
Definition grid2vtk.c:16
static void get_edge_points(double *a, double *b, enum yac_edge_type edge_type, double(**points)[3], size_t *points_array_size, size_t *num_points, size_t num_points_per_edge)
Definition grid2vtk.c:159
void yac_write_basic_grid_data_to_file(struct yac_basic_grid_data *grid, char const *name)
Definition grid2vtk.c:54
void yac_write_basic_grid_to_file(struct yac_basic_grid *grid, char const *name)
Definition grid2vtk.c:283
void yac_write_grid_cells_to_file(struct yac_grid_cell *cells, size_t num_cells, char *name, size_t num_points_per_edge)
Definition grid2vtk.c:210
yac_edge_type
Definition grid_cell.h:12
@ YAC_LAT_CIRCLE_EDGE
latitude circle
Definition grid_cell.h:14
#define xcalloc(nmemb, size)
Definition ppm_xfuncs.h:64
#define xmalloc(size)
Definition ppm_xfuncs.h:66
struct sin_cos_angle inc_angle
angle between the middle point and the boundary of the spherical cap
Definition geometry.h:61
double base_vector[3]
Definition geometry.h:59
double sq_crd
Definition geometry.h:64
yac_coordinate_pointer vertex_coordinates
size_t num_corners
Definition grid_cell.h:21
enum yac_edge_type * edge_type
Definition grid_cell.h:20
double(* coordinates_xyz)[3]
Definition grid_cell.h:19
#define MAX(a, b)
void yac_vtk_write_cell_scalars_uint(YAC_VTK_FILE *vtk_file, unsigned *scalars, unsigned num_cells, char *name)
Definition vtk_output.c:243
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
general routines for writing vtk files
static struct user_input_data_points ** points
Definition yac.c:120
static size_t num_points
Definition yac.c:121
double(* yac_coordinate_pointer)[3]
Definition yac_types.h:19