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