1 /*****************************************************************************
2 * Copyright (c) 2019 FrontISTR Commons
3 * This software is released under the MIT License, see LICENSE.txt
4 *****************************************************************************/
5 /*----------------------------------------------------------------------
6 * Subroutines in this file on isosurface generation for hexahedra by Marching
7 *Cubes is based
8 * on the revision of Dr. Yuriko Takeshima's codes when she was working part
9 *time in RIST
10 *----------------------------------------------------------------------*/
11
12 #include "hecmw_vis_surface_compute.h"
13
14 #include <math.h>
15 #include "hecmw_vis_mem_util.h"
16 #include "hecmw_vis_case_table.h"
17 #include "hecmw_vis_tetra_intersect.h"
18 #include "hecmw_vis_patch_const.h"
19 #include "hecmw_malloc.h"
20
HECMW_vis_surface_compute(Surface * sff,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int * bdflag,int * sum_v,int * sum_t,int * tvertex,int * tpatch,double * minc,double * maxc,Result * result,int sf_i,int mynode,HECMW_Comm VIS_COMM)21 int HECMW_vis_surface_compute(Surface *sff, struct hecmwST_local_mesh *mesh,
22 struct hecmwST_result_data *data, int *bdflag,
23 int *sum_v, int *sum_t, int *tvertex, int *tpatch,
24 double *minc, double *maxc, Result *result,
25 int sf_i, int mynode, HECMW_Comm VIS_COMM)
26
27 {
28 /* in_volume */
29 int i, j, k, mm;
30 Point **CS_verts_head;
31 Polygon **CS_polys_head;
32 int alpha_index, beta_index;
33 int sum_verts;
34 int sum_table;
35 Point **CS_verts_tail, **CS_verts_refer;
36 Polygon **CS_polys_tail;
37 Cube_polygons alpha_cube, beta_cube;
38
39 int n_elem, tmp_int;
40 Cell *cell;
41 Tetra *tetra;
42
43 int disamb_flag;
44
45 Polygon *CS_polys, *CS_polys_tmp;
46 Point *CS_verts, *CS_verts_tmp;
47 int sum_polys, poly_num;
48 int aplist_size, bplist_size, cplist_size;
49 int num_nh_verts,
50 num_nh_patch; /* vertex on patches generated in non-hexahedra */
51 int flag_hexa, flag_tetra;
52 Hash_vertex *vertex_hash_table;
53 Tetra_point *tetra_point;
54 Head_patch_tetra *head_patch_tetra;
55 Tetra_point *p1, *p2;
56 Patch_tetra *t1, *t2;
57 Hash_vertex *h1, *h2;
58 int minus_patch;
59 double tmp_data[9];
60 int tn_component, c_base;
61 double tmp;
62
63 sum_verts = 0;
64 num_nh_verts = 0;
65 num_nh_patch = 0;
66 sum_polys = 0;
67
68 n_elem = mesh->n_elem;
69 tn_component = 0;
70 for (i = 0; i < data->nn_component; i++) tn_component += data->nn_dof[i];
71
72 /* prism = (Prism *)HECMW_malloc(sizeof(Prism));
73 */
74 if (mesh->elem_type[0] > 300) {
75 flag_hexa = 0;
76 flag_tetra = 0;
77
78 disamb_flag = 1;
79 /* fprintf(stderr, "sff inf: color=%d color_sub=%d con=%lf %lf %lf% lf\n",
80 sff->color_comp, sff->color_subcomp, sff->cont_equ[6],
81 sff->cont_equ[7], sff->cont_equ[8],sff->cont_equ[9]);
82
83 sprintf(test_file, "test_ucd.%d.%d.inp", time_step, mynode);
84 write_mesh_display(test_file, mesh, data);
85 HECMW_Barrier(VIS_COMM);
86 */
87 /*
88 for(i = 0; i < mesh->ne_internal; i++) {
89 tmp_int=mesh->elem_internal_list[i]-1;
90 */
91 for (i = 0; i < mesh->n_elem; i++)
92 if (mesh->elem_ID[i * 2 + 1] == mynode) {
93 tmp_int = i;
94 if ((mesh->elem_type[tmp_int] == 361) ||
95 (mesh->elem_type[tmp_int] == 362) ||
96 (mesh->elem_type[tmp_int] == 351) ||
97 (mesh->elem_type[tmp_int] == 352)) {
98 if (flag_hexa == 0) {
99 flag_hexa = 1;
100
101 CS_polys_head = (Polygon **)HECMW_malloc(sizeof(Polygon *));
102 CS_verts_head = (Point **)HECMW_calloc(TABLE_SIZE, sizeof(Point *));
103 sum_table = TABLE_SIZE;
104 CS_verts_tail = (Point **)HECMW_calloc(sum_table, sizeof(Point *));
105 CS_verts_refer = (Point **)HECMW_calloc(sum_table, sizeof(Point *));
106 for (j = 0; j < sum_table; j++) {
107 if ((CS_verts_refer[j] = CS_verts_tail[j] = CS_verts_head[j] =
108 alloc_verts(VERTEX_PACK)) == NULL)
109 HECMW_vis_memory_exit("verts for hexahedra");
110 }
111
112 CS_polys_tail = (Polygon **)HECMW_malloc(sizeof(Polygon *));
113 if ((*CS_polys_tail = *CS_polys_head =
114 alloc_polygons(POLYGON_PACK)) == NULL)
115 HECMW_vis_memory_exit("CS_polys_tail");
116
117 alpha_cube.isosurf = (int **)HECMW_malloc(sizeof(int *));
118 beta_cube.isosurf = (int **)HECMW_malloc(sizeof(int *));
119 cell = (Cell *)HECMW_malloc(sizeof(Cell));
120 }
121 get_data(sff, mesh, data, tmp_int, cell, tn_component);
122 if (sff->surface_style == 2) {
123 alpha_index = make_tile(mesh, cell, sff->iso_value, 0, &alpha_cube,
124 disamb_flag);
125 beta_index = make_tile(mesh, cell, sff->iso_value, 1, &beta_cube,
126 disamb_flag);
127 if (alpha_index && beta_index) {
128 if ((alpha_index + beta_index) == HEX_NODE_INDEX) {
129 if (!merge_vol_iso(0, cell, sff->iso_value, &alpha_cube, 0,
130 &beta_cube, bdflag[i], &sum_verts,
131 CS_verts_tail, CS_verts_refer, CS_verts_head,
132 CS_polys_tail)) {
133 return 0;
134 }
135 }
136 }
137
138 } else if (sff->surface_style == 3) {
139 alpha_index = make_tile(mesh, cell, 0, 0, &alpha_cube, disamb_flag);
140 beta_index = make_tile(mesh, cell, 0, 1, &beta_cube, disamb_flag);
141 if (alpha_index && beta_index) {
142 if ((alpha_index + beta_index) == HEX_NODE_INDEX) {
143 if (!merge_vol_iso(0, cell, 0, &alpha_cube, 0, &beta_cube,
144 bdflag[i], &sum_verts, CS_verts_tail,
145 CS_verts_refer, CS_verts_head,
146 CS_polys_tail)) {
147 return 0;
148 }
149 }
150 }
151 }
152 } else if ((mesh->elem_type[tmp_int] == 341) ||
153 (mesh->elem_type[tmp_int] == 342)) {
154 /* tetrahedra */
155 if (flag_tetra == 0) {
156 flag_tetra = 1;
157 /* initialize */
158 tetra = (Tetra *)HECMW_malloc(sizeof(Tetra));
159 vertex_hash_table = (Hash_vertex *)HECMW_calloc(
160 mesh->n_node * 2, sizeof(Hash_vertex));
161 tetra_point = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
162 head_patch_tetra =
163 (Head_patch_tetra *)HECMW_malloc(sizeof(Head_patch_tetra));
164 if ((vertex_hash_table == NULL) || (tetra == NULL) ||
165 (tetra_point == NULL) || (head_patch_tetra == NULL))
166 HECMW_vis_memory_exit("initialize tetra");
167 for (j = 0; j < mesh->n_node * 2; j++) {
168 vertex_hash_table[j].ident = 0;
169 vertex_hash_table[j].next_vertex = NULL;
170 }
171 tetra_point->ident = 0;
172 head_patch_tetra->num_patch = 0;
173 }
174 get_tetra_data(sff, mesh, data, tmp_int, tetra, tn_component);
175 if (sff->surface_style == 3) sff->iso_value = 0.0;
176 find_intersection_tetra(tetra, sff->iso_value, tetra_point,
177 head_patch_tetra, vertex_hash_table);
178 }
179 }
180 if (flag_tetra == 1) {
181 num_nh_verts = tetra_point->ident;
182 num_nh_patch = head_patch_tetra->num_patch;
183 HECMW_free(tetra);
184 for (j = 0; j < mesh->n_node * 2; j++) {
185 h1 = vertex_hash_table[j].next_vertex;
186 for (i = 0; i < vertex_hash_table[j].ident; i++) {
187 h2 = h1->next_vertex;
188 HECMW_free(h1);
189 h1 = h2;
190 }
191 }
192 HECMW_free(vertex_hash_table);
193 }
194
195 if (flag_hexa > 0) {
196 mfree(alpha_cube.isosurf);
197 mfree(beta_cube.isosurf);
198 mfree(CS_verts_tail);
199 mfree(CS_verts_refer);
200 mfree(CS_polys_tail);
201
202 *sum_v = sum_verts;
203 *sum_t = sum_table;
204 HECMW_free(cell);
205
206 CS_polys = *CS_polys_head;
207
208 sum_polys = 0;
209
210 aplist_size = bplist_size = cplist_size = 1;
211 while (CS_polys->plist != NULL) {
212 switch (CS_polys->type) {
213 case 0:
214 aplist_size += CS_polys->plist[0] + 1;
215 sum_polys++;
216 break;
217 case 1:
218 bplist_size += CS_polys->plist[0] + 1;
219 sum_polys++;
220 break;
221 case 2:
222 cplist_size += (CS_polys->plist[0] - 2) * 4;
223 sum_polys += CS_polys->plist[0] - 2;
224 break;
225 }
226 CS_polys = CS_polys->nextpolygon;
227 }
228 }
229 /*sum_polys=sum_polys*2;*/
230 if ((sum_verts + num_nh_verts) > 0) {
231 result[sf_i].n_vertex = sum_verts + num_nh_verts;
232 result[sf_i].n_patch = sum_polys + num_nh_patch;
233 result[sf_i].vertex =
234 (double *)HECMW_calloc(result[sf_i].n_vertex * 3, sizeof(double));
235 result[sf_i].patch =
236 (int *)HECMW_calloc(result[sf_i].n_patch * 3, sizeof(int));
237 result[sf_i].color =
238 (double *)HECMW_calloc(result[sf_i].n_vertex, sizeof(double));
239 if ((result[sf_i].vertex == NULL) || (result[sf_i].patch == NULL) ||
240 (result[sf_i].color == NULL))
241 HECMW_vis_memory_exit("result");
242 }
243 mm = 0;
244 if (sum_verts > 0) {
245 /* make main vertex table of CS */
246 for (i = 0; i < sum_table; i++) {
247 CS_verts_tmp = CS_verts = CS_verts_head[i];
248 j = 0;
249 while (CS_verts->ident != 0) {
250 /* verts_geom[CS_verts->ident] = CS_verts->geom;
251 verts_field[CS_verts->ident] = CS_verts->field;
252 verts_color[CS_verts->ident] = CS_verts->cdata;
253
254 fprintf(vfile, " %d %lf %lf %lf\n", *tvertex+CS_verts->ident,
255 verts_geom[CS_verts->ident].x,
256 verts_geom[CS_verts->ident].y,
257 verts_geom[CS_verts->ident].z);
258 */
259 result[sf_i].vertex[mm * 3] = CS_verts->geom.x;
260 result[sf_i].vertex[mm * 3 + 1] = CS_verts->geom.y;
261 result[sf_i].vertex[mm * 3 + 2] = CS_verts->geom.z;
262 result[sf_i].color[mm] = CS_verts->cdata;
263 mm++;
264 CS_verts = CS_verts->nextpoint;
265 if (!((++j) % VERTEX_PACK)) {
266 mfree(CS_verts_tmp);
267 CS_verts_tmp = CS_verts;
268 }
269 }
270 mfree(CS_verts_tmp);
271 }
272 mfree(CS_verts_head);
273
274 /* make polygon table and decide vertex ID
275 of each object (alpha,beta,cross) */
276 /* trilist=(Triangle *)HECMW_calloc(sum_polys+1,sizeof(Triangle));
277 */
278 CS_polys_tmp = CS_polys = *CS_polys_head;
279 i = 0;
280 poly_num = 1;
281 minus_patch = 0;
282
283 while (CS_polys->plist != NULL) {
284 switch (CS_polys->type) {
285 case 0:
286 /* fprintf(pfile, "%d ", poly_num+(*tpatch));
287 */
288 k = -1;
289 for (j = 1; j <= CS_polys->plist[0]; j++) {
290 k++;
291 /* fprintf(pfile, "%d ",
292 (*tvertex)+CS_polys->plist[j]);
293
294 trilist[poly_num -1].vertex[k]=CS_polys->plist[j];
295 */
296 if ((sff->output_type == 1) || (sff->output_type == 2))
297 result[sf_i].patch[(poly_num - 1) * 3 + k] =
298 *tvertex + CS_polys->plist[j];
299 else
300 result[sf_i].patch[(poly_num - 1) * 3 + k] = CS_polys->plist[j];
301
302 tmp_data[k * 3] =
303 result[sf_i].vertex[(CS_polys->plist[j] - 1) * 3];
304 tmp_data[k * 3 + 1] =
305 result[sf_i].vertex[(CS_polys->plist[j] - 1) * 3 + 1];
306 tmp_data[k * 3 + 2] =
307 result[sf_i].vertex[(CS_polys->plist[j] - 1) * 3 + 2];
308 }
309 if (((fabs(tmp_data[0] - tmp_data[3]) < EPSILON) &&
310 (fabs(tmp_data[1] - tmp_data[4]) < EPSILON) &&
311 (fabs(tmp_data[2] - tmp_data[5]) < EPSILON)) ||
312 ((fabs(tmp_data[0] - tmp_data[6]) < EPSILON) &&
313 (fabs(tmp_data[1] - tmp_data[7]) < EPSILON) &&
314 (fabs(tmp_data[2] - tmp_data[8]) < EPSILON)) ||
315 ((fabs(tmp_data[6] - tmp_data[3]) < EPSILON) &&
316 (fabs(tmp_data[7] - tmp_data[4]) < EPSILON) &&
317 (fabs(tmp_data[8] - tmp_data[5]) < EPSILON)))
318 minus_patch++;
319 else
320 poly_num++;
321
322 break;
323 }
324
325 CS_polys = CS_polys->nextpolygon;
326 if (!((++i) % POLYGON_PACK)) {
327 mfree(CS_polys_tmp->plist);
328
329 mfree(CS_polys_tmp);
330 CS_polys_tmp = CS_polys;
331 }
332 }
333 /*#ifdef DEBUG
334 fprintf(stderr, "the previous patch num is %d now is %d\n",
335 result[sf_i].n_patch, result[sf_i].n_patch-minus_patch);
336 #endif
337 */
338 result[sf_i].n_patch -= minus_patch;
339
340 mfree(CS_polys_tmp);
341 }
342 if (num_nh_verts > 0) {
343 p1 = tetra_point->nextpoint;
344 for (i = 0; i < tetra_point->ident; i++) {
345 for (j = 0; j < 3; j++)
346 result[sf_i].vertex[(sum_verts + p1->ident) * 3 + j] = p1->geom[j];
347 result[sf_i].color[sum_verts + p1->ident] = p1->cdata;
348 p2 = p1;
349 p1 = p1->nextpoint;
350 HECMW_free(p2);
351 }
352 HECMW_free(tetra_point);
353 }
354 if (num_nh_patch > 0) {
355 t1 = head_patch_tetra->patch_link;
356 for (i = 0; i < head_patch_tetra->num_patch; i++) {
357 for (j = 0; j < 3; j++) {
358 if (sff->output_type == 3)
359 result[sf_i].patch[(sum_polys - minus_patch + i) * 3 + j] =
360 t1->patch[j] + 1 + sum_verts;
361 else
362 result[sf_i].patch[(sum_polys - minus_patch + i) * 3 + j] =
363 *tvertex + t1->patch[j] + 1 + sum_verts;
364 }
365 t2 = t1;
366 t1 = t1->next_patch;
367 HECMW_free(t2);
368 }
369 HECMW_free(head_patch_tetra);
370 }
371
372 /*
373 fprintf(stderr, "On surface %d PE %d: n_vertex is %d n_patch is %d\n", sf_i,
374 mynode, result[sf_i].n_vertex, result[sf_i].n_patch);
375 */
376 if (result[sf_i].n_vertex > 0) {
377 *minc = *maxc = result[sf_i].color[0];
378 for (i = 1; i <= result[sf_i].n_vertex; i++) {
379 if (result[sf_i].color[i - 1] < (*minc))
380 (*minc) = result[sf_i].color[i - 1];
381 if (result[sf_i].color[i - 1] > (*maxc))
382 (*maxc) = result[sf_i].color[i - 1];
383 }
384 /*
385 #ifdef DEBUG
386 fprintf(stderr,"On surface %d PE %d: minimum color=%lf maximum color=%lf\n",
387 sf_i, mynode, *minc,*maxc);
388 #endif
389 */
390 }
391 } /* endof if elem_type>300 */
392 else if ((mesh->elem_type[0] > 200) && (mesh->elem_type[0] < 300)) {
393 result[sf_i].n_patch = 0;
394 for (i = 0; i < n_elem; i++) {
395 if (mesh->elem_type[i] == 231)
396 result[sf_i].n_patch++;
397 else if (mesh->elem_type[i] == 241)
398 result[sf_i].n_patch += 2;
399 }
400 result[sf_i].n_vertex = mesh->n_node;
401 result[sf_i].vertex =
402 (double *)HECMW_calloc(mesh->n_node * 3, sizeof(double));
403 result[sf_i].color = (double *)HECMW_calloc(mesh->n_node, sizeof(double));
404 result[sf_i].patch =
405 (int *)HECMW_calloc(result[sf_i].n_patch * 3, sizeof(int));
406 if ((result[sf_i].vertex == NULL) || (result[sf_i].color == NULL) ||
407 (result[sf_i].patch == NULL))
408 HECMW_vis_memory_exit("result: vertex, color and patch");
409 for (i = 0; i < mesh->n_node; i++) {
410 for (j = 0; j < 3; j++)
411 result[sf_i].vertex[i * 3 + j] = mesh->node[i * 3 + j];
412 }
413 if (data->nn_dof[sff->color_comp] == 1) {
414 c_base = 0;
415 for (i = 0; i < sff->color_comp; i++) c_base += data->nn_dof[i];
416 } else if (data->nn_dof[sff->color_comp] > 1) {
417 c_base = 0;
418 for (i = 0; i < sff->color_comp; i++) c_base += data->nn_dof[i];
419 }
420
421 if ((data->nn_dof[sff->color_comp] > 1) && (sff->color_subcomp == 0)) {
422 for (i = 0; i < mesh->n_node; i++) {
423 result[sf_i].color[i] = 0.0;
424 for (j = 0; j < data->nn_dof[sff->color_comp]; j++) {
425 tmp = data->node_val_item[c_base + i * tn_component + j];
426 result[sf_i].color[i] += tmp * tmp;
427 }
428 result[sf_i].color[i] = sqrt(result[sf_i].color[i]);
429 }
430 }
431
432 else if (data->nn_dof[sff->color_comp] > 1) {
433 for (i = 0; i < mesh->n_node; i++) {
434 result[sf_i].color[i] = data->node_val_item[c_base + i * tn_component +
435 (sff->color_subcomp - 1)];
436 }
437 } else if (data->nn_dof[sff->color_comp] == 1) {
438 for (i = 0; i < mesh->n_node; i++) {
439 result[sf_i].color[i] = data->node_val_item[c_base + i * tn_component];
440 }
441 }
442 poly_num = 0;
443 for (i = 0; i < n_elem; i++) {
444 if (mesh->elem_type[i] == 231) {
445 for (j = 0; j < 3; j++)
446 result->patch[poly_num * 3 + j] =
447 mesh->elem_node_item[mesh->elem_node_index[i] + j] + *tvertex;
448 poly_num++;
449 } else if (mesh->elem_type[i] == 241) {
450 for (j = 0; j < 3; j++)
451 result->patch[poly_num * 3 + j] =
452 mesh->elem_node_item[mesh->elem_node_index[i] + j] + *tvertex;
453 poly_num++;
454 result->patch[poly_num * 3 + 0] =
455 mesh->elem_node_item[mesh->elem_node_index[i] + 0] + *tvertex;
456 result->patch[poly_num * 3 + 1] =
457 mesh->elem_node_item[mesh->elem_node_index[i] + 2] + *tvertex;
458 result->patch[poly_num * 3 + 2] =
459 mesh->elem_node_item[mesh->elem_node_index[i] + 3] + *tvertex;
460 poly_num++;
461 }
462 }
463 /*
464 fprintf(stderr, "n_vertex is %d n_patch is %d\n", result[sf_i].n_vertex,
465 result[sf_i].n_patch);
466 */
467 if (result[sf_i].n_vertex > 0) {
468 *minc = *maxc = result[sf_i].color[0];
469 for (i = 1; i <= result[sf_i].n_vertex; i++) {
470 if (result[sf_i].color[i - 1] < (*minc))
471 (*minc) = result[sf_i].color[i - 1];
472 if (result[sf_i].color[i - 1] > (*maxc))
473 (*maxc) = result[sf_i].color[i - 1];
474 }
475 /*
476 #ifdef DEBUG
477 fprintf(stderr,"On surface %d PE %d: minimum color=%lf maximum color=%lf\n",
478 sf_i, mynode, *minc,*maxc);
479 #endif
480 */
481 }
482 }
483
484 (*tvertex) += result[sf_i].n_vertex;
485 (*tpatch) += result[sf_i].n_patch;
486 /* if(sum_verts>0) {
487 *minv=mincolor;
488 *maxv=maxcolor;
489 }
490
491 mfree(verts_info);
492 mfree(verts_geom);
493 mfree(verts_field);
494 mfree(verts_color);
495 mfree(trilist);
496 */
497
498 return (1);
499 }
500
HECMW_vis_chk_bounds(struct hecmwST_local_mesh * mesh,int * bdflag)501 int HECMW_vis_chk_bounds(struct hecmwST_local_mesh *mesh, int *bdflag) {
502 int i;
503 int n_elem;
504
505 n_elem = mesh->n_elem;
506 for (i = 0; i < n_elem; i++) {
507 *(bdflag + i) = 64;
508 }
509 return 1;
510 }
511
512 #ifdef old
chk_node_data(struct visual_buf * v,int s_comp,int c_comp)513 int chk_node_data(struct visual_buf *v, int s_comp, int c_comp) {
514 int *global_node_id;
515 double *shape_data;
516 double *color_data;
517 int imp_node;
518 int s_base, c_base;
519 int i, j, k;
520 int ne;
521 int n_export_node, export_base;
522 int n_import_node, import_base;
523
524 HECMW_Status stat;
525
526 s_base = c_base = 0;
527
528 for (i = 0; i < s_comp; i++) {
529 s_base += v->node.n_free[i] * v->node.n;
530 }
531 for (i = 0; i < c_comp; i++) {
532 c_base += v->node.n_free[i] * v->node.n;
533 }
534
535 export_base = 0;
536 import_base = 0;
537 for (i = 0; i < v->mesh->n_neighbor_pe; i++) {
538 ne = v->mesh->neighbor_pe[i];
539 n_export_node = v->mesh->export_index[i + 1] - v->mesh->export_index[i];
540 global_node_id = (int *)HECMW_calloc(n_export_node, sizeof(int));
541 shape_data = (double *)HECMW_calloc(n_export_node, sizeof(double));
542 color_data = (double *)HECMW_calloc(n_export_node, sizeof(double));
543 for (j = 0; j < n_export_node; j++) {
544 global_node_id[j] =
545 v->mesh->global_node_id[v->mesh->export_node[export_base + j] - 1];
546 shape_data[j] = (double)*(v->node.data + s_base +
547 v->mesh->export_node[export_base + j] - 1);
548 color_data[j] = (double)*(v->node.data + c_base +
549 v->mesh->export_node[export_base + j] - 1);
550 }
551 HECMW_Send(&n_export_node, 1, HECMW_INT, ne, 0, geofem_app_comm);
552 HECMW_Send(global_node_id, n_export_node, HECMW_INT, ne, 0,
553 geofem_app_comm);
554 HECMW_Send(shape_data, n_export_node, HECMW_DOUBLE, ne, 0, geofem_app_comm);
555 HECMW_Send(color_data, n_export_node, HECMW_DOUBLE, ne, 0, geofem_app_comm);
556
557 HECMW_free(global_node_id);
558 HECMW_free(shape_data);
559 HECMW_free(color_data);
560 export_base = v->mesh->export_index[i + 1];
561
562 n_import_node = v->mesh->import_index[i + 1] - v->mesh->import_index[i];
563 HECMW_Recv(&imp_node, 1, HECMW_INT, ne, HECMW_ANY_TAG, geofem_app_comm,
564 &stat);
565
566 if (imp_node != n_import_node) {
567 GEOFEM_abort(711, "chk_node_data error\n");
568 }
569 global_node_id = (int *)HECMW_calloc(n_import_node, sizeof(int));
570 shape_data = (double *)HECMW_calloc(n_import_node, sizeof(double));
571 color_data = (double *)HECMW_calloc(n_import_node, sizeof(double));
572 HECMW_Recv(global_node_id, n_import_node, HECMW_INT, ne, HECMW_ANY_TAG,
573 geofem_app_comm, &stat);
574 HECMW_Recv(shape_data, n_import_node, HECMW_DOUBLE, ne, HECMW_ANY_TAG,
575 geofem_app_comm, &stat);
576 HECMW_Recv(color_data, n_import_node, HECMW_DOUBLE, ne, HECMW_ANY_TAG,
577 geofem_app_comm, &stat);
578 for (k = import_base; k < import_base + n_import_node; k++) {
579 for (j = 0; j < n_import_node; j++) {
580 if (v->mesh->global_node_id[v->mesh->import_node[k] - 1] ==
581 global_node_id[j]) {
582 v->node.data[s_base + v->mesh->import_node[k] - 1] = shape_data[j];
583 v->node.data[c_base + v->mesh->import_node[k] - 1] = color_data[j];
584 break;
585 }
586 }
587 }
588 HECMW_free(global_node_id);
589 HECMW_free(shape_data);
590 HECMW_free(color_data);
591 import_base = v->mesh->import_index[i + 1];
592 }
593
594 return 1;
595 }
596 #endif
597
598 /*
599 int chk_bounds(struct visual_buf *vol, int *bdflag)
600 {
601 int i, j, k, l;
602 int bd1, bd2;
603 int n_elem;
604
605 n_elem = vol->mesh->n_elem;
606 for (i = 0; i < n_elem; i++) {
607 for (j = 0; j < i; j++) {
608 if (*(bdflag + j) < HEX_FACE_INDEX) {
609 bd1 = bd2 = 0;
610 for (k = 0; k < HEX_N_NODE; k++) {
611 for (l = 0; l < HEX_N_NODE; l++) {
612 if ((*(vol->mesh->elem + i + n_elem * k))
613 == (*(vol->mesh->elem + j + n_elem * l))) {
614 bd1 += 1 << k;
615 bd2 += 1 << l;
616 }
617 }
618 }
619 switch (bd1) {
620 case 15: *(bdflag + i) += 1; break;
621 case 102: *(bdflag + i) += 2; break;
622 case 240: *(bdflag + i) += 4; break;
623 case 153: *(bdflag + i) += 8; break;
624 case 204: *(bdflag + i) += 16; break;
625 case 51: *(bdflag + i) += 32; break;
626 default: break;
627 }
628 switch (bd2) {
629 case 15: *(bdflag + j) += 1; break;
630 case 102: *(bdflag + j) += 2; break;
631 case 240: *(bdflag + j) += 4; break;
632 case 153: *(bdflag + j) += 8; break;
633 case 204: *(bdflag + j) += 16; break;
634 case 51: *(bdflag + j) += 32; break;
635 default: break;
636 }
637 }
638
639 if (*(bdflag + i) == HEX_FACE_INDEX) break;
640 }
641 }
642
643 return 1;
644
645 }
646
647 */
find_isoline(int isonum,int sum_polys,double mincolor,double maxcolor,double * vcoord,int * plist,double * vcolor,Isohead * isohead)648 void find_isoline(int isonum, int sum_polys, double mincolor, double maxcolor,
649 double *vcoord, int *plist, double *vcolor,
650 Isohead *isohead) {
651 int i, j, k;
652 double c[3];
653 double isocolor, deltac;
654 Fgeom g[3];
655
656 if (isonum <= 0) {
657 fprintf(stderr, "isonumber input wrong\n");
658 exit(0);
659 }
660 deltac = (maxcolor - mincolor) / (isonum + 1);
661 for (k = 0; k < isonum; k++) {
662 isocolor = mincolor + (k + 1) * deltac;
663 for (i = 0; i < sum_polys; i++) {
664 for (j = 0; j < 3; j++) {
665 c[j] = vcolor[plist[i * 3 + j] - 1];
666 g[j].x = vcoord[(plist[i * 3 + j] - 1) * 3];
667 g[j].y = vcoord[(plist[i * 3 + j] - 1) * 3 + 1];
668 g[j].z = vcoord[(plist[i * 3 + j] - 1) * 3 + 2];
669 }
670 line_find(isocolor, c, g, k, isohead);
671 }
672 }
673 }
674
find_line_segment(double f[3][3],double c[3],double isocolor,double iso_p[6])675 int find_line_segment(double f[3][3], double c[3], double isocolor,
676 double iso_p[6]) {
677 int j, m, j1;
678 int pnum, flag;
679
680 flag = 0;
681 pnum = -1;
682 for (j = 0; j < 3; j++) {
683 j1 = j + 1;
684 if (j1 == 3) j1 = 0;
685 if ((fabs(c[j] - c[j1]) < 0.0000001) &&
686 (fabs(c[j] - isocolor) < 0.0000001)) {
687 /*this edge is isoline*/
688 flag = 1;
689
690 for (m = 0; m < 3; m++) iso_p[m] = f[j][m];
691 for (m = 0; m < 3; m++) iso_p[m + 3] = f[j1][m];
692 return 1;
693 } else if (((c[j] >= isocolor) && (c[j1] < isocolor)) ||
694 ((c[j] < isocolor) && (c[j1] >= isocolor))) {
695 pnum++;
696 for (m = 0; m < 3; m++)
697 iso_p[pnum * 3 + m] =
698 f[j][m] + (isocolor - c[j]) / (c[j1] - c[j]) * (f[j1][m] - f[j][m]);
699 }
700 }
701 if (pnum == 1) flag = 1;
702
703 return flag;
704 }
705
line_find(double isocolor,double c[3],Fgeom g[3],int k,Isohead * isohead)706 void line_find(double isocolor, double c[3], Fgeom g[3], int k,
707 Isohead *isohead) {
708 Isoline *p1, *p2;
709 int j, m, j1;
710 int pnum;
711 Fgeom isopoint[2];
712
713 pnum = -1;
714 for (j = 0; j < 2; j++) {
715 isopoint[j].x = isopoint[j].y = isopoint[j].z = 0.0;
716 }
717 for (j = 0; j < 3; j++) {
718 j1 = j + 1;
719 if (j1 == 3) j1 = 0;
720 if ((fabs(c[j] - c[j1]) < 0.0000001) &&
721 (fabs(c[j] - isocolor) < 0.0000001)) {
722 /*this edge is isoline*/
723 isohead[k].linenum++;
724 p1 = isohead[k].nextline;
725 p2 = (Isoline *)HECMW_malloc(sizeof(Isoline));
726
727 p2->point[0].x = g[j].x;
728 p2->point[0].y = g[j].y;
729 p2->point[0].z = g[j].z;
730 p2->point[1].x = g[j1].x;
731 p2->point[1].y = g[j1].y;
732 p2->point[1].z = g[j1].z;
733
734 p2->nextline = p1;
735 isohead[k].nextline = p2;
736 return;
737 }
738
739 else if (((c[j] >= isocolor) && (c[j1] < isocolor)) ||
740 ((c[j] < isocolor) && (c[j1] >= isocolor))) {
741 pnum++;
742 isopoint[pnum].x =
743 g[j].x + (isocolor - c[j]) / (c[j1] - c[j]) * (g[j1].x - g[j].x);
744 isopoint[pnum].y =
745 g[j].y + (isocolor - c[j]) / (c[j1] - c[j]) * (g[j1].y - g[j].y);
746 isopoint[pnum].z =
747 g[j].z + (isocolor - c[j]) / (c[j1] - c[j]) * (g[j1].z - g[j].z);
748 }
749 }
750 if (pnum == 1) {
751 isohead[k].linenum++;
752 p1 = isohead[k].nextline;
753 p2 = (Isoline *)HECMW_malloc(sizeof(Isoline));
754 for (m = 0; m < 2; m++) {
755 p2->point[m].x = isopoint[m].x;
756 p2->point[m].y = isopoint[m].y;
757 p2->point[m].z = isopoint[m].z;
758 }
759 p2->nextline = p1;
760 isohead[k].nextline = p2;
761 }
762
763 return;
764 }
765
isoline_free(int isonum,Isohead * isohead)766 void isoline_free(int isonum, Isohead *isohead) {
767 Isoline *p1, *p2;
768 int i;
769 for (i = 0; i < isonum; i++) {
770 p1 = isohead[i].nextline;
771 while (p1 != NULL) {
772 p2 = p1;
773 p1 = p1->nextline;
774 HECMW_free(p2);
775 }
776 }
777 HECMW_free(isohead);
778 return;
779 }
780