FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_vis_new_refine.c
Go to the documentation of this file.
1/*****************************************************************************
2 * Copyright (c) 2019 FrontISTR Commons
3 * This software is released under the MIT License, see LICENSE.txt
4 *****************************************************************************/
5
7
8#include <math.h>
10#include "hecmw_vis_comm_util.h"
12#include "hecmw_vis_ray_trace.h"
13#include "hecmw_malloc.h"
14
15/* static void free_headpatch(Head_surfacep_info *head_patch, int nx, int ny,
16 * int nz); */
17static void free_cubehead(Cube_head *cubehead, int voxn);
18
19void transform_face_node(int face, int node[4]) {
20 switch (face) {
21 case 0:
22 node[0] = 0;
23 node[1] = 4;
24 node[2] = 7;
25 node[3] = 3;
26 break;
27 case 1:
28 node[0] = 1;
29 node[1] = 5;
30 node[2] = 6;
31 node[3] = 2;
32 break;
33 case 2:
34 node[0] = 0;
35 node[1] = 1;
36 node[2] = 5;
37 node[3] = 4;
38 break;
39 case 3:
40 node[0] = 3;
41 node[1] = 2;
42 node[2] = 6;
43 node[3] = 7;
44 break;
45 case 4:
46 node[0] = 0;
47 node[1] = 1;
48 node[2] = 2;
49 node[3] = 3;
50 break;
51 case 5:
52 node[0] = 4;
53 node[1] = 5;
54 node[2] = 6;
55 node[3] = 7;
56 break;
57 }
58 return;
59}
60
61static void transform_face_node_351(int face, int node[4]) {
62 switch (face) {
63 case 0:
64 node[0] = 0;
65 node[1] = 1;
66 node[2] = 4;
67 node[3] = 3;
68 break;
69 case 1:
70 node[0] = 1;
71 node[1] = 2;
72 node[2] = 5;
73 node[3] = 4;
74 break;
75 case 2:
76 node[0] = 2;
77 node[1] = 0;
78 node[2] = 3;
79 node[3] = 5;
80 break;
81 case 3:
82 node[0] = 0;
83 node[1] = 2;
84 node[2] = 1;
85 node[3] = 1;
86 break;
87 case 4:
88 node[0] = 3;
89 node[1] = 4;
90 node[2] = 5;
91 node[3] = 5;
92 break;
93 }
94 return;
95}
96
97static void value2_compute(struct hecmwST_local_mesh *mesh, double *node1,
98 double x, double y, double z, int elem_ID,
99 double *field) {
100 int i, j, k;
101 double dist[MAX_N_NODE], coef[MAX_N_NODE];
102 double coord_x, coord_y, coord_z;
103 double v_data[MAX_N_NODE];
104 int nodeID;
105 double sum_coef;
106 int return_flag;
107
108 *field = 0.0;
109 return_flag = 0;
110
111 i = elem_ID;
112 for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
113 j++) {
114 dist[j] = 0.0;
115 nodeID = mesh->elem_node_item[mesh->elem_node_index[i] + j] - 1;
116 coord_x = mesh->node[nodeID * 3];
117 coord_y = mesh->node[nodeID * 3 + 1];
118 coord_z = mesh->node[nodeID * 3 + 2];
119 v_data[j] = node1[nodeID];
120
121 dist[j] = (x - coord_x) * (x - coord_x) + (y - coord_y) * (y - coord_y) +
122 (z - coord_z) * (z - coord_z);
123 if (fabs(dist[j]) < EPSILON) {
124 *field = v_data[j];
125 return_flag = 1;
126 }
127 }
128 if (return_flag == 0) {
129 sum_coef = 0.0;
130 for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
131 j++) {
132 coef[j] = 1.0;
133 for (k = 0; k < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
134 k++) {
135 if (j != k) coef[j] = coef[j] * dist[k];
136 }
137 sum_coef += coef[j];
138 }
139
140 for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
141 j++) {
142 *field = *field + v_data[j] * coef[j] / sum_coef;
143 }
144 }
145 return;
146}
147
148static int judge_inner_voxel_361(double vv[MAX_N_NODE * 3], double point[3]) {
149 double fp[4][3], n_f[3], cen_point[3], sign_f, sign_p, f_cen_point[3], c_c[3],
150 p_c[3];
151 double n_norm, n_c, n_p;
152 int i, j, m, inside, node[4];
153
154 /*find the center point of the element */
155 for (j = 0; j < 3; j++) {
156 cen_point[j] = 0.0;
157 for (i = 0; i < 8; i++) cen_point[j] += vv[i * 3 + j];
158 cen_point[j] /= 8.0;
159 }
160 inside = 1;
161 /*for each face of the element */
162 for (m = 0; m < 6; m++) {
163 transform_face_node(m, node);
164 for (j = 0; j < 3; j++)
165 for (i = 0; i < 4; i++) {
166 fp[i][j] = vv[node[i] * 3 + j];
167 }
168 /*find the normal of the face */
169 n_f[0] = (fp[1][1] - fp[0][1]) * (fp[3][2] - fp[0][2]) -
170 (fp[3][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
171 n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[3][2] - fp[0][2]) +
172 (fp[3][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
173 n_f[2] = (fp[1][0] - fp[0][0]) * (fp[3][1] - fp[0][1]) -
174 (fp[3][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
175 n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
176 if (fabs(n_norm) > EPSILON) {
177 for (j = 0; j < 3; j++) n_f[j] /= n_norm;
178 }
179 /*selce the direction point to inside the element */
180 for (j = 0; j < 3; j++)
181 f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j] + fp[3][j]) / 4.0;
182 for (j = 0; j < 3; j++) c_c[j] = cen_point[j] - f_cen_point[j];
183 n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
184 if (fabs(n_c) > EPSILON) {
185 for (j = 0; j < 3; j++) c_c[j] /= n_c;
186 }
187 sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
188 if (sign_f < -EPSILON) {
189 for (j = 0; j < 3; j++) n_f[j] = -n_f[j];
190 }
191 for (j = 0; j < 3; j++) p_c[j] = point[j] - f_cen_point[j];
192 n_p = sqrt(p_c[0] * p_c[0] + p_c[1] * p_c[1] + p_c[2] * p_c[2]);
193 if (fabs(n_p) > EPSILON) {
194 for (j = 0; j < 3; j++) p_c[j] /= n_p;
195 }
196 sign_p = p_c[0] * n_f[0] + p_c[1] * n_f[1] + p_c[2] * n_f[2];
197 if (sign_p < -EPSILON) {
198 inside = 0;
199 break;
200 }
201 }
202 return (inside);
203}
204
205static int judge_inner_voxel_351(double vv[MAX_N_NODE * 3], double point[3]) {
206 double fp[4][3], n_f[3], cen_point[3], sign_f, sign_p, f_cen_point[3], c_c[3],
207 p_c[3];
208 double n_norm, n_c, n_p;
209 int i, j, m, inside, node[4];
210
211 /*find the center point of the element */
212 for (j = 0; j < 3; j++) {
213 cen_point[j] = 0.0;
214 for (i = 0; i < 6; i++) cen_point[j] += vv[i * 3 + j];
215 cen_point[j] /= 6.0;
216 }
217 inside = 1;
218 /*for each face of the element */
219 for (m = 0; m < 3; m++) {
220 transform_face_node_351(m, node);
221 for (j = 0; j < 3; j++)
222 for (i = 0; i < 4; i++) {
223 fp[i][j] = vv[node[i] * 3 + j];
224 }
225 /*find the normal of the face */
226 n_f[0] = (fp[1][1] - fp[0][1]) * (fp[3][2] - fp[0][2]) -
227 (fp[3][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
228 n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[3][2] - fp[0][2]) +
229 (fp[3][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
230 n_f[2] = (fp[1][0] - fp[0][0]) * (fp[3][1] - fp[0][1]) -
231 (fp[3][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
232 n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
233 if (fabs(n_norm) > EPSILON) {
234 for (j = 0; j < 3; j++) n_f[j] /= n_norm;
235 }
236 /*selce the direction point to inside the element */
237 for (j = 0; j < 3; j++)
238 f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j] + fp[3][j]) / 4.0;
239 for (j = 0; j < 3; j++) c_c[j] = cen_point[j] - f_cen_point[j];
240 n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
241 if (fabs(n_c) > EPSILON) {
242 for (j = 0; j < 3; j++) c_c[j] /= n_c;
243 }
244 sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
245 if (sign_f < -EPSILON) {
246 for (j = 0; j < 3; j++) n_f[j] = -n_f[j];
247 }
248 for (j = 0; j < 3; j++) p_c[j] = point[j] - f_cen_point[j];
249 n_p = sqrt(p_c[0] * p_c[0] + p_c[1] * p_c[1] + p_c[2] * p_c[2]);
250 if (fabs(n_p) > EPSILON) {
251 for (j = 0; j < 3; j++) p_c[j] /= n_p;
252 }
253 sign_p = p_c[0] * n_f[0] + p_c[1] * n_f[1] + p_c[2] * n_f[2];
254 if (sign_p < -EPSILON) {
255 inside = 0;
256 break;
257 }
258 }
259 for (m = 3; m < 5; m++) {
260 for (j = 0; j < 3; j++)
261 for (i = 0; i < 3; i++) {
262 fp[i][j] = vv[node[i] * 3 + j];
263 }
264 /*find the normal of the face */
265 n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
266 (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
267 n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
268 (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
269 n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
270 (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
271 n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
272 if (fabs(n_norm) > EPSILON) {
273 for (j = 0; j < 3; j++) n_f[j] /= n_norm;
274 }
275 /*selce the direction point to inside the element */
276 for (j = 0; j < 3; j++)
277 f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j]) / 3.0;
278 for (j = 0; j < 3; j++) c_c[j] = cen_point[j] - f_cen_point[j];
279 n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
280 if (fabs(n_c) > EPSILON) {
281 for (j = 0; j < 3; j++) c_c[j] /= n_c;
282 }
283 sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
284 if (sign_f < -EPSILON) {
285 for (j = 0; j < 3; j++) n_f[j] = -n_f[j];
286 }
287 for (j = 0; j < 3; j++) p_c[j] = point[j] - f_cen_point[j];
288 n_p = sqrt(p_c[0] * p_c[0] + p_c[1] * p_c[1] + p_c[2] * p_c[2]);
289 if (fabs(n_p) > EPSILON) {
290 for (j = 0; j < 3; j++) p_c[j] /= n_p;
291 }
292 sign_p = p_c[0] * n_f[0] + p_c[1] * n_f[1] + p_c[2] * n_f[2];
293 if (sign_p < -EPSILON) {
294 inside = 0;
295 break;
296 }
297 }
298
299 return (inside);
300}
301
302static int judge_inner_voxel_341(double vv[MAX_N_NODE * 3], double point[3]) {
303 double fp[4][3], n_f[3], cen_point[3], sign_f, sign_p, f_cen_point[3], c_c[3],
304 p_c[3];
305 double n_norm, n_c, n_p;
306 int i, j, m, inside, node[4];
307
308 /*find the center point of the element */
309 for (j = 0; j < 3; j++) {
310 cen_point[j] = 0.0;
311 for (i = 0; i < 4; i++) cen_point[j] += vv[i * 3 + j];
312 cen_point[j] /= 4.0;
313 }
314 inside = 1;
315 /*for each face of the element */
316 for (m = 0; m < 4; m++) {
317 if (m == 0) {
318 node[0] = 0;
319 node[1] = 2;
320 node[2] = 1;
321 } else if (m == 1) {
322 node[0] = 0;
323 node[1] = 1;
324 node[2] = 3;
325 } else if (m == 2) {
326 node[0] = 1;
327 node[1] = 2;
328 node[2] = 3;
329 } else if (m == 3) {
330 node[0] = 2;
331 node[1] = 0;
332 node[2] = 3;
333 }
334
335 for (j = 0; j < 3; j++)
336 for (i = 0; i < 3; i++) {
337 fp[i][j] = vv[node[i] * 3 + j];
338 }
339 /*find the normal of the face */
340 n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
341 (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
342 n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
343 (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
344 n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
345 (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
346 n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
347 if (fabs(n_norm) > EPSILON) {
348 for (j = 0; j < 3; j++) n_f[j] /= n_norm;
349 }
350 /*selce the direction point to inside the element */
351 for (j = 0; j < 3; j++)
352 f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j]) / 3.0;
353 for (j = 0; j < 3; j++) c_c[j] = cen_point[j] - f_cen_point[j];
354 n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
355 if (fabs(n_c) > EPSILON) {
356 for (j = 0; j < 3; j++) c_c[j] /= n_c;
357 }
358 sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
359 if (sign_f < -EPSILON) {
360 for (j = 0; j < 3; j++) n_f[j] = -n_f[j];
361 }
362 for (j = 0; j < 3; j++) p_c[j] = point[j] - f_cen_point[j];
363 n_p = sqrt(p_c[0] * p_c[0] + p_c[1] * p_c[1] + p_c[2] * p_c[2]);
364 if (fabs(n_p) > EPSILON) {
365 for (j = 0; j < 3; j++) p_c[j] /= n_p;
366 }
367 sign_p = p_c[0] * n_f[0] + p_c[1] * n_f[1] + p_c[2] * n_f[2];
368 if (sign_p < -EPSILON) {
369 inside = 0;
370 break;
371 }
372 }
373 return (inside);
374}
375
376/*
377void find2_minmax(In_surface *surface, int i, double minmax[6])
378{
379 int j;
380
381
382 for(j=0;j<3;j++)
383 minmax[j*2]=minmax[j*2+1]=surface->vertex[(surface->patch[i*3]-1)*3+j];
384 for(j=0;j<3;j++) {
385 if(surface->vertex[(surface->patch[i*3+1]-1)*3+j]<minmax[j*2])
386 minmax[j*2]=surface->vertex[(surface->patch[i*3+1]-1)*3+j];
387 if(surface->vertex[(surface->patch[i*3+1]-1)*3+j]>minmax[j*2+1])
388 minmax[j*2+1]=surface->vertex[(surface->patch[i*3+1]-1)*3+j];
389 if(surface->vertex[(surface->patch[i*3+2]-1)*3+j]<minmax[j*2])
390 minmax[j*2]=surface->vertex[(surface->patch[i*3+2]-1)*3+j];
391 if(surface->vertex[(surface->patch[i*3+2]-1)*3+j]>minmax[j*2+1])
392 minmax[j*2+1]=surface->vertex[(surface->patch[i*3+2]-1)*3+j];
393 }
394
395 return;
396}
397
398 */
399
400void refinement(struct hecmwST_local_mesh *mesh, double *node1, int n_voxel,
401 double *voxel_dxyz, double *voxel_orig_xyz, int *level,
403 double *extent, int my_rank, HECMW_Comm VIS_COMM,
404 int *empty_flag, double *var)
405
406{
407 int i, j;
408
409 int pe_size;
410 Cube_head *cubehead;
411 Cube_point *p1;
412 int m1, m2, n1, n2, n3, m;
413 int nx, ny, nz, flag;
414 int i1, j1, k1, i2, j2, k2;
415 double dx, dy, dz, minx, miny, minz, maxx, maxy, maxz, x, y, z;
416 double field;
417 int *code;
418 double *value;
419 int vox_id, point_num;
420 int pe_id;
421 HECMW_Status stat;
422 double vv[MAX_N_NODE * 3], in_point[3];
423 int inside, nodeID;
424
425 HECMW_Comm_size(VIS_COMM, &pe_size);
426 cubehead = (Cube_head *)HECMW_calloc(n_voxel, sizeof(Cube_head));
427 if (cubehead == NULL) HECMW_vis_memory_exit("In refinement: cubehead");
428 for (i = 0; i < n_voxel; i++) {
429 cubehead[i].point_num = 0;
430 cubehead[i].cube_link = NULL;
431 }
432
433 for (m1 = 0; m1 < n_voxel; m1++)
434 for (m2 = 0; m2 < voxel_n_neighbor_pe[m1]; m2++) {
435 if (voxel_neighbor_pe[m1][m2] == my_rank) {
436 nx = level[m1 * 3];
437 ny = level[m1 * 3 + 1];
438 nz = level[m1 * 3 + 2];
439 dx = voxel_dxyz[m1 * 3] / (double)nx;
440 dy = voxel_dxyz[m1 * 3 + 1] / (double)ny;
441 dz = voxel_dxyz[m1 * 3 + 2] / (double)nz;
442 /* fprintf(stderr, "in voxel %d: nx, ny, nz=%d %d %d dx,
443 * dy, dz=%lf %lf %lf\n", m1,nx, ny, nz, dx, dy, dz);
444 */
445 minx = voxel_orig_xyz[m1 * 3];
446 maxx = voxel_orig_xyz[m1 * 3] + voxel_dxyz[m1 * 3];
447 miny = voxel_orig_xyz[m1 * 3 + 1];
448 maxy = voxel_orig_xyz[m1 * 3 + 1] + voxel_dxyz[m1 * 3 + 1];
449 minz = voxel_orig_xyz[m1 * 3 + 2];
450 maxz = voxel_orig_xyz[m1 * 3 + 2] + voxel_dxyz[m1 * 3 + 2];
451 for (i = 0; i < mesh->n_elem; i++) {
452 if (mesh->elem_type[i] < 400) {
453 flag = 1;
454 if ((extent[i * 6] > maxx) || (extent[i * 6 + 1] < minx) ||
455 (extent[i * 6 + 2] > maxy) || (extent[i * 6 + 3] < miny) ||
456 (extent[i * 6 + 4] > maxz) || (extent[i * 6 + 5] < minz))
457 flag = 0;
458 if (flag == 1) {
459 i1 = (int)((extent[i * 6] - EPSILON - minx) / dx);
460 j1 = (int)((extent[i * 6 + 2] - EPSILON - miny) / dy);
461 k1 = (int)((extent[i * 6 + 4] - EPSILON - minz) / dz);
462 i2 = (int)((extent[i * 6 + 1] + EPSILON - minx) / dx);
463 j2 = (int)((extent[i * 6 + 3] + EPSILON - miny) / dy);
464 k2 = (int)((extent[i * 6 + 5] + EPSILON - minz) / dz);
465 if (i1 < 0) i1 = 0;
466 if (i1 > nx) i1 = i2 + 1;
467 if (j1 < 0) j1 = 0;
468 if (j1 > ny) j1 = j2 + 1;
469 if (k1 < 0) k1 = 0;
470 if (k1 > nz) k1 = k2 + 1;
471 if (i2 < 0) i2 = i1 - 1;
472 if (i2 > nx) i2 = nx;
473 if (j2 < 0) j2 = j1 - 1;
474 if (j2 > ny) j2 = ny;
475 if (k2 < 0) k2 = k1 - 1;
476 if (k2 > nz) k2 = nz;
477 for (n3 = k1; n3 <= k2; n3++)
478 for (n2 = j1; n2 <= j2; n2++)
479 for (n1 = i1; n1 <= i2; n1++) {
480 x = minx + dx * n1;
481 y = miny + dy * n2;
482 z = minz + dz * n3;
483 /* if((x>=extent[i].x_min-EPSILON) &&
484 (x<=extent[i].x_max+EPSILON) &&
485 (y>=extent[i].y_min-EPSILON) && (y<=extent[i].y_max+EPSILON) &&
486 (z>=extent[i].z_min-EPSILON) && (z<=extent[i].z_max+EPSILON)) {
487 */
488 in_point[0] = x;
489 in_point[1] = y;
490 in_point[2] = z;
491 for (m = mesh->elem_node_index[i];
492 m < mesh->elem_node_index[i + 1]; m++) {
493 nodeID = mesh->elem_node_item[m] - 1;
494 vv[(m - mesh->elem_node_index[i]) * 3] =
495 mesh->node[nodeID * 3];
496 vv[(m - mesh->elem_node_index[i]) * 3 + 1] =
497 mesh->node[nodeID * 3 + 1];
498 vv[(m - mesh->elem_node_index[i]) * 3 + 2] =
499 mesh->node[nodeID * 3 + 2];
500 }
501 if ((mesh->elem_type[i] == 361) ||
502 (mesh->elem_type[i] == 362))
503 inside = judge_inner_voxel_361(vv, in_point);
504 else if ((mesh->elem_type[i] == 341) ||
505 (mesh->elem_type[i] == 342))
506 inside = judge_inner_voxel_341(vv, in_point);
507 else if ((mesh->elem_type[i] == 351) ||
508 (mesh->elem_type[i] == 352))
509 inside = judge_inner_voxel_351(vv, in_point);
510
511 /*------------need add judge inner for tetrahedras later
512!!!!!!!!!!!!!!!!!!!!
513-----------------*/
514 if (inside == 1) {
515 cubehead[m1].point_num++;
516 p1 = (Cube_point *)HECMW_malloc(sizeof(Cube_point));
517 if (p1 == NULL) HECMW_vis_memory_exit("new_refine: p1");
518 p1->code[0] = n1;
519 p1->code[1] = n2;
520 p1->code[2] = n3;
521
522 value2_compute(mesh, node1, x, y, z, i, &field);
523 p1->field = field;
524 p1->next_point = cubehead[m1].cube_link;
525 cubehead[m1].cube_link = p1;
526 } /*end if(xyz) */
527 } /*end of n1, n2, n3 loop */
528 } /*end of flag=1 */
529 } /*end of element type */
530 } /* end of element loop */
531 } /* end of if */
532 } /* end of m1, m2 loop */
533 /* pe_vox=(int *)HECMW_calloc(pe_size, sizeof(int));
534 if(pe_vox==NULL) {
535 fprintf(stderr, "There is no enough memory for pe_vox\n");
536 exit(0);
537 }
538 pe_vox[0]=7; pe_vox[1]=6; pe_vox[2]=5; pe_vox[3]=4; pe_vox[4]=3;
539 pe_vox[5]=2; pe_vox[6]=1; pe_vox[7]=0;
540 */
541 HECMW_Barrier(VIS_COMM);
542
543 /* start send cube points for each voxel in this PE */
544
545 /* if (nflag == 0) {
546if ((sta1 = (HECMW_Status *)HECMW_calloc(HECMW_STATUS_SIZE,
547sizeof(HECMW_Status)))
548 == NULL) {
549fprintf(stderr, "Not enough memory\n");
550exit(1);
551}
552if ((sta2 = (HECMW_Status *)HECMW_calloc(HECMW_STATUS_SIZE,
553sizeof(HECMW_Status)))
554 == NULL) {
555fprintf(stderr, "Not enough memory\n");
556exit(1);
557}
558if ((req1 = (HECMW_Request *)HECMW_calloc(pe_size-1, sizeof(HECMW_Request)))
559 == NULL) {
560fprintf(stderr, "Not enough memory\n");
561exit(1);
562}
563if ((req2 = (HECMW_Request *)HECMW_calloc(pe_size-1, sizeof(HECMW_Request)))
564 == NULL) {
565fprintf(stderr, "Not enough memory\n");
566exit(1);
567}
568nflag = 1;
569}
570 */
571 for (m1 = 0; m1 < n_voxel; m1++) {
572 pe_id = m1 % pe_size;
573 /* fprintf(stderr, "m1=%d pe_size=%d\n", m1, pe_size);
574 */
575 /* pe_id=pe_vox[m1];
576 */
577 if (pe_id != my_rank) {
578 if (cubehead[m1].point_num == 0) {
579 HECMW_Send(&m1, 1, HECMW_INT, pe_id, 0, VIS_COMM);
580 HECMW_Send(&(cubehead[m1].point_num), 1, HECMW_INT, pe_id, 0, VIS_COMM);
581
582 /* HECMW_Isend(&m1,1,HECMW_INT, pe_id, 0,
583VIS_COMM,&req1[pe_id]);
584HECMW_Isend(&(cubehead[m1].point_num),1,HECMW_INT, pe_id, 0,
585VIS_COMM,&req1[pe_id]);
586 */
587 } else if (cubehead[m1].point_num > 0) {
588 code = (int *)HECMW_calloc(3 * cubehead[m1].point_num, sizeof(int));
589 value = (double *)HECMW_calloc(cubehead[m1].point_num, sizeof(double));
590 if ((code == NULL) || (value == NULL))
591 HECMW_vis_memory_exit("new_refine: code, value");
592 p1 = cubehead[m1].cube_link;
593 for (i = 0; i < cubehead[m1].point_num; i++) {
594 for (j = 0; j < 3; j++) code[i * 3 + j] = p1->code[j];
595 value[i] = p1->field;
596 p1 = p1->next_point;
597 }
598 HECMW_Send(&m1, 1, HECMW_INT, pe_id, 0, VIS_COMM);
599 HECMW_Send(&(cubehead[m1].point_num), 1, HECMW_INT, pe_id, 0, VIS_COMM);
600 HECMW_Send(code, 3 * cubehead[m1].point_num, HECMW_INT, pe_id, 0,
601 VIS_COMM);
602 HECMW_Send(value, cubehead[m1].point_num, HECMW_DOUBLE, pe_id, 0,
603 VIS_COMM);
604
605 /* HECMW_Isend(&m1, 1, HECMW_INT, pe_id, 0, VIS_COMM,&req1[pe_id]);
606 HECMW_Isend(&(cubehead[m1].point_num),1,HECMW_INT, pe_id, 0, VIS_COMM,&req1[pe_id]);
607 HECMW_Isend(code,3*cubehead[m1].point_num,HECMW_INT, pe_id, 0, VIS_COMM,&req1[pe_id]);
608 HECMW_Isend(value, cubehead[m1].point_num,HECMW_DOUBLE, pe_id, 0, VIS_COMM,&req1[pe_id]);
609 HECMW_Isend(gradient, 3*cubehead[m1].point_num,HECMW_DOUBLE, pe_id, 0, VIS_COMM,&req1[pe_id]);
610 */ HECMW_free(code);
611 HECMW_free(value);
612 } /*end of else ifnum>0*/
613 /* fprintf(stderr, "on pe %d finish send vox %d to pe
614 * %d\n", my_rank, m1, pe_id);
615 */
616 } /* end of if rank!=mype */
617 else if (pe_id == my_rank) {
618 /* max_level = vox->info[m1].level[0];
619 if (max_level < vox->info[m1].level[1])
620 max_level = vox->info[m1].level[1];
621 if (max_level < vox->info[m1].level[2])
622 max_level = vox->info[m1].level[2];
623 */
624 /* nx=ny=nz= (int)pow(2.0,(double)max_level);
625 */
626 nx = level[m1 * 3];
627 ny = level[m1 * 3 + 1];
628 nz = level[m1 * 3 + 2];
629 /* if(mynode==0)
630fprintf(stderr, "Final nx, ny, nz in refinement is %d %d %d\n", nx, ny, nz);
631 */
632
633 for (i = 0; i < (nx + 1) * (ny + 1) * (nz + 1); i++) empty_flag[i] = 0;
634 /*first copy the part result on this pe into result file */
635 if (cubehead[m1].point_num > 0) {
636 p1 = cubehead[m1].cube_link;
637 for (i = 0; i < cubehead[m1].point_num; i++) {
638 i1 = p1->code[0];
639 j1 = p1->code[1];
640 k1 = p1->code[2];
641 if ((k1 * (ny + 1) * (nx + 1) + j1 * (nx + 1) + i1 < 0) ||
642 (k1 * (ny + 1) * (nx + 1) + j1 * (nx + 1) + i1 >
643 (nx + 1) * (ny + 1) * (nz + 1))) {
644 fprintf(stderr, "There is wrong with code code=%d %d %d\n", i1, j1,
645 k1);
647 "ERROR: HEC-MW-VIS-E2004: voxel search error in new_refine");
648 }
649 /* result[(k1*(ny+1)*(nx+1)+j1*(nx+1)+i1)*4]=p1->field;
650 */ var[k1 * (ny + 1) * (nx + 1) +
651 j1 * (nx + 1) + i1] = p1->field;
652 /*
653 for(j=0;j<3;j++)
654 result[(k1*(ny+1)*(nx+1)+j1*(nx+1)+i1)*4+j+1]=p1->grad[j];
655 empty_flag[k1*(ny+1)*(nx+1)+j1*(nx+1)+i1]=1;
656 */
657 empty_flag[k1 * (ny + 1) * (nx + 1) + j1 * (nx + 1) + i1] = 1;
658 p1 = p1->next_point;
659 }
660 }
661 /*second copy the part result from other pes into result file */
662
663 for (i = 0; i < pe_size; i++)
664 if (i != my_rank) {
665 HECMW_Recv(&vox_id, 1, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM, &stat);
666 /*
667 HECMW_Irecv(&vox_id, 1, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM, &req2[i]);
668 */ /*
669 fprintf(stderr, "vox=%d pe=%d from pe
670 %d\n", m1, my_rank, i);
671 */
672 if (vox_id != m1)
674 "ERROR: HEC-MW-VIS-E2005: data communication error in "
675 "new_refine");
676 HECMW_Recv(&point_num, 1, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM,
677 &stat);
678
679 /* HECMW_Irecv(&point_num, 1, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM, &req2[i]);
680 */ if (point_num > 0) {
681 code = (int *)HECMW_calloc(3 * point_num, sizeof(int));
682 value = (double *)HECMW_calloc(point_num, sizeof(double));
683 if ((code == NULL) || (value == NULL))
684 HECMW_vis_memory_exit("new_refine: code, value");
685 HECMW_Recv(code, point_num * 3, HECMW_INT, i, HECMW_ANY_TAG,
686 VIS_COMM, &stat);
687 HECMW_Recv(value, point_num, HECMW_DOUBLE, i, HECMW_ANY_TAG,
688 VIS_COMM, &stat);
689 /*
690 HECMW_Irecv(code, point_num*3, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM, &req2[i]);
691 HECMW_Irecv(value, point_num, HECMW_DOUBLE, i, HECMW_ANY_TAG, VIS_COMM,&req2[i]);
692 HECMW_Irecv(gradient, point_num*3, HECMW_DOUBLE, i, HECMW_ANY_TAG, VIS_COMM, &req2[i]);
693 */ for (j = 0; j < point_num;
694 j++) {
695 i1 = code[j * 3];
696 j1 = code[j * 3 + 1];
697 k1 = code[j * 3 + 2];
698 /* result[(k1*(ny+1)*(nx+1)+j1*(nx+1)+i1)*4]=value[j];
699for(k=0;k<3;k++)
700 result[(k1*(ny+1)*(nx+1)+j1*(nx+1)+i1)*4+k+1]=gradient[j*3+k];
701empty_flag[k1*(ny+1)*(nx+1)+j1*(nx+1)+i1]=1;
702 */
703 var[k1 * (ny + 1) * (nx + 1) + j1 * (nx + 1) + i1] = value[j];
704 empty_flag[k1 * (ny + 1) * (nx + 1) + j1 * (nx + 1) + i1] = 1;
705 }
706 HECMW_free(code);
707 HECMW_free(value);
708 } /* end of ifpoint_num>0*/
709 } /* end of pe loop */
710 /* fprintf(stderr, "Finish to
711 * receive\n");
712 */
713 /* third, organize to octree */
714 /*
715HECMW_free(sta1);
716HECMW_free(sta2);
717HECMW_free(req1);
718HECMW_free(req2);
719 */
720 /* fourth, output result file */
721 /* sprintf(filename, "%s-%d.%d", outfile, node->iter, m1);
722if ((fp = fopen(filename, "w")) == NULL) {
723fprintf(stderr, "output file error\n");
724exit(1);
725}
726fprintf(fp, " %lf %lf %lf\n", vox->info[m1].dx, vox->info[m1].dy,
727vox->info[m1].dz);
728fprintf(fp, " %d\n", max_level);
729fprintf(fp, " %d\n", cont->n_var);
730for (i = 0; i < cont->n_var; i++) {
731fprintf(fp, " %s\n", cont->name[i]);
732}
733fprintf(fp, " %d %d %d\n", 1, 1, 1);
734fprintf(fp, " %lf %lf %lf\n", vox->info[m1].orig_x, vox->info[m1].orig_y,
735vox->info[m1].orig_z);
736fprintf(fp, "%d %d %d\n", nx, ny, nz);
737for(i=0;i<(nx+1)*(ny+1)*(nz+1);i++) {
738fprintf(fp, "%d\n", empty_flag[i]);
739if(empty_flag[i]==0)
740fprintf(fp, "%d\n", i);
741if(empty_flag[i]==1)
742fprintf(fp, "%lf %lf %lf %lf\n", result[i*4], result[i*4+1], result[i*4+2],
743result[i*4+3]);
744}
745fclose(fp);
746 */
747
748 } /*end of if=my_rank*/
749 HECMW_Barrier(VIS_COMM);
750 }
751
752 /* HECMW_free application memory */
753 free_cubehead(cubehead, n_voxel);
754
755 return;
756}
757
758#if 0
759static void free_headpatch(Head_surfacep_info *head_patch, int nx, int ny, int nz)
760{
761 int i;
762 Surfacep_info *p1, *p2;
763
764 for(i=0;i<nx*ny*nz;i++) {
765 if(head_patch[i].num_of_patch>0) {
766 p1=head_patch[i].next_patch;
767 while(p1!=NULL) {
768 p2=p1;
769 p1=p1->next_patch;
770 HECMW_free(p2);
771 }
772 }
773 }
774 HECMW_free(head_patch);
775 return;
776}
777#endif
778
779static void free_cubehead(Cube_head *cubehead, int voxn) {
780 int i;
781 Cube_point *p1, *p2;
782
783 for (i = 0; i < voxn; i++) {
784 if (cubehead[i].point_num > 0) {
785 p1 = cubehead[i].cube_link;
786 while (p1 != NULL) {
787 p2 = p1;
788 p1 = p1->next_point;
789 HECMW_free(p2);
790 }
791 }
792 }
793 HECMW_free(cubehead);
794 return;
795}
int HECMW_Comm_size(HECMW_Comm comm, int *size)
Definition: hecmw_comm.c:37
int HECMW_Send(void *buffer, int count, HECMW_Datatype datatype, int dest, int tag, HECMW_Comm comm)
Definition: hecmw_comm.c:193
int HECMW_Barrier(HECMW_Comm comm)
Definition: hecmw_comm.c:95
int HECMW_Recv(void *buffer, int count, HECMW_Datatype datatype, int source, int tag, HECMW_Comm comm, HECMW_Status *status)
Definition: hecmw_comm.c:235
#define HECMW_INT
Definition: hecmw_config.h:48
MPI_Status HECMW_Status
Definition: hecmw_config.h:36
#define HECMW_DOUBLE
Definition: hecmw_config.h:50
MPI_Comm HECMW_Comm
Definition: hecmw_config.h:30
int HECMW_ANY_TAG
struct hecmwST_local_mesh * mesh
Definition: hecmw_repart.h:71
#define NULL
#define HECMW_calloc(nmemb, size)
Definition: hecmw_malloc.h:21
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
#define HECMW_malloc(size)
Definition: hecmw_malloc.h:20
#define EPSILON
void HECMW_vis_print_exit(char *var)
void HECMW_vis_memory_exit(char *var)
void refinement(struct hecmwST_local_mesh *mesh, double *node1, int n_voxel, double *voxel_dxyz, double *voxel_orig_xyz, int *level, int *voxel_n_neighbor_pe, int **voxel_neighbor_pe, double *extent, int my_rank, HECMW_Comm VIS_COMM, int *empty_flag, double *var)
void transform_face_node(int face, int node[4])
int ** voxel_neighbor_pe
int * voxel_n_neighbor_pe
double * voxel_dxyz
double * extent
double * voxel_orig_xyz
int * level
#define MAX_N_NODE
struct cube_pointer_struct * cube_link
struct cube_pointer_struct * next_point
Surfacep_info * next_patch
struct surfacep_info * next_patch