YAC 3.7.1
Yet Another Coupler
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
generate_cubed_sphere.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 <stdlib.h>
7#include <math.h>
8
10#include "generate_reg2d.h"
11#include "geometry.h"
12#include "ppm/ppm_xfuncs.h"
13
14#ifdef YAC_NETCDF_ENABLED
15#include "io_utils.h"
16#include <netcdf.h>
17#endif
18
19// This routine is based on Mathlab code provided by Mike Hobson and written by
20// Mike Rezny (both from MetOffice)
22 unsigned n, unsigned * num_cells, unsigned * num_vertices,
23 double ** x_vertices, double ** y_vertices, double ** z_vertices,
24 unsigned ** cell_to_vertex, unsigned ** face_id) {
25
26 YAC_ASSERT((n >= 1) && (n <= 40000), "invalid number of linear subdivisions")
27
28 *num_cells = n * n * 6;
29 *num_vertices = *num_cells + 2;
30
31 // allocation of output variables
32 *x_vertices = xmalloc(*num_vertices * sizeof(**x_vertices));
33 *y_vertices = xmalloc(*num_vertices * sizeof(**y_vertices));
34 *z_vertices = xmalloc(*num_vertices * sizeof(**z_vertices));
35
37
38 *face_id = xmalloc(*num_cells * sizeof(**face_id));
39
40 // allocation of temporary coordinate variables
41 unsigned num_edge_coords = n - 1;
42 double * cube_edge_x_vertices = xmalloc(12 * num_edge_coords * sizeof(*cube_edge_x_vertices));
43 double * cube_edge_y_vertices = xmalloc(12 * num_edge_coords * sizeof(*cube_edge_y_vertices));
44 double * cube_edge_z_vertices = xmalloc(12 * num_edge_coords * sizeof(*cube_edge_z_vertices));
45
46 unsigned num_inner_coords = num_edge_coords * num_edge_coords;
47 double * cube_inner_x_vertices = xmalloc(num_inner_coords * sizeof(*cube_inner_x_vertices));
48 double * cube_inner_y_vertices = xmalloc(num_inner_coords * sizeof(*cube_inner_y_vertices));
49 double * cube_inner_z_vertices = xmalloc(num_inner_coords * sizeof(*cube_inner_z_vertices));
50
51 double * temp_x_vertices = xmalloc((n + 1) * (n + 1) * sizeof(*temp_x_vertices));
52 double * temp_y_vertices = xmalloc((n + 1) * (n + 1) * sizeof(*temp_y_vertices));
53 double * temp_z_vertices = xmalloc((n + 1) * (n + 1) * sizeof(*temp_z_vertices));
54
55 {
56 double * theta = xmalloc((n + 1) * sizeof(*theta));
57
58 {
59 double d = (M_PI_2 / (double)n);
60 for (unsigned i = 0; i < n + 1; ++i)
61 theta[i] = tan(- M_PI_4 + d * (double)i);
62 }
63
64 for (unsigned i = 0; i < n + 1; ++i) {
65 for (unsigned j = 0; j < n + 1; ++j) {
66
67 double scale =
68 sqrt(1.0 / (theta[i] * theta[i] + theta[j] * theta[j] + 1.0));
69
70 temp_x_vertices[i * (n + 1) + j] = theta[j] * scale;
71 temp_y_vertices[i * (n + 1) + j] = theta[i] * scale;
72 temp_z_vertices[i * (n + 1) + j] = - scale;
73 }
74 }
75 free(theta);
76 }
77
78 // Store the coordinates for the 4 vertices on face 1
79 (*x_vertices)[0] = temp_x_vertices[ 0 * (n + 1) + 0];
80 (*y_vertices)[0] = temp_y_vertices[ 0 * (n + 1) + 0];
81 (*z_vertices)[0] = temp_z_vertices[ 0 * (n + 1) + 0];
82 (*x_vertices)[1] = temp_x_vertices[ n * (n + 1) + 0];
83 (*y_vertices)[1] = temp_y_vertices[ n * (n + 1) + 0];
84 (*z_vertices)[1] = temp_z_vertices[ n * (n + 1) + 0];
85 (*x_vertices)[2] = temp_x_vertices[ 0 * (n + 1) + n];
86 (*y_vertices)[2] = temp_y_vertices[ 0 * (n + 1) + n];
87 (*z_vertices)[2] = temp_z_vertices[ 0 * (n + 1) + n];
88 (*x_vertices)[3] = temp_x_vertices[ n * (n + 1) + n];
89 (*y_vertices)[3] = temp_y_vertices[ n * (n + 1) + n];
90 (*z_vertices)[3] = temp_z_vertices[ n * (n + 1) + n];
91 // Store the coordinates for the 4 vertices on face 2
92 (*x_vertices)[4] = temp_x_vertices[ 0 * (n + 1) + 0];
93 (*y_vertices)[4] = temp_y_vertices[ 0 * (n + 1) + 0];
94 (*z_vertices)[4] = -temp_z_vertices[ 0 * (n + 1) + 0];
95 (*x_vertices)[5] = temp_x_vertices[ n * (n + 1) + 0];
96 (*y_vertices)[5] = temp_y_vertices[ n * (n + 1) + 0];
97 (*z_vertices)[5] = -temp_z_vertices[ n * (n + 1) + 0];
98 (*x_vertices)[6] = temp_x_vertices[ 0 * (n + 1) + n];
99 (*y_vertices)[6] = temp_y_vertices[ 0 * (n + 1) + n];
100 (*z_vertices)[6] = -temp_z_vertices[ 0 * (n + 1) + n];
101 (*x_vertices)[7] = temp_x_vertices[ n * (n + 1) + n];
102 (*y_vertices)[7] = temp_y_vertices[ n * (n + 1) + n];
103 (*z_vertices)[7] = -temp_z_vertices[ n * (n + 1) + n];
104
105 // Store the coordinates for the edges
106 for (unsigned i = 0; i < num_edge_coords; ++i) {
107 cube_edge_x_vertices[ 0 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + 0];
108 cube_edge_y_vertices[ 0 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + 0];
109 cube_edge_z_vertices[ 0 * num_edge_coords + i] = temp_z_vertices[(1 + i) * (n + 1) + 0];
110 cube_edge_x_vertices[ 1 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
111 cube_edge_y_vertices[ 1 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
112 cube_edge_z_vertices[ 1 * num_edge_coords + i] = temp_z_vertices[0 * (n + 1) + (1 + i)];
113 cube_edge_x_vertices[ 2 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
114 cube_edge_y_vertices[ 2 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
115 cube_edge_z_vertices[ 2 * num_edge_coords + i] = temp_z_vertices[n * (n + 1) + (1 + i)];
116 cube_edge_x_vertices[ 3 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + n];
117 cube_edge_y_vertices[ 3 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + n];
118 cube_edge_z_vertices[ 3 * num_edge_coords + i] = temp_z_vertices[(1 + i) * (n + 1) + n];
119 cube_edge_x_vertices[ 4 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + 0];
120 cube_edge_y_vertices[ 4 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + 0];
121 cube_edge_z_vertices[ 4 * num_edge_coords + i] = -temp_z_vertices[(1 + i) * (n + 1) + 0];
122 cube_edge_x_vertices[ 5 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
123 cube_edge_y_vertices[ 5 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
124 cube_edge_z_vertices[ 5 * num_edge_coords + i] = -temp_z_vertices[0 * (n + 1) + (1 + i)];
125 cube_edge_x_vertices[ 6 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
126 cube_edge_y_vertices[ 6 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
127 cube_edge_z_vertices[ 6 * num_edge_coords + i] = -temp_z_vertices[n * (n + 1) + (1 + i)];
128 cube_edge_x_vertices[ 7 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + n];
129 cube_edge_y_vertices[ 7 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + n];
130 cube_edge_z_vertices[ 7 * num_edge_coords + i] = -temp_z_vertices[(1 + i) * (n + 1) + n];
131 cube_edge_x_vertices[ 8 * num_edge_coords + i] = temp_z_vertices[0 * (n + 1) + (1 + i)];
132 cube_edge_y_vertices[ 8 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
133 cube_edge_z_vertices[ 8 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
134 cube_edge_x_vertices[ 9 * num_edge_coords + i] = temp_z_vertices[n * (n + 1) + (1 + i)];
135 cube_edge_y_vertices[ 9 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
136 cube_edge_z_vertices[ 9 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
137 cube_edge_x_vertices[10 * num_edge_coords + i] = -temp_z_vertices[0 * (n + 1) + (1 + i)];
138 cube_edge_y_vertices[10 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
139 cube_edge_z_vertices[10 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
140 cube_edge_x_vertices[11 * num_edge_coords + i] = -temp_z_vertices[n * (n + 1) + (1 + i)];
141 cube_edge_y_vertices[11 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
142 cube_edge_z_vertices[11 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
143 }
144
145 // Move the 12 edges to the final vertices array
146 unsigned Estart = 9 - 1;
147 unsigned Eend = Estart + 12 * num_edge_coords;
148 for (unsigned i = 0; i < 12 * num_edge_coords; ++i) {
149 (*x_vertices)[i + Estart] = cube_edge_x_vertices[i];
150 (*y_vertices)[i + Estart] = cube_edge_y_vertices[i];
151 (*z_vertices)[i + Estart] = cube_edge_z_vertices[i];
152 }
153
154 free(cube_edge_x_vertices);
155 free(cube_edge_y_vertices);
156 free(cube_edge_z_vertices);
157
158 // store the internal vertices for face 1
159 for (unsigned i = 0; i < num_edge_coords; ++i) {
160 for (unsigned j = 0; j < num_edge_coords; ++j) {
161
162 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
163 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
164 cube_inner_z_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
165 }
166 }
167
168 // Move face 1 to final Vertices array
169 unsigned Fstart = Eend;
170 for (unsigned i = 0; i < num_inner_coords; ++i) {
171 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
172 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
173 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
174 }
175
176 // store the internal vertices for face 2
177 for (unsigned i = 0; i < num_edge_coords; ++i) {
178 for (unsigned j = 0; j < num_edge_coords; ++j) {
179
180 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
181 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
182 cube_inner_z_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
183 }
184 }
185
186 // Move face 2 to final Vertices array
187 Fstart += num_inner_coords;
188 for (unsigned i = 0; i < num_inner_coords; ++i) {
189 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
190 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
191 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
192 }
193
194 // store the internal vertices for face 3
195 for (unsigned i = 0; i < num_edge_coords; ++i) {
196 for (unsigned j = 0; j < num_edge_coords; ++j) {
197
198 cube_inner_x_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
199 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
200 cube_inner_z_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
201 }
202 }
203
204 // Move face 3 to final Vertices array
205 Fstart += num_inner_coords;
206 for (unsigned i = 0; i < num_inner_coords; ++i) {
207 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
208 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
209 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
210 }
211
212 // store the internal vertices for face 4
213 for (unsigned i = 0; i < num_edge_coords; ++i) {
214 for (unsigned j = 0; j < num_edge_coords; ++j) {
215
216 cube_inner_x_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
217 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
218 cube_inner_z_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
219 }
220 }
221
222 // Move face 4 to final Vertices array
223 Fstart += num_inner_coords;
224 for (unsigned i = 0; i < num_inner_coords; ++i) {
225 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
226 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
227 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
228 }
229
230 // store the internal vertices for face 5
231 for (unsigned i = 0; i < num_edge_coords; ++i) {
232 for (unsigned j = 0; j < num_edge_coords; ++j) {
233
234 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
235 cube_inner_y_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
236 cube_inner_z_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
237 }
238 }
239
240 // Move face 5 to final Vertices array
241 Fstart += num_inner_coords;
242 for (unsigned i = 0; i < num_inner_coords; ++i) {
243 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
244 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
245 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
246 }
247
248 // store the internal vertices for face 6
249 for (unsigned i = 0; i < num_edge_coords; ++i) {
250 for (unsigned j = 0; j < num_edge_coords; ++j) {
251
252 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
253 cube_inner_y_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
254 cube_inner_z_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
255 }
256 }
257
258 // Move face 6 to final Vertices array
259 Fstart += num_inner_coords;
260 for (unsigned i = 0; i < num_inner_coords; ++i) {
261 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
262 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
263 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
264 }
265
266 free(temp_x_vertices);
267 free(temp_y_vertices);
268 free(temp_z_vertices);
269 free(cube_inner_x_vertices);
270 free(cube_inner_y_vertices);
271 free(cube_inner_z_vertices);
272
273 for (unsigned i = 0; i < 6; ++i)
274 for (unsigned j = 0; j < n * n; ++j)
275 (*face_id)[i * n * n + j] = i + 1;
276
277 unsigned * edge_vertices = xmalloc(num_edge_coords * 12 * sizeof(*edge_vertices));
278
279 Fstart = Estart + 12 * num_edge_coords;
280
281 for (unsigned i = 0; i < num_edge_coords * 12; ++i) edge_vertices[i] = i + Estart;
282
283 unsigned * F = xmalloc((n + 1) * (n + 1) * sizeof(*F));
284 unsigned cell_to_vertex_offset = 0;
285
286// Calculate the indices for all the vertices on face 1
287// v2-----e3-----v4
288// | |
289// e1 f1 e4
290// | |
291// v1-----e2-----v3
292
293 F[0] = 0; F[n] = 2; F[n * (n + 1) + 0] = 1; F[n * (n + 1) + n] = 3;
294 for (unsigned i = 0; i < num_edge_coords; ++i) {
295 F[i + 1] = edge_vertices[1 * num_edge_coords + i];
296 F[(i + 1) * (n + 1)] = edge_vertices[0 * num_edge_coords + i];
297 F[(i + 1) * (n + 1) + n] = edge_vertices[3 * num_edge_coords + i];
298 F[n * (n + 1) + i + 1] = edge_vertices[2 * num_edge_coords + i];
299 }
300
301 for (unsigned i = 0; i < n-1; ++i)
302 for (unsigned j = 0; j < n-1; ++j)
303 F[(i + 1) * (n + 1) + (j + 1)] = Fstart + i * (n - 1) + j;
304
305 for (unsigned i = 0; i < n; ++i) {
306 for (unsigned j = 0; j < n; ++j) {
307 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
308 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
309 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
310 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
311 }
312 }
313
314// Calculate the indices for all the vertices on face 2
315// v6-----e7-----v8
316// | |
317// e5 f2 e8
318// | |
319// v5-----e6-----v7
320
321 F[0] = 4; F[n] = 6; F[n * (n + 1) + 0] = 5; F[n * (n + 1) + n] = 7;
322 for (unsigned i = 0; i < num_edge_coords; ++i) {
323 F[i + 1] = edge_vertices[5 * num_edge_coords + i];
324 F[(i + 1) * (n + 1)] = edge_vertices[4 * num_edge_coords + i];
325 F[(i + 1) * (n + 1) + n] = edge_vertices[7 * num_edge_coords + i];
326 F[n * (n + 1) + i + 1] = edge_vertices[6 * num_edge_coords + i];
327 }
328
329 for (unsigned i = 0; i < n-1; ++i)
330 for (unsigned j = 0; j < n-1; ++j)
331 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
332
333 for (unsigned i = 0; i < n; ++i) {
334 for (unsigned j = 0; j < n; ++j) {
335 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
336 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
337 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
338 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
339 }
340 }
341
342// Calculate the indices for all the vertices on face 3
343// v2-----e10----v6
344// | |
345// e1 f3 e5
346// | |
347// v1-----e9-----v5
348
349 F[0] = 0; F[n] = 4; F[n * (n + 1) + 0] = 1; F[n * (n + 1) + n] = 5;
350 for (unsigned i = 0; i < num_edge_coords; ++i) {
351 F[i + 1] = edge_vertices[8 * num_edge_coords + i];
352 F[(i + 1) * (n + 1)] = edge_vertices[0 * num_edge_coords + i];
353 F[(i + 1) * (n + 1) + n] = edge_vertices[4 * num_edge_coords + i];
354 F[n * (n + 1) + i + 1] = edge_vertices[9 * num_edge_coords + i];
355 }
356
357 for (unsigned i = 0; i < n-1; ++i)
358 for (unsigned j = 0; j < n-1; ++j)
359 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
360
361 for (unsigned i = 0; i < n; ++i) {
362 for (unsigned j = 0; j < n; ++j) {
363 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
364 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
365 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
366 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
367 }
368 }
369
370// Calculate the indices for all the vertices on face 4
371// v4-----e12----v8
372// | |
373// e4 f4 e8
374// | |
375// v3-----e11----v7
376
377 F[0] = 2; F[n] = 6; F[n * (n + 1) + 0] = 3; F[n * (n + 1) + n] = 7;
378 for (unsigned i = 0; i < num_edge_coords; ++i) {
379 F[i + 1] = edge_vertices[10 * num_edge_coords + i];
380 F[(i + 1) * (n + 1)] = edge_vertices[ 3 * num_edge_coords + i];
381 F[(i + 1) * (n + 1) + n] = edge_vertices[ 7 * num_edge_coords + i];
382 F[n * (n + 1) + i + 1] = edge_vertices[11 * num_edge_coords + i];
383 }
384
385 for (unsigned i = 0; i < n-1; ++i)
386 for (unsigned j = 0; j < n-1; ++j)
387 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
388
389 for (unsigned i = 0; i < n; ++i) {
390 for (unsigned j = 0; j < n; ++j) {
391 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
392 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
393 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
394 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
395 }
396 }
397
398// Calculate the indices for all the vertices on face 5
399// v5-----e6-----v7
400// | |
401// e9 f5 e11
402// | |
403// v1-----e2-----v3
404
405 F[0] = 0; F[n] = 2; F[n * (n + 1) + 0] = 4; F[n * (n + 1) + n] = 6;
406 for (unsigned i = 0; i < num_edge_coords; ++i) {
407 F[i + 1] = edge_vertices[ 1 * num_edge_coords + i];
408 F[(i + 1) * (n + 1)] = edge_vertices[ 8 * num_edge_coords + i];
409 F[(i + 1) * (n + 1) + n] = edge_vertices[10 * num_edge_coords + i];
410 F[n * (n + 1) + i + 1] = edge_vertices[ 5 * num_edge_coords + i];
411 }
412
413 for (unsigned i = 0; i < n-1; ++i)
414 for (unsigned j = 0; j < n-1; ++j)
415 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
416
417 for (unsigned i = 0; i < n; ++i) {
418 for (unsigned j = 0; j < n; ++j) {
419 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
420 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
421 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
422 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
423 }
424 }
425
426// Calculate the indices for all the vertices on face 6
427// v6-----e7-----v8
428// | |
429// e10 f6 e12
430// | |
431// v2-----e3-----v4
432
433 F[0] = 1; F[n] = 3; F[n * (n + 1) + 0] = 5; F[n * (n + 1) + n] = 7;
434 for (unsigned i = 0; i < num_edge_coords; ++i) {
435 F[i + 1] = edge_vertices[ 2 * num_edge_coords + i];
436 F[(i + 1) * (n + 1)] = edge_vertices[ 9 * num_edge_coords + i];
437 F[(i + 1) * (n + 1) + n] = edge_vertices[11 * num_edge_coords + i];
438 F[n * (n + 1) + i + 1] = edge_vertices[ 6 * num_edge_coords + i];
439 }
440
441 for (unsigned i = 0; i < n-1; ++i)
442 for (unsigned j = 0; j < n-1; ++j)
443 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
444
445 for (unsigned i = 0; i < n; ++i) {
446 for (unsigned j = 0; j < n; ++j) {
447 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
448 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
449 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
450 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
451 }
452 }
453
454 free(edge_vertices);
455 free(F);
456}
457
459
460 unsigned num_cells, num_vertices;
461 double * x_vertices, * y_vertices, * z_vertices;
462 unsigned * cell_to_vertex;
463 unsigned * face_id;
464
466 n, &num_cells, &num_vertices, &x_vertices, &y_vertices, &z_vertices,
467 &cell_to_vertex, &face_id);
468
469 free(face_id);
470
473 for (unsigned i = 0; i < num_cells; ++i) num_vertices_per_cell[i] = 4;
474
478 x_vertices, y_vertices, (int *)cell_to_vertex);
479 grid_data.cell_ids =
480 xmalloc((size_t)num_cells * sizeof(*(grid_data.cell_ids)));
481 grid_data.vertex_ids =
482 xmalloc((size_t)num_vertices * sizeof(*(grid_data.vertex_ids)));
483 grid_data.edge_ids =
484 xmalloc(grid_data.num_edges * sizeof(*(grid_data.edge_ids)));
485 for (size_t i = 0; i < (size_t)num_cells; ++i)
486 grid_data.cell_ids[i] = (yac_int)i;
487 for (size_t i = 0; i < (size_t)num_vertices; ++i) {
488 grid_data.vertex_ids[i] = (yac_int)i;
489 grid_data.vertex_coordinates[i][0] = x_vertices[i];
490 grid_data.vertex_coordinates[i][1] = y_vertices[i];
491 grid_data.vertex_coordinates[i][2] = z_vertices[i];
492 }
493 for (size_t i = 0; i < grid_data.num_edges; ++i)
494 grid_data.edge_ids[i] = (yac_int)i;
496 free(x_vertices);
497 free(y_vertices);
498 free(z_vertices);
499 free(cell_to_vertex);
500
501 return grid_data;
502}
503
505 char const * name, size_t n) {
506
507 return
509 name,
510 yac_generate_cubed_sphere_grid((unsigned)n));
511}
512
513static void decompose_domain_simple(unsigned n, int size, int * cell_owner) {
514
515 unsigned nbr_cells = n * n * 6;
516
517 for (unsigned i = 0; i < nbr_cells; ++i)
518 cell_owner[i] = (i * size + size - 1) / nbr_cells;
519}
520
521static void decompose_domain_2d(unsigned n, int size, int * cell_owner) {
522
523 // distribute processes among all six sides of the cubed sphere
524 int proc_start[6 + 1];
525 for (int i = 0; i < 6 + 1; ++i)
526 proc_start[i] = (size * i) / 6;
527
528 // for each side of the cube
529 for (unsigned i = 0, offset = 0; i < 6; ++i) {
530
531 int num_procs[2];
532 int curr_num_procs = proc_start[i+1] - proc_start[i];
533
534 yac_generate_reg2d_decomp((int[2]){(int)n,(int)n}, curr_num_procs, num_procs);
535
536 int local_start_x[num_procs[0] + 1];
537 int local_start_y[num_procs[1] + 1];
538
539 for (int j = 0; j < num_procs[0] + 1; ++j)
540 local_start_x[j] = (n * j) / num_procs[0];
541 for (int j = 0; j < num_procs[1] + 1; ++j)
542 local_start_y[j] = (n * j) / num_procs[1];
543
544 for (int p_y = 0; p_y < num_procs[1]; ++p_y) {
545 for (int row = local_start_y[p_y]; row < local_start_y[p_y+1]; ++row) {
546 for (int p_x = 0; p_x < num_procs[0]; ++p_x) {
547 for (int column = local_start_x[p_x]; column < local_start_x[p_x+1];
548 ++column) {
549
550 cell_owner[offset++] = proc_start[i] + p_x + num_procs[0] * p_y;
551 }
552 }
553 }
554 }
555 }
556}
557
558static void decompose_domain(unsigned n, int size, int * cell_owner) {
559
560 if (size <= 6) {
561 decompose_domain_simple(n, size, cell_owner);
562 } else {
563 decompose_domain_2d(n, size, cell_owner);
564 }
565}
566
568 unsigned n, unsigned * nbr_vertices, unsigned * nbr_cells,
569 unsigned ** num_vertices_per_cell, unsigned ** cell_to_vertex,
570 double ** x_vertices, double ** y_vertices, double ** x_cells,
571 double ** y_cells, int ** global_cell_id, int ** cell_core_mask,
572 int ** global_corner_id, int ** corner_core_mask, int rank, int size) {
573
574 double * x_vertices_3d, * y_vertices_3d, * z_vertices_3d;
575 unsigned * dummy;
576
577 // generate global grid
578 yac_generate_cubed_sphere_grid_information(n, nbr_cells, nbr_vertices,
579 &x_vertices_3d, &y_vertices_3d,
580 &z_vertices_3d, cell_to_vertex,
581 &dummy);
582 *num_vertices_per_cell = xmalloc(*nbr_cells * sizeof(**num_vertices_per_cell));
583 for (unsigned i = 0; i < *nbr_cells; ++i) (*num_vertices_per_cell)[i] = 4;
584 free(dummy);
585
586 double * x_cell_3d = xmalloc(*nbr_cells * sizeof(*x_cell_3d));
587 double * y_cell_3d = xmalloc(*nbr_cells * sizeof(*y_cell_3d));
588 double * z_cell_3d = xmalloc(*nbr_cells * sizeof(*z_cell_3d));
589
590 for (unsigned i = 0; i < *nbr_cells; ++i) {
591
592 x_cell_3d[i] = y_cell_3d[i] = z_cell_3d[i] = 0;
593
594 for (unsigned j = 0; j < 4; ++j) {
595
596 x_cell_3d[i] += x_vertices_3d[(*cell_to_vertex)[4 * i + j]];
597 y_cell_3d[i] += y_vertices_3d[(*cell_to_vertex)[4 * i + j]];
598 z_cell_3d[i] += z_vertices_3d[(*cell_to_vertex)[4 * i + j]];
599 }
600
601 double scale = 1.0 / sqrt(x_cell_3d[i] * x_cell_3d[i] +
602 y_cell_3d[i] * y_cell_3d[i] +
603 z_cell_3d[i] * z_cell_3d[i]);
604
605 x_cell_3d[i] *= scale;
606 y_cell_3d[i] *= scale;
607 z_cell_3d[i] *= scale;
608 }
609
610 int * cell_is_on_rank = xmalloc(*nbr_cells * sizeof(*cell_is_on_rank));
611
612 decompose_domain(n, size, cell_is_on_rank);
613
614 // mask for required vertices and cells
615 int * required_vertices = xcalloc(*nbr_vertices, sizeof(*required_vertices));
616 int * required_cells = xcalloc(*nbr_cells, sizeof(*required_cells));
617
618 // mark all local cells and their vertices as required
619 for (unsigned i = 0, offset = 0; i < *nbr_cells;
620 offset += (*num_vertices_per_cell)[i++]) {
621
622 if (cell_is_on_rank[i] != rank) continue;
623
624 cell_is_on_rank[i] = -1;
625
626 required_cells[i] = 1;
627 for (unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j)
628 required_vertices[(*cell_to_vertex)[offset+j]] = 1;
629 }
630
631 // mark all halo cells as required and generate cell_is_on_rank
632 for (unsigned i = 0, offset = 0; i < *nbr_cells;
633 offset += (*num_vertices_per_cell)[i++]) {
634
635 if (!required_cells[i]) {
636
637 for (unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
638
639 if (required_vertices[(*cell_to_vertex)[offset+j]]) {
640
641 required_cells[i] = 1;
642 break;
643 }
644 }
645
646 }
647 }
648
649 // mask for halo vertices
650 int * vertex_is_on_rank = xmalloc(*nbr_vertices * sizeof(*vertex_is_on_rank));
651
652 // mark all halo vertices as required and generate vertex_is_on_rank
653 for (unsigned i = 0; i < *nbr_vertices; ++i)
654 if (required_vertices[i])
655 vertex_is_on_rank[i] = -1;
656 for (unsigned i = 0, offset = 0; i < *nbr_cells;
657 offset += (*num_vertices_per_cell)[i++]) {
658
659 if (required_cells[i] && cell_is_on_rank[i] != -1) {
660
661 for (unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
662
663 if (!required_vertices[(*cell_to_vertex)[offset+j]]) {
664
665 required_vertices[(*cell_to_vertex)[offset+j]] = 1;
666 vertex_is_on_rank[(*cell_to_vertex)[offset+j]] = cell_is_on_rank[i];
667 }
668 }
669 }
670 }
671
672 // count the number of cells and vertices
673 int part_num_vertices = 0;
674 int part_num_cells = 0;
675 for (unsigned i = 0; i < *nbr_vertices; ++i)
676 if (required_vertices[i])
677 part_num_vertices++;
678 for (unsigned i = 0; i < *nbr_cells; ++i)
679 if(required_cells[i])
680 part_num_cells++;
681
682 *global_cell_id = xmalloc(part_num_cells * sizeof(**global_cell_id));
683 *cell_core_mask = xmalloc(part_num_cells * sizeof(**cell_core_mask));
684 *global_corner_id = xmalloc(part_num_vertices * sizeof(**global_corner_id));
685 *corner_core_mask = xmalloc(part_num_vertices * sizeof(**corner_core_mask));
686
687 *x_cells = xmalloc(part_num_cells * sizeof(**x_cells));
688 *y_cells = xmalloc(part_num_cells * sizeof(**y_cells));
689 *x_vertices = xmalloc(part_num_vertices * sizeof(**x_vertices));
690 *y_vertices = xmalloc(part_num_vertices * sizeof(**y_vertices));
691
692 // generate final vertex data
693 part_num_vertices = 0;
694 int * global_to_local_vertex = xmalloc(*nbr_vertices * sizeof(*global_to_local_vertex));
695 for (unsigned i = 0; i < *nbr_vertices; ++i) {
696
697 if (required_vertices[i]) {
698
699 (*global_corner_id)[part_num_vertices] = i;
700 (*corner_core_mask)[part_num_vertices] = vertex_is_on_rank[i] == -1;
701 double p[3] = {x_vertices_3d[i], y_vertices_3d[i], z_vertices_3d[i]};
702 XYZtoLL(p, (*x_vertices) + part_num_vertices, (*y_vertices) + part_num_vertices);
703 global_to_local_vertex[i] = part_num_vertices;
704 part_num_vertices++;
705 }
706 }
707
708 free(vertex_is_on_rank);
709 *nbr_vertices = part_num_vertices;
710 free(required_vertices);
711
712 // generate final cell data
713 int num_cell_vertex_dependencies = 0;
714 part_num_cells = 0;
715 for (unsigned i = 0, offset = 0; i < *nbr_cells;
716 offset += (*num_vertices_per_cell)[i++]) {
717
718 if (required_cells[i]) {
719
720 (*global_cell_id)[part_num_cells] = i;
721 (*cell_core_mask)[part_num_cells] = (cell_is_on_rank[i] == -1);
722 double middle_point[3] = {x_cell_3d[i],y_cell_3d[i],z_cell_3d[i]};
723
724 for (unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
725 (*cell_to_vertex)[num_cell_vertex_dependencies+j] =
726 global_to_local_vertex[(*cell_to_vertex)[offset+j]];
727 }
728
729 XYZtoLL(middle_point, (*x_cells)+part_num_cells, (*y_cells)+part_num_cells);
730
731 num_cell_vertex_dependencies += (*num_vertices_per_cell)[i];
732
733 (*num_vertices_per_cell)[part_num_cells] = (*num_vertices_per_cell)[i];
734
735 part_num_cells++;
736 }
737 }
738
739 free(x_cell_3d);
740 free(y_cell_3d);
741 free(z_cell_3d);
742 free(x_vertices_3d);
743 free(y_vertices_3d);
744 free(z_vertices_3d);
745
746 *num_vertices_per_cell = xrealloc(*num_vertices_per_cell, part_num_cells *
747 sizeof(**num_vertices_per_cell));
748 *cell_to_vertex = xrealloc(*cell_to_vertex, num_cell_vertex_dependencies *
749 sizeof(**cell_to_vertex));
750 free(cell_is_on_rank);
751 *nbr_cells = part_num_cells;
752 free(required_cells);
753 free(global_to_local_vertex);
754}
755
756#ifdef YAC_NETCDF_ENABLED
757void yac_write_cubed_sphere_grid(unsigned n, char const * filename) {
758
759 // generate basic grid information
760 unsigned num_cells, num_vertices;
761 double * x_vertices, * y_vertices, * z_vertices;
762 unsigned * vertices_of_cell, * dummy;
764 n, &num_cells, &num_vertices, &x_vertices, &y_vertices, &z_vertices,
765 &vertices_of_cell, &dummy);
766 free(dummy);
767
768 // create file and define dimensions and variables
769 size_t nv = 4, ne = 4;
770 int ncid;
771 int dim_cell_id, dim_vertex_id, dim_nv_id, dim_ne_id;
772 int var_vlon_id, var_vlat_id, var_clon_id, var_clat_id, var_mask_id,
773 var_v2c_id, var_c2v_id;
774 yac_nc_create(filename, NC_CLOBBER, &ncid);
775 YAC_HANDLE_ERROR(nc_def_dim(ncid, "cell", (size_t)num_cells, &dim_cell_id));
776 YAC_HANDLE_ERROR(nc_def_dim(ncid, "vertex", (size_t)num_vertices, &dim_vertex_id));
777 YAC_HANDLE_ERROR(nc_def_dim(ncid, "nv", 4, &dim_nv_id));
778 YAC_HANDLE_ERROR(nc_def_dim(ncid, "ne", 4, &dim_ne_id));
780 nc_def_var(ncid, "vlon", NC_DOUBLE, 1, &dim_vertex_id, &var_vlon_id));
782 nc_def_var(ncid, "vlat", NC_DOUBLE, 1, &dim_vertex_id, &var_vlat_id));
784 nc_def_var(ncid, "clon", NC_DOUBLE, 1, &dim_cell_id, &var_clon_id));
786 nc_def_var(ncid, "clat", NC_DOUBLE, 1, &dim_cell_id, &var_clat_id));
788 nc_def_var(
789 ncid, "cell_sea_land_mask", NC_INT, 1, &dim_cell_id, &var_mask_id));
791 nc_def_var(
792 ncid, "vertex_of_cell", NC_INT, 2, (int[]){dim_nv_id, dim_cell_id},
793 &var_c2v_id));
795 nc_def_var(
796 ncid, "cells_of_vertex", NC_INT, 2, (int[]){dim_ne_id, dim_vertex_id},
797 &var_v2c_id));
798 YAC_HANDLE_ERROR(nc_enddef(ncid));
799
800 // generate and write grid data
801
802 double * lon_buffer =
803 xmalloc(MAX(num_cells, num_vertices) * sizeof(*lon_buffer));
804 double * lat_buffer =
805 xmalloc(MAX(num_cells, num_vertices) * sizeof(*lat_buffer));
806
807 for (unsigned i = 0; i < num_vertices; ++i) {
808 double coord_3d[3] = {x_vertices[i], y_vertices[i], z_vertices[i]};
809 XYZtoLL(coord_3d, lon_buffer + i, lat_buffer + i);
810 }
811 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_vlon_id, lon_buffer));
812 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_vlat_id, lat_buffer));
813
814 for (unsigned i = 0; i < num_cells; ++i) {
815 double coord_3d[3] = {0.0, 0.0, 0.0};
816 for (unsigned j = 0; j < 4; ++j) {
817 unsigned vertex_idx = vertices_of_cell[4 * i + j];
818 coord_3d[0] += x_vertices[vertex_idx];
819 coord_3d[1] += y_vertices[vertex_idx];
820 coord_3d[2] += z_vertices[vertex_idx];
821 }
822 normalise_vector(coord_3d);
823 XYZtoLL(coord_3d, lon_buffer + i, lat_buffer + i);
824 }
825 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_clon_id, lon_buffer));
826 YAC_HANDLE_ERROR(nc_put_var_double(ncid, var_clat_id, lat_buffer));
827
828 free(lat_buffer);
829 free(lon_buffer);
830 free(z_vertices);
831 free(y_vertices);
832 free(x_vertices);
833
834 int * int_buffer =
835 xmalloc(MAX(nv, ne) * MAX(num_cells, num_vertices) * sizeof(*int_buffer));
836
837 for (unsigned i = 0; i < num_cells; ++i) int_buffer[i] = 1;
838 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_mask_id, int_buffer));
839
840 for (unsigned i = 0; i < ne * num_vertices; ++i)
841 int_buffer[i] = INT_MAX;
842 for (unsigned i = 0; i < num_cells; ++i) {
843 int cell_id = (int)(i + 1);
844 for (unsigned j = 0; j < ne; ++j) {
845 unsigned vertex_idx = vertices_of_cell[ne * i + j];
846 for (unsigned k = 0; k < ne; ++k) {
847 if (int_buffer[k * num_vertices + vertex_idx] == INT_MAX) {
848 int_buffer[k * num_vertices + vertex_idx] = cell_id;
849 break;
850 }
851 }
852 }
853 }
854 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_v2c_id, int_buffer));
855
856 for (unsigned i = 0; i < num_cells; ++i)
857 for (unsigned j = 0; j < nv; ++j)
858 int_buffer[j * num_cells + i] =
859 vertices_of_cell[nv * i + j] + 1;
860 YAC_HANDLE_ERROR(nc_put_var_int(ncid, var_c2v_id, int_buffer));
861
862 free(int_buffer);
863 free(vertices_of_cell);
864
865 YAC_HANDLE_ERROR(nc_close(ncid));
866}
867#else
868void yac_write_cubed_sphere_grid(unsigned n, char const * filename) {
869
870 UNUSED(n);
871 UNUSED(filename);
872 die(
873 "ERROR(yac_write_cubed_sphere_grid): "
874 "YAC is built without the NetCDF support");
875}
876#endif // YAC_NETCDF_ENABLED
#define YAC_ASSERT(exp, msg)
struct yac_basic_grid * yac_basic_grid_new(char const *name, struct yac_basic_grid_data grid_data)
Definition basic_grid.c:50
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)
#define UNUSED(x)
Definition core.h:73
void yac_generate_part_cube_grid_information(unsigned n, unsigned *nbr_vertices, unsigned *nbr_cells, unsigned **num_vertices_per_cell, unsigned **cell_to_vertex, double **x_vertices, double **y_vertices, double **x_cells, double **y_cells, int **global_cell_id, int **cell_core_mask, int **global_corner_id, int **corner_core_mask, int rank, int size)
void yac_write_cubed_sphere_grid(unsigned n, char const *filename)
static void decompose_domain(unsigned n, int size, int *cell_owner)
struct yac_basic_grid * yac_generate_cubed_sphere_basic_grid(char const *name, size_t n)
static void decompose_domain_2d(unsigned n, int size, int *cell_owner)
static void decompose_domain_simple(unsigned n, int size, int *cell_owner)
struct yac_basic_grid_data yac_generate_cubed_sphere_grid(unsigned n)
void yac_generate_cubed_sphere_grid_information(unsigned n, unsigned *num_cells, unsigned *num_vertices, double **x_vertices, double **y_vertices, double **z_vertices, unsigned **cell_to_vertex, unsigned **face_id)
void yac_generate_reg2d_decomp(int num_points[2], int total_num_procs, int *num_procs)
static void normalise_vector(double v[])
Definition geometry.h:650
void yac_nc_create(const char *path, int cmode, int *ncidp)
Definition io_utils.c:367
static void XYZtoLL(double const p_in[], double *lon, double *lat)
add versions of standard API functions not returning on error
#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
int * cell_to_vertex
size_t num_cells[2]
int * cell_core_mask
#define YAC_HANDLE_ERROR(exp)
Definition toy_output.c:13
char const * name
Definition toy_scrip.c:114
double(* p)(double lon, double lat)
Definition toy_scrip.c:119
#define MAX(a, b)
#define die(msg)
Definition yac_assert.h:12
Xt_int yac_int
Definition yac_types.h:15