FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 
6 #include "hecmw_vis_new_refine.h"
7 
8 #include <math.h>
9 #include "hecmw_vis_mem_util.h"
10 #include "hecmw_vis_comm_util.h"
11 #include "hecmw_vis_resampling.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); */
17 static void free_cubehead(Cube_head *cubehead, int voxn);
18 
19 void 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 
61 static 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 
97 static 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 
148 static 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 
205 static 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 
302 static 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 /*
377 void 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 
400 void refinement(struct hecmwST_local_mesh *mesh, double *node1, int n_voxel,
401  double *voxel_dxyz, double *voxel_orig_xyz, int *level,
402  int *voxel_n_neighbor_pe, int **voxel_neighbor_pe,
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) {
546 if ((sta1 = (HECMW_Status *)HECMW_calloc(HECMW_STATUS_SIZE,
547 sizeof(HECMW_Status)))
548  == NULL) {
549 fprintf(stderr, "Not enough memory\n");
550 exit(1);
551 }
552 if ((sta2 = (HECMW_Status *)HECMW_calloc(HECMW_STATUS_SIZE,
553 sizeof(HECMW_Status)))
554  == NULL) {
555 fprintf(stderr, "Not enough memory\n");
556 exit(1);
557 }
558 if ((req1 = (HECMW_Request *)HECMW_calloc(pe_size-1, sizeof(HECMW_Request)))
559  == NULL) {
560 fprintf(stderr, "Not enough memory\n");
561 exit(1);
562 }
563 if ((req2 = (HECMW_Request *)HECMW_calloc(pe_size-1, sizeof(HECMW_Request)))
564  == NULL) {
565 fprintf(stderr, "Not enough memory\n");
566 exit(1);
567 }
568 nflag = 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,
583 VIS_COMM,&req1[pe_id]);
584 HECMW_Isend(&(cubehead[m1].point_num),1,HECMW_INT, pe_id, 0,
585 VIS_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)
630 fprintf(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);
646  HECMW_vis_print_exit(
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)
673  HECMW_vis_print_exit(
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];
699 for(k=0;k<3;k++)
700  result[(k1*(ny+1)*(nx+1)+j1*(nx+1)+i1)*4+k+1]=gradient[j*3+k];
701 empty_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  /*
715 HECMW_free(sta1);
716 HECMW_free(sta2);
717 HECMW_free(req1);
718 HECMW_free(req2);
719  */
720  /* fourth, output result file */
721  /* sprintf(filename, "%s-%d.%d", outfile, node->iter, m1);
722 if ((fp = fopen(filename, "w")) == NULL) {
723 fprintf(stderr, "output file error\n");
724 exit(1);
725 }
726 fprintf(fp, " %lf %lf %lf\n", vox->info[m1].dx, vox->info[m1].dy,
727 vox->info[m1].dz);
728 fprintf(fp, " %d\n", max_level);
729 fprintf(fp, " %d\n", cont->n_var);
730 for (i = 0; i < cont->n_var; i++) {
731 fprintf(fp, " %s\n", cont->name[i]);
732 }
733 fprintf(fp, " %d %d %d\n", 1, 1, 1);
734 fprintf(fp, " %lf %lf %lf\n", vox->info[m1].orig_x, vox->info[m1].orig_y,
735 vox->info[m1].orig_z);
736 fprintf(fp, "%d %d %d\n", nx, ny, nz);
737 for(i=0;i<(nx+1)*(ny+1)*(nz+1);i++) {
738 fprintf(fp, "%d\n", empty_flag[i]);
739 if(empty_flag[i]==0)
740 fprintf(fp, "%d\n", i);
741 if(empty_flag[i]==1)
742 fprintf(fp, "%lf %lf %lf %lf\n", result[i*4], result[i*4+1], result[i*4+2],
743 result[i*4+3]);
744 }
745 fclose(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
759 static 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 
779 static 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 }
hecmw_malloc.h
hecmw_vis_mem_util.h
HECMW_DOUBLE
#define HECMW_DOUBLE
Definition: hecmw_config.h:50
cube_head_struct::cube_link
struct cube_pointer_struct * cube_link
Definition: hecmw_vis_resampling.h:83
cube_head_struct
Definition: hecmw_vis_resampling.h:81
HECMW_Send
int HECMW_Send(void *buffer, int count, HECMW_Datatype datatype, int dest, int tag, HECMW_Comm comm)
Definition: hecmw_comm.c:193
hecmw_vis_ray_trace.h
hecmwST_local_mesh::elem_node_item
int * elem_node_item
Definition: hecmw_struct.h:196
HECMW_malloc
#define HECMW_malloc(size)
Definition: hecmw_malloc.h:20
hecmw_vis_resampling.h
mesh
struct hecmwST_local_mesh * mesh
Definition: hecmw_repart.h:71
refinement
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)
Definition: hecmw_vis_new_refine.c:400
surfacep_info
Definition: hecmw_vis_resampling.h:28
hecmwST_local_mesh
Definition: hecmw_struct.h:139
voxel_n_neighbor_pe
int * voxel_n_neighbor_pe
Definition: hecmw_vis_pvr_main.c:26
hecmwST_local_mesh::elem_type
int * elem_type
Definition: hecmw_struct.h:191
HECMW_vis_print_exit
void HECMW_vis_print_exit(char *var)
Definition: hecmw_vis_mem_util.c:21
hecmwST_local_mesh::n_elem
int n_elem
Definition: hecmw_struct.h:184
voxel_dxyz
double * voxel_dxyz
Definition: hecmw_vis_pvr_main.c:25
level
int * level
Definition: hecmw_vis_pvr_main.c:26
HECMW_calloc
#define HECMW_calloc(nmemb, size)
Definition: hecmw_malloc.h:21
voxel_orig_xyz
double * voxel_orig_xyz
Definition: hecmw_vis_pvr_main.c:25
hecmwST_local_mesh::node
double * node
Definition: hecmw_struct.h:170
HECMW_vis_memory_exit
void HECMW_vis_memory_exit(char *var)
Definition: hecmw_vis_mem_util.c:12
extent
double * extent
Definition: hecmw_vis_pvr_main.c:27
HECMW_Status
MPI_Status HECMW_Status
Definition: hecmw_config.h:36
hecmw_vis_new_refine.h
head_surface_info
Definition: hecmw_vis_resampling.h:34
hecmw_vis_comm_util.h
hecmwST_local_mesh::elem_node_index
int * elem_node_index
Definition: hecmw_struct.h:195
HECMW_ANY_TAG
int HECMW_ANY_TAG
Definition: hecmw_dlb_comm_util.h:7
cube_pointer_struct::field
double field
Definition: hecmw_vis_resampling.h:75
voxel_neighbor_pe
int ** voxel_neighbor_pe
Definition: hecmw_vis_pvr_main.c:26
surfacep_info::next_patch
struct surfacep_info * next_patch
Definition: hecmw_vis_resampling.h:31
cube_pointer_struct
Definition: hecmw_vis_resampling.h:73
HECMW_Comm
MPI_Comm HECMW_Comm
Definition: hecmw_config.h:30
HECMW_Comm_size
int HECMW_Comm_size(HECMW_Comm comm, int *size)
Definition: hecmw_comm.c:37
EPSILON
#define EPSILON
Definition: hecmw_vis_color_mapping.c:10
HECMW_INT
#define HECMW_INT
Definition: hecmw_config.h:48
cube_pointer_struct::code
int code[3]
Definition: hecmw_vis_resampling.h:74
HECMW_Recv
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
NULL
#define NULL
Definition: hecmw_io_nastran.c:30
head_surface_info::next_patch
Surfacep_info * next_patch
Definition: hecmw_vis_resampling.h:36
cube_head_struct::point_num
int point_num
Definition: hecmw_vis_resampling.h:82
HECMW_free
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
transform_face_node
void transform_face_node(int face, int node[4])
Definition: hecmw_vis_new_refine.c:19
cube_pointer_struct::next_point
struct cube_pointer_struct * next_point
Definition: hecmw_vis_resampling.h:77
HECMW_Barrier
int HECMW_Barrier(HECMW_Comm comm)
Definition: hecmw_comm.c:95
MAX_N_NODE
#define MAX_N_NODE
Definition: hecmw_vis_ray_trace.h:54