21 unsigned n,
unsigned *
num_cells,
unsigned * num_vertices,
22 double ** x_vertices,
double ** y_vertices,
double ** z_vertices,
25 YAC_ASSERT((n >= 1) && (n <= 40000),
"invalid number of linear subdivisions")
31 *x_vertices =
xmalloc(*num_vertices *
sizeof(**x_vertices));
32 *y_vertices =
xmalloc(*num_vertices *
sizeof(**y_vertices));
33 *z_vertices =
xmalloc(*num_vertices *
sizeof(**z_vertices));
40 unsigned num_edge_coords = n - 1;
41 double * cube_edge_x_vertices =
xmalloc(12 * num_edge_coords *
sizeof(*cube_edge_x_vertices));
42 double * cube_edge_y_vertices =
xmalloc(12 * num_edge_coords *
sizeof(*cube_edge_y_vertices));
43 double * cube_edge_z_vertices =
xmalloc(12 * num_edge_coords *
sizeof(*cube_edge_z_vertices));
45 unsigned num_inner_coords = num_edge_coords * num_edge_coords;
46 double * cube_inner_x_vertices =
xmalloc(num_inner_coords *
sizeof(*cube_inner_x_vertices));
47 double * cube_inner_y_vertices =
xmalloc(num_inner_coords *
sizeof(*cube_inner_y_vertices));
48 double * cube_inner_z_vertices =
xmalloc(num_inner_coords *
sizeof(*cube_inner_z_vertices));
50 double * temp_x_vertices =
xmalloc((n + 1) * (n + 1) *
sizeof(*temp_x_vertices));
51 double * temp_y_vertices =
xmalloc((n + 1) * (n + 1) *
sizeof(*temp_y_vertices));
52 double * temp_z_vertices =
xmalloc((n + 1) * (n + 1) *
sizeof(*temp_z_vertices));
55 double * theta =
xmalloc((n + 1) *
sizeof(*theta));
58 double d = (M_PI_2 / (double)n);
59 for (
unsigned i = 0; i < n + 1; ++i)
60 theta[i] = tan(- M_PI_4 + d * (
double)i);
63 for (
unsigned i = 0; i < n + 1; ++i) {
64 for (
unsigned j = 0; j < n + 1; ++j) {
67 sqrt(1.0 / (theta[i] * theta[i] + theta[j] * theta[j] + 1.0));
69 temp_x_vertices[i * (n + 1) + j] = theta[j] * scale;
70 temp_y_vertices[i * (n + 1) + j] = theta[i] * scale;
71 temp_z_vertices[i * (n + 1) + j] = - scale;
78 (*x_vertices)[0] = temp_x_vertices[ 0 * (n + 1) + 0];
79 (*y_vertices)[0] = temp_y_vertices[ 0 * (n + 1) + 0];
80 (*z_vertices)[0] = temp_z_vertices[ 0 * (n + 1) + 0];
81 (*x_vertices)[1] = temp_x_vertices[ n * (n + 1) + 0];
82 (*y_vertices)[1] = temp_y_vertices[ n * (n + 1) + 0];
83 (*z_vertices)[1] = temp_z_vertices[ n * (n + 1) + 0];
84 (*x_vertices)[2] = temp_x_vertices[ 0 * (n + 1) + n];
85 (*y_vertices)[2] = temp_y_vertices[ 0 * (n + 1) + n];
86 (*z_vertices)[2] = temp_z_vertices[ 0 * (n + 1) + n];
87 (*x_vertices)[3] = temp_x_vertices[ n * (n + 1) + n];
88 (*y_vertices)[3] = temp_y_vertices[ n * (n + 1) + n];
89 (*z_vertices)[3] = temp_z_vertices[ n * (n + 1) + n];
91 (*x_vertices)[4] = temp_x_vertices[ 0 * (n + 1) + 0];
92 (*y_vertices)[4] = temp_y_vertices[ 0 * (n + 1) + 0];
93 (*z_vertices)[4] = -temp_z_vertices[ 0 * (n + 1) + 0];
94 (*x_vertices)[5] = temp_x_vertices[ n * (n + 1) + 0];
95 (*y_vertices)[5] = temp_y_vertices[ n * (n + 1) + 0];
96 (*z_vertices)[5] = -temp_z_vertices[ n * (n + 1) + 0];
97 (*x_vertices)[6] = temp_x_vertices[ 0 * (n + 1) + n];
98 (*y_vertices)[6] = temp_y_vertices[ 0 * (n + 1) + n];
99 (*z_vertices)[6] = -temp_z_vertices[ 0 * (n + 1) + n];
100 (*x_vertices)[7] = temp_x_vertices[ n * (n + 1) + n];
101 (*y_vertices)[7] = temp_y_vertices[ n * (n + 1) + n];
102 (*z_vertices)[7] = -temp_z_vertices[ n * (n + 1) + n];
105 for (
unsigned i = 0; i < num_edge_coords; ++i) {
106 cube_edge_x_vertices[ 0 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + 0];
107 cube_edge_y_vertices[ 0 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + 0];
108 cube_edge_z_vertices[ 0 * num_edge_coords + i] = temp_z_vertices[(1 + i) * (n + 1) + 0];
109 cube_edge_x_vertices[ 1 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
110 cube_edge_y_vertices[ 1 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
111 cube_edge_z_vertices[ 1 * num_edge_coords + i] = temp_z_vertices[0 * (n + 1) + (1 + i)];
112 cube_edge_x_vertices[ 2 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
113 cube_edge_y_vertices[ 2 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
114 cube_edge_z_vertices[ 2 * num_edge_coords + i] = temp_z_vertices[n * (n + 1) + (1 + i)];
115 cube_edge_x_vertices[ 3 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + n];
116 cube_edge_y_vertices[ 3 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + n];
117 cube_edge_z_vertices[ 3 * num_edge_coords + i] = temp_z_vertices[(1 + i) * (n + 1) + n];
118 cube_edge_x_vertices[ 4 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + 0];
119 cube_edge_y_vertices[ 4 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + 0];
120 cube_edge_z_vertices[ 4 * num_edge_coords + i] = -temp_z_vertices[(1 + i) * (n + 1) + 0];
121 cube_edge_x_vertices[ 5 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
122 cube_edge_y_vertices[ 5 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
123 cube_edge_z_vertices[ 5 * num_edge_coords + i] = -temp_z_vertices[0 * (n + 1) + (1 + i)];
124 cube_edge_x_vertices[ 6 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
125 cube_edge_y_vertices[ 6 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
126 cube_edge_z_vertices[ 6 * num_edge_coords + i] = -temp_z_vertices[n * (n + 1) + (1 + i)];
127 cube_edge_x_vertices[ 7 * num_edge_coords + i] = temp_x_vertices[(1 + i) * (n + 1) + n];
128 cube_edge_y_vertices[ 7 * num_edge_coords + i] = temp_y_vertices[(1 + i) * (n + 1) + n];
129 cube_edge_z_vertices[ 7 * num_edge_coords + i] = -temp_z_vertices[(1 + i) * (n + 1) + n];
130 cube_edge_x_vertices[ 8 * num_edge_coords + i] = temp_z_vertices[0 * (n + 1) + (1 + i)];
131 cube_edge_y_vertices[ 8 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
132 cube_edge_z_vertices[ 8 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
133 cube_edge_x_vertices[ 9 * num_edge_coords + i] = temp_z_vertices[n * (n + 1) + (1 + i)];
134 cube_edge_y_vertices[ 9 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
135 cube_edge_z_vertices[ 9 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
136 cube_edge_x_vertices[10 * num_edge_coords + i] = -temp_z_vertices[0 * (n + 1) + (1 + i)];
137 cube_edge_y_vertices[10 * num_edge_coords + i] = temp_y_vertices[0 * (n + 1) + (1 + i)];
138 cube_edge_z_vertices[10 * num_edge_coords + i] = temp_x_vertices[0 * (n + 1) + (1 + i)];
139 cube_edge_x_vertices[11 * num_edge_coords + i] = -temp_z_vertices[n * (n + 1) + (1 + i)];
140 cube_edge_y_vertices[11 * num_edge_coords + i] = temp_y_vertices[n * (n + 1) + (1 + i)];
141 cube_edge_z_vertices[11 * num_edge_coords + i] = temp_x_vertices[n * (n + 1) + (1 + i)];
145 unsigned Estart = 9 - 1;
146 unsigned Eend = Estart + 12 * num_edge_coords;
147 for (
unsigned i = 0; i < 12 * num_edge_coords; ++i) {
148 (*x_vertices)[i + Estart] = cube_edge_x_vertices[i];
149 (*y_vertices)[i + Estart] = cube_edge_y_vertices[i];
150 (*z_vertices)[i + Estart] = cube_edge_z_vertices[i];
153 free(cube_edge_x_vertices);
154 free(cube_edge_y_vertices);
155 free(cube_edge_z_vertices);
158 for (
unsigned i = 0; i < num_edge_coords; ++i) {
159 for (
unsigned j = 0; j < num_edge_coords; ++j) {
161 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
162 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
163 cube_inner_z_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
168 unsigned Fstart = Eend;
169 for (
unsigned i = 0; i < num_inner_coords; ++i) {
170 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
171 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
172 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
176 for (
unsigned i = 0; i < num_edge_coords; ++i) {
177 for (
unsigned j = 0; j < num_edge_coords; ++j) {
179 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
180 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
181 cube_inner_z_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
186 Fstart += num_inner_coords;
187 for (
unsigned i = 0; i < num_inner_coords; ++i) {
188 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
189 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
190 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
194 for (
unsigned i = 0; i < num_edge_coords; ++i) {
195 for (
unsigned j = 0; j < num_edge_coords; ++j) {
197 cube_inner_x_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
198 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
199 cube_inner_z_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
204 Fstart += num_inner_coords;
205 for (
unsigned i = 0; i < num_inner_coords; ++i) {
206 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
207 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
208 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
212 for (
unsigned i = 0; i < num_edge_coords; ++i) {
213 for (
unsigned j = 0; j < num_edge_coords; ++j) {
215 cube_inner_x_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
216 cube_inner_y_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
217 cube_inner_z_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
222 Fstart += num_inner_coords;
223 for (
unsigned i = 0; i < num_inner_coords; ++i) {
224 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
225 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
226 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
230 for (
unsigned i = 0; i < num_edge_coords; ++i) {
231 for (
unsigned j = 0; j < num_edge_coords; ++j) {
233 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
234 cube_inner_y_vertices[i * num_edge_coords + j] = temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
235 cube_inner_z_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
240 Fstart += num_inner_coords;
241 for (
unsigned i = 0; i < num_inner_coords; ++i) {
242 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
243 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
244 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
248 for (
unsigned i = 0; i < num_edge_coords; ++i) {
249 for (
unsigned j = 0; j < num_edge_coords; ++j) {
251 cube_inner_x_vertices[i * num_edge_coords + j] = temp_x_vertices[(1 + i) * (n + 1) + (1 + j)];
252 cube_inner_y_vertices[i * num_edge_coords + j] = -temp_z_vertices[(1 + i) * (n + 1) + (1 + j)];
253 cube_inner_z_vertices[i * num_edge_coords + j] = temp_y_vertices[(1 + i) * (n + 1) + (1 + j)];
258 Fstart += num_inner_coords;
259 for (
unsigned i = 0; i < num_inner_coords; ++i) {
260 (*x_vertices)[i + Fstart] = cube_inner_x_vertices[i];
261 (*y_vertices)[i + Fstart] = cube_inner_y_vertices[i];
262 (*z_vertices)[i + Fstart] = cube_inner_z_vertices[i];
265 free(temp_x_vertices);
266 free(temp_y_vertices);
267 free(temp_z_vertices);
268 free(cube_inner_x_vertices);
269 free(cube_inner_y_vertices);
270 free(cube_inner_z_vertices);
272 for (
unsigned i = 0; i < 6; ++i)
273 for (
unsigned j = 0; j < n * n; ++j)
274 (*face_id)[i * n * n + j] = i + 1;
276 unsigned * edge_vertices =
xmalloc(num_edge_coords * 12 *
sizeof(*edge_vertices));
278 Fstart = Estart + 12 * num_edge_coords;
280 for (
unsigned i = 0; i < num_edge_coords * 12; ++i) edge_vertices[i] = i + Estart;
282 unsigned * F =
xmalloc((n + 1) * (n + 1) *
sizeof(*F));
283 unsigned cell_to_vertex_offset = 0;
292 F[0] = 0; F[n] = 2; F[n * (n + 1) + 0] = 1; F[n * (n + 1) + n] = 3;
293 for (
unsigned i = 0; i < num_edge_coords; ++i) {
294 F[i + 1] = edge_vertices[1 * num_edge_coords + i];
295 F[(i + 1) * (n + 1)] = edge_vertices[0 * num_edge_coords + i];
296 F[(i + 1) * (n + 1) + n] = edge_vertices[3 * num_edge_coords + i];
297 F[n * (n + 1) + i + 1] = edge_vertices[2 * num_edge_coords + i];
300 for (
unsigned i = 0; i < n-1; ++i)
301 for (
unsigned j = 0; j < n-1; ++j)
302 F[(i + 1) * (n + 1) + (j + 1)] = Fstart + i * (n - 1) + j;
304 for (
unsigned i = 0; i < n; ++i) {
305 for (
unsigned j = 0; j < n; ++j) {
306 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
307 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
308 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
309 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
320 F[0] = 4; F[n] = 6; F[n * (n + 1) + 0] = 5; F[n * (n + 1) + n] = 7;
321 for (
unsigned i = 0; i < num_edge_coords; ++i) {
322 F[i + 1] = edge_vertices[5 * num_edge_coords + i];
323 F[(i + 1) * (n + 1)] = edge_vertices[4 * num_edge_coords + i];
324 F[(i + 1) * (n + 1) + n] = edge_vertices[7 * num_edge_coords + i];
325 F[n * (n + 1) + i + 1] = edge_vertices[6 * num_edge_coords + i];
328 for (
unsigned i = 0; i < n-1; ++i)
329 for (
unsigned j = 0; j < n-1; ++j)
330 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
332 for (
unsigned i = 0; i < n; ++i) {
333 for (
unsigned j = 0; j < n; ++j) {
334 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
335 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
336 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
337 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
348 F[0] = 0; F[n] = 4; F[n * (n + 1) + 0] = 1; F[n * (n + 1) + n] = 5;
349 for (
unsigned i = 0; i < num_edge_coords; ++i) {
350 F[i + 1] = edge_vertices[8 * num_edge_coords + i];
351 F[(i + 1) * (n + 1)] = edge_vertices[0 * num_edge_coords + i];
352 F[(i + 1) * (n + 1) + n] = edge_vertices[4 * num_edge_coords + i];
353 F[n * (n + 1) + i + 1] = edge_vertices[9 * num_edge_coords + i];
356 for (
unsigned i = 0; i < n-1; ++i)
357 for (
unsigned j = 0; j < n-1; ++j)
358 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
360 for (
unsigned i = 0; i < n; ++i) {
361 for (
unsigned j = 0; j < n; ++j) {
362 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
363 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
364 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
365 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
376 F[0] = 2; F[n] = 6; F[n * (n + 1) + 0] = 3; F[n * (n + 1) + n] = 7;
377 for (
unsigned i = 0; i < num_edge_coords; ++i) {
378 F[i + 1] = edge_vertices[10 * num_edge_coords + i];
379 F[(i + 1) * (n + 1)] = edge_vertices[ 3 * num_edge_coords + i];
380 F[(i + 1) * (n + 1) + n] = edge_vertices[ 7 * num_edge_coords + i];
381 F[n * (n + 1) + i + 1] = edge_vertices[11 * num_edge_coords + i];
384 for (
unsigned i = 0; i < n-1; ++i)
385 for (
unsigned j = 0; j < n-1; ++j)
386 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
388 for (
unsigned i = 0; i < n; ++i) {
389 for (
unsigned j = 0; j < n; ++j) {
390 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
391 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
392 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
393 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
404 F[0] = 0; F[n] = 2; F[n * (n + 1) + 0] = 4; F[n * (n + 1) + n] = 6;
405 for (
unsigned i = 0; i < num_edge_coords; ++i) {
406 F[i + 1] = edge_vertices[ 1 * num_edge_coords + i];
407 F[(i + 1) * (n + 1)] = edge_vertices[ 8 * num_edge_coords + i];
408 F[(i + 1) * (n + 1) + n] = edge_vertices[10 * num_edge_coords + i];
409 F[n * (n + 1) + i + 1] = edge_vertices[ 5 * num_edge_coords + i];
412 for (
unsigned i = 0; i < n-1; ++i)
413 for (
unsigned j = 0; j < n-1; ++j)
414 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
416 for (
unsigned i = 0; i < n; ++i) {
417 for (
unsigned j = 0; j < n; ++j) {
418 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
419 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
420 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
421 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
432 F[0] = 1; F[n] = 3; F[n * (n + 1) + 0] = 5; F[n * (n + 1) + n] = 7;
433 for (
unsigned i = 0; i < num_edge_coords; ++i) {
434 F[i + 1] = edge_vertices[ 2 * num_edge_coords + i];
435 F[(i + 1) * (n + 1)] = edge_vertices[ 9 * num_edge_coords + i];
436 F[(i + 1) * (n + 1) + n] = edge_vertices[11 * num_edge_coords + i];
437 F[n * (n + 1) + i + 1] = edge_vertices[ 6 * num_edge_coords + i];
440 for (
unsigned i = 0; i < n-1; ++i)
441 for (
unsigned j = 0; j < n-1; ++j)
442 F[(i + 1) * (n + 1) + (j + 1)] += (n - 1) * (n - 1);
444 for (
unsigned i = 0; i < n; ++i) {
445 for (
unsigned j = 0; j < n; ++j) {
446 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 0)];
447 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 0)];
448 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 1) * (n + 1) + (j + 1)];
449 (*cell_to_vertex)[cell_to_vertex_offset++] = F[(i + 0) * (n + 1) + (j + 1)];
567 unsigned n,
unsigned * nbr_vertices,
unsigned * nbr_cells,
569 double ** x_vertices,
double ** y_vertices,
double ** x_cells,
571 int ** global_corner_id,
int ** corner_core_mask,
int rank,
int size) {
573 double * x_vertices_3d, * y_vertices_3d, * z_vertices_3d;
578 &x_vertices_3d, &y_vertices_3d,
581 *num_vertices_per_cell =
xmalloc(*nbr_cells *
sizeof(**num_vertices_per_cell));
582 for (
unsigned i = 0; i < *nbr_cells; ++i) (*num_vertices_per_cell)[i] = 4;
585 double * x_cell_3d =
xmalloc(*nbr_cells *
sizeof(*x_cell_3d));
586 double * y_cell_3d =
xmalloc(*nbr_cells *
sizeof(*y_cell_3d));
587 double * z_cell_3d =
xmalloc(*nbr_cells *
sizeof(*z_cell_3d));
589 for (
unsigned i = 0; i < *nbr_cells; ++i) {
591 x_cell_3d[i] = y_cell_3d[i] = z_cell_3d[i] = 0;
593 for (
unsigned j = 0; j < 4; ++j) {
595 x_cell_3d[i] += x_vertices_3d[(*cell_to_vertex)[4 * i + j]];
596 y_cell_3d[i] += y_vertices_3d[(*cell_to_vertex)[4 * i + j]];
597 z_cell_3d[i] += z_vertices_3d[(*cell_to_vertex)[4 * i + j]];
600 double scale = 1.0 / sqrt(x_cell_3d[i] * x_cell_3d[i] +
601 y_cell_3d[i] * y_cell_3d[i] +
602 z_cell_3d[i] * z_cell_3d[i]);
604 x_cell_3d[i] *= scale;
605 y_cell_3d[i] *= scale;
606 z_cell_3d[i] *= scale;
609 int * cell_is_on_rank =
xmalloc(*nbr_cells *
sizeof(*cell_is_on_rank));
614 int * required_vertices =
xcalloc(*nbr_vertices,
sizeof(*required_vertices));
615 int * required_cells =
xcalloc(*nbr_cells,
sizeof(*required_cells));
618 for (
unsigned i = 0, offset = 0; i < *nbr_cells;
619 offset += (*num_vertices_per_cell)[i++]) {
621 if (cell_is_on_rank[i] != rank)
continue;
623 cell_is_on_rank[i] = -1;
625 required_cells[i] = 1;
626 for (
unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j)
631 for (
unsigned i = 0, offset = 0; i < *nbr_cells;
632 offset += (*num_vertices_per_cell)[i++]) {
634 if (!required_cells[i]) {
636 for (
unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
640 required_cells[i] = 1;
649 int * vertex_is_on_rank =
xmalloc(*nbr_vertices *
sizeof(*vertex_is_on_rank));
652 for (
unsigned i = 0; i < *nbr_vertices; ++i)
653 if (required_vertices[i])
654 vertex_is_on_rank[i] = -1;
655 for (
unsigned i = 0, offset = 0; i < *nbr_cells;
656 offset += (*num_vertices_per_cell)[i++]) {
658 if (required_cells[i] && cell_is_on_rank[i] != -1) {
660 for (
unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
664 required_vertices[(*cell_to_vertex)[offset+j]] = 1;
665 vertex_is_on_rank[(*cell_to_vertex)[offset+j]] = cell_is_on_rank[i];
672 int part_num_vertices = 0;
673 int part_num_cells = 0;
674 for (
unsigned i = 0; i < *nbr_vertices; ++i)
675 if (required_vertices[i])
677 for (
unsigned i = 0; i < *nbr_cells; ++i)
678 if(required_cells[i])
681 *global_cell_id =
xmalloc(part_num_cells *
sizeof(**global_cell_id));
683 *global_corner_id =
xmalloc(part_num_vertices *
sizeof(**global_corner_id));
684 *corner_core_mask =
xmalloc(part_num_vertices *
sizeof(**corner_core_mask));
686 *x_cells =
xmalloc(part_num_cells *
sizeof(**x_cells));
687 *y_cells =
xmalloc(part_num_cells *
sizeof(**y_cells));
688 *x_vertices =
xmalloc(part_num_vertices *
sizeof(**x_vertices));
689 *y_vertices =
xmalloc(part_num_vertices *
sizeof(**y_vertices));
692 part_num_vertices = 0;
693 int * global_to_local_vertex =
xmalloc(*nbr_vertices *
sizeof(*global_to_local_vertex));
694 for (
unsigned i = 0; i < *nbr_vertices; ++i) {
696 if (required_vertices[i]) {
698 (*global_corner_id)[part_num_vertices] = i;
699 (*corner_core_mask)[part_num_vertices] = vertex_is_on_rank[i] == -1;
700 double p[3] = {x_vertices_3d[i], y_vertices_3d[i], z_vertices_3d[i]};
701 XYZtoLL(
p, (*x_vertices) + part_num_vertices, (*y_vertices) + part_num_vertices);
702 global_to_local_vertex[i] = part_num_vertices;
707 free(vertex_is_on_rank);
708 *nbr_vertices = part_num_vertices;
709 free(required_vertices);
712 int num_cell_vertex_dependencies = 0;
714 for (
unsigned i = 0, offset = 0; i < *nbr_cells;
715 offset += (*num_vertices_per_cell)[i++]) {
717 if (required_cells[i]) {
719 (*global_cell_id)[part_num_cells] = i;
720 (*cell_core_mask)[part_num_cells] = (cell_is_on_rank[i] == -1);
721 double middle_point[3] = {x_cell_3d[i],y_cell_3d[i],z_cell_3d[i]};
723 for (
unsigned j = 0; j < (*num_vertices_per_cell)[i]; ++j) {
724 (*cell_to_vertex)[num_cell_vertex_dependencies+j] =
725 global_to_local_vertex[(*cell_to_vertex)[offset+j]];
728 XYZtoLL(middle_point, (*x_cells)+part_num_cells, (*y_cells)+part_num_cells);
730 num_cell_vertex_dependencies += (*num_vertices_per_cell)[i];
732 (*num_vertices_per_cell)[part_num_cells] = (*num_vertices_per_cell)[i];
745 *num_vertices_per_cell =
xrealloc(*num_vertices_per_cell, part_num_cells *
746 sizeof(**num_vertices_per_cell));
749 free(cell_is_on_rank);
750 *nbr_cells = part_num_cells;
751 free(required_cells);
752 free(global_to_local_vertex);
760 double * x_vertices, * y_vertices, * z_vertices;
761 unsigned * vertices_of_cell, * dummy;
763 n, &
num_cells, &num_vertices, &x_vertices, &y_vertices, &z_vertices,
764 &vertices_of_cell, &dummy);
768 size_t nv = 4, ne = 4;
770 int dim_cell_id, dim_vertex_id, dim_nv_id, dim_ne_id;
771 int var_vlon_id, var_vlat_id, var_clon_id, var_clat_id, var_mask_id,
772 var_v2c_id, var_c2v_id;
773 yac_nc_create(filename, NC_CLOBBER | NC_64BIT_OFFSET, &ncid);
775 YAC_HANDLE_ERROR(nc_def_dim(ncid,
"vertex", (
size_t)num_vertices, &dim_vertex_id));
779 nc_def_var(ncid,
"vlon", NC_DOUBLE, 1, &dim_vertex_id, &var_vlon_id));
781 nc_def_var(ncid,
"vlat", NC_DOUBLE, 1, &dim_vertex_id, &var_vlat_id));
783 nc_def_var(ncid,
"clon", NC_DOUBLE, 1, &dim_cell_id, &var_clon_id));
785 nc_def_var(ncid,
"clat", NC_DOUBLE, 1, &dim_cell_id, &var_clat_id));
788 ncid,
"cell_sea_land_mask", NC_INT, 1, &dim_cell_id, &var_mask_id));
791 ncid,
"vertex_of_cell", NC_INT, 2, (
int[]){dim_nv_id, dim_cell_id},
795 ncid,
"cells_of_vertex", NC_INT, 2, (
int[]){dim_ne_id, dim_vertex_id},
801 double * lon_buffer =
803 double * lat_buffer =
806 for (
unsigned i = 0; i < num_vertices; ++i) {
807 double coord_3d[3] = {x_vertices[i], y_vertices[i], z_vertices[i]};
808 XYZtoLL(coord_3d, lon_buffer + i, lat_buffer + i);
813 for (
unsigned i = 0; i <
num_cells; ++i) {
814 double coord_3d[3] = {0.0, 0.0, 0.0};
815 for (
unsigned j = 0; j < 4; ++j) {
816 unsigned vertex_idx = vertices_of_cell[4 * i + j];
817 coord_3d[0] += x_vertices[vertex_idx];
818 coord_3d[1] += y_vertices[vertex_idx];
819 coord_3d[2] += z_vertices[vertex_idx];
822 XYZtoLL(coord_3d, lon_buffer + i, lat_buffer + i);
836 for (
unsigned i = 0; i <
num_cells; ++i) int_buffer[i] = 1;
839 for (
unsigned i = 0; i < ne * num_vertices; ++i)
840 int_buffer[i] = INT_MAX;
841 for (
unsigned i = 0; i <
num_cells; ++i) {
842 int cell_id = (int)(i + 1);
843 for (
unsigned j = 0; j < ne; ++j) {
844 unsigned vertex_idx = vertices_of_cell[ne * i + j];
845 for (
unsigned k = 0; k < ne; ++k) {
846 if (int_buffer[k * num_vertices + vertex_idx] == INT_MAX) {
847 int_buffer[k * num_vertices + vertex_idx] = cell_id;
856 for (
unsigned j = 0; j < nv; ++j)
858 vertices_of_cell[nv * i + j] + 1;
862 free(vertices_of_cell);