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_connectivity_build.h"
7 
8 #include <math.h>
9 #include "hecmw_vis_mem_util.h"
10 #include "hecmw_vis_comm_util.h"
11 #include "hecmw_malloc.h"
12 
find_index_connectivity(struct hecmwST_local_mesh * mesh,int * index_connect)13 void find_index_connectivity(struct hecmwST_local_mesh *mesh,
14                              int *index_connect) {
15   int i;
16 
17   index_connect[0] = 0;
18   for (i = 0; i < mesh->n_elem; i++) {
19     if ((mesh->elem_type[i] == 341) || (mesh->elem_type[i] == 342) ||
20         (mesh->elem_type[i] == 3414))
21       index_connect[i + 1] = index_connect[i] + 4;
22     else if ((mesh->elem_type[i] == 351) || (mesh->elem_type[i] == 352))
23       index_connect[i + 1] = index_connect[i] + 5;
24     else if ((mesh->elem_type[i] == 361) || (mesh->elem_type[i] == 362))
25       index_connect[i + 1] = index_connect[i] + 6;
26     else if (mesh->elem_type[i] > 400)
27       index_connect[i + 1] = index_connect[i] + 0;
28   }
29   return;
30 }
31 
find_index_a_connect(struct hecmwST_local_mesh * mesh,int num_export,int pe_no,int * export_element,int * index_a_connect)32 void find_index_a_connect(struct hecmwST_local_mesh *mesh, int num_export,
33                           int pe_no, int *export_element,
34                           int *index_a_connect) {
35   int i;
36 
37   index_a_connect[0] = 0;
38   for (i = 0; i < num_export; i++) {
39     if ((mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] == 341) ||
40         (mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] == 342) ||
41         (mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] == 3414))
42       index_a_connect[i + 1] = index_a_connect[i] + 4;
43     else if ((mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] ==
44               351) ||
45              (mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] == 352))
46       index_a_connect[i + 1] = index_a_connect[i] + 5;
47     else if ((mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] ==
48               361) ||
49              (mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] == 362))
50       index_a_connect[i + 1] = index_a_connect[i] + 6;
51     else if (mesh->elem_type[export_element[pe_no * mesh->n_elem + i]] > 400)
52       index_a_connect[i + 1] = index_a_connect[i] + 0;
53   }
54   return;
55 }
56 
generate_face(int index_face_tetra[5],int face_tetra[12],int index_face_prism[6],int face_prism[18],int index_face_hexa[7],int face_hexa[24])57 void generate_face(int index_face_tetra[5], int face_tetra[12],
58                    int index_face_prism[6], int face_prism[18],
59                    int index_face_hexa[7], int face_hexa[24]) {
60   int i;
61   index_face_tetra[0] = 0;
62   for (i = 0; i < 4; i++) index_face_tetra[i + 1] = (i + 1) * 3;
63   face_tetra[0]                                   = 0;
64   face_tetra[1]                                   = 2;
65   face_tetra[2]                                   = 1;
66   face_tetra[3]                                   = 0;
67   face_tetra[4]                                   = 1;
68   face_tetra[5]                                   = 3;
69   face_tetra[6]                                   = 1;
70   face_tetra[7]                                   = 2;
71   face_tetra[8]                                   = 3;
72   face_tetra[9]                                   = 2;
73   face_tetra[10]                                  = 0;
74   face_tetra[11]                                  = 3;
75 
76   index_face_prism[0] = 0;
77   for (i = 0; i < 3; i++) index_face_prism[i + 1] = (i + 1) * 4;
78   index_face_prism[4]                             = 15;
79   index_face_prism[5]                             = 18;
80   face_prism[0]                                   = 0;
81   face_prism[1]                                   = 1;
82   face_prism[2]                                   = 4;
83   face_prism[3]                                   = 3;
84   face_prism[4]                                   = 1;
85   face_prism[5]                                   = 2;
86   face_prism[6]                                   = 5;
87   face_prism[7]                                   = 4;
88   face_prism[8]                                   = 2;
89   face_prism[9]                                   = 0;
90   face_prism[10]                                  = 3;
91   face_prism[11]                                  = 5;
92   face_prism[12]                                  = 0;
93   face_prism[13]                                  = 2;
94   face_prism[14]                                  = 1;
95   face_prism[15]                                  = 3;
96   face_prism[16]                                  = 4;
97   face_prism[17]                                  = 5;
98 
99   index_face_hexa[0] = 0;
100   for (i = 0; i < 6; i++) index_face_hexa[i + 1] = (i + 1) * 4;
101   face_hexa[0]                                   = 0;
102   face_hexa[1]                                   = 4;
103   face_hexa[2]                                   = 7;
104   face_hexa[3]                                   = 3;
105   face_hexa[4]                                   = 1;
106   face_hexa[5]                                   = 2;
107   face_hexa[6]                                   = 6;
108   face_hexa[7]                                   = 5;
109   face_hexa[8]                                   = 0;
110   face_hexa[9]                                   = 1;
111   face_hexa[10]                                  = 5;
112   face_hexa[11]                                  = 4;
113   face_hexa[12]                                  = 2;
114   face_hexa[13]                                  = 3;
115   face_hexa[14]                                  = 7;
116   face_hexa[15]                                  = 6;
117   face_hexa[16]                                  = 3;
118   face_hexa[17]                                  = 2;
119   face_hexa[18]                                  = 1;
120   face_hexa[19]                                  = 0;
121   face_hexa[20]                                  = 4;
122   face_hexa[21]                                  = 5;
123   face_hexa[22]                                  = 6;
124   face_hexa[23]                                  = 7;
125   return;
126 }
127 
add_to_hash(int elemID,int faceID,int hashID,Hash_table * h_table)128 void add_to_hash(int elemID, int faceID, int hashID, Hash_table *h_table) {
129   Hash_table *p1, *p2;
130 
131   if ((p1 = (Hash_table *)HECMW_malloc(sizeof(Hash_table))) == NULL)
132     HECMW_vis_memory_exit("hash_table: p1");
133   h_table[hashID].elemID++;
134   p2                        = h_table[hashID].next_elem;
135   h_table[hashID].next_elem = p1;
136   p1->elemID                = elemID;
137   p1->faceID                = faceID;
138   p1->next_elem             = p2;
139   return;
140 }
141 
build_hash_table(struct hecmwST_local_mesh * mesh,int * index_connect,Hash_table * h_table,int index_face_tetra[5],int face_tetra[12],int index_face_prism[6],int face_prism[18],int index_face_hexa[7],int face_hexa[24])142 void build_hash_table(struct hecmwST_local_mesh *mesh, int *index_connect,
143                       Hash_table *h_table, int index_face_tetra[5],
144                       int face_tetra[12], int index_face_prism[6],
145                       int face_prism[18], int index_face_hexa[7],
146                       int face_hexa[24]) {
147   int i, j, k, n_elem, node[MAX_N_NODE], nodesum;
148 
149   n_elem = mesh->n_elem;
150   for (i = 0; i < n_elem; i++) {
151     for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
152          j++)
153       node[j] = mesh->elem_node_item[mesh->elem_node_index[i] + j];
154     if ((mesh->elem_type[i] == 341) || (mesh->elem_type[i] == 342) ||
155         (mesh->elem_type[i] == 3414)) {
156       for (j = 0; j < 4; j++) {
157         nodesum = 0;
158         for (k = index_face_tetra[j]; k < index_face_tetra[j + 1]; k++)
159           nodesum += node[face_tetra[k]];
160         add_to_hash(i, j, nodesum, h_table);
161       }
162     } else if ((mesh->elem_type[i] == 351) || (mesh->elem_type[i] == 352)) {
163       for (j = 0; j < 5; j++) {
164         nodesum = 0;
165         for (k = index_face_prism[j]; k < index_face_prism[j + 1]; k++)
166           nodesum += node[face_prism[k]];
167         add_to_hash(i, j, nodesum, h_table);
168       }
169     } else if ((mesh->elem_type[i] == 361) || (mesh->elem_type[i] == 362)) {
170       for (j = 0; j < 6; j++) {
171         nodesum = 0;
172         for (k = index_face_hexa[j]; k < index_face_hexa[j + 1]; k++)
173           nodesum += node[face_hexa[k]];
174         add_to_hash(i, j, nodesum, h_table);
175       }
176     }
177   }
178   return;
179 }
180 
is_equal_array(int n[4],int nn[4],int num)181 int is_equal_array(int n[4], int nn[4], int num) {
182   int tmp;
183   int i, j;
184   int flag = 1;
185   for (i = 0; i < num - 1; i++)
186     for (j = i + 1; j < num; j++) {
187       if (n[i] < n[j]) {
188         tmp  = n[i];
189         n[i] = n[j];
190         n[j] = tmp;
191       }
192       if (nn[i] < nn[j]) {
193         tmp   = nn[i];
194         nn[i] = nn[j];
195         nn[j] = tmp;
196       }
197     }
198   for (i = 0; i < num; i++) {
199     if (n[i] != nn[i]) flag = 0;
200   }
201   return (flag);
202 }
203 
is_connect(int elemID1,int faceID1,int elemID2,int faceID2,struct hecmwST_local_mesh * mesh,int index_face_tetra[5],int face_tetra[12],int index_face_prism[6],int face_prism[18],int index_face_hexa[7],int face_hexa[24])204 int is_connect(int elemID1, int faceID1, int elemID2, int faceID2,
205                struct hecmwST_local_mesh *mesh, int index_face_tetra[5],
206                int face_tetra[12], int index_face_prism[6], int face_prism[18],
207                int index_face_hexa[7], int face_hexa[24]) {
208   int j, node1[MAX_N_NODE], node2[MAX_N_NODE], face1[4], face2[4], k, flag,
209       f_num1, f_num2;
210   flag = 0;
211   for (j = 0;
212        j < mesh->elem_node_index[elemID1 + 1] - mesh->elem_node_index[elemID1];
213        j++)
214     node1[j] = mesh->elem_node_item[mesh->elem_node_index[elemID1] + j];
215   for (j = 0;
216        j < mesh->elem_node_index[elemID2 + 1] - mesh->elem_node_index[elemID2];
217        j++)
218     node2[j] = mesh->elem_node_item[mesh->elem_node_index[elemID2] + j];
219 
220   if ((mesh->elem_type[elemID1] == 341) || (mesh->elem_type[elemID1] == 342) ||
221       (mesh->elem_type[elemID1] == 3414)) {
222     f_num1 = index_face_tetra[faceID1 + 1] - index_face_tetra[faceID1];
223     for (k     = 0; k < f_num1; k++)
224       face1[k] = node1[face_tetra[index_face_tetra[faceID1] + k]];
225   } else if ((mesh->elem_type[elemID1] == 351) ||
226              (mesh->elem_type[elemID1] == 352)) {
227     f_num1 = index_face_prism[faceID1 + 1] - index_face_prism[faceID1];
228     for (k     = 0; k < f_num1; k++)
229       face1[k] = node1[face_prism[index_face_prism[faceID1] + k]];
230   } else if ((mesh->elem_type[elemID1] == 361) ||
231              (mesh->elem_type[elemID1] == 362)) {
232     f_num1 = index_face_hexa[faceID1 + 1] - index_face_hexa[faceID1];
233     for (k     = 0; k < f_num1; k++)
234       face1[k] = node1[face_hexa[index_face_hexa[faceID1] + k]];
235   }
236   if ((mesh->elem_type[elemID2] == 341) || (mesh->elem_type[elemID2] == 342) ||
237       (mesh->elem_type[elemID2] == 3414)) {
238     f_num2 = index_face_tetra[faceID2 + 1] - index_face_tetra[faceID2];
239     for (k     = 0; k < f_num2; k++)
240       face2[k] = node2[face_tetra[index_face_tetra[faceID2] + k]];
241   } else if ((mesh->elem_type[elemID2] == 351) ||
242              (mesh->elem_type[elemID2] == 352)) {
243     f_num2 = index_face_prism[faceID2 + 1] - index_face_prism[faceID2];
244     for (k     = 0; k < f_num2; k++)
245       face2[k] = node2[face_prism[index_face_prism[faceID2] + k]];
246   } else if ((mesh->elem_type[elemID2] == 361) ||
247              (mesh->elem_type[elemID2] == 541) ||
248              (mesh->elem_type[elemID2] == 362) ||
249              (mesh->elem_type[elemID2] == 542)) {
250     f_num2 = index_face_hexa[faceID2 + 1] - index_face_hexa[faceID2];
251     for (k     = 0; k < f_num2; k++)
252       face2[k] = node2[face_hexa[index_face_hexa[faceID2] + k]];
253   }
254   flag                       = 1;
255   if (f_num1 != f_num2) flag = 0;
256   if (flag == 1) {
257     if (!is_equal_array(face1, face2, f_num1)) {
258       flag = 0;
259     }
260   }
261   return (flag);
262 }
263 
find_to_hash(int elemID1,int faceID1,int hashID,Hash_table * h_table,struct hecmwST_local_mesh * mesh,int tmp_connect[2],int index_face_tetra[5],int face_tetra[12],int index_face_prism[6],int face_prism[18],int index_face_hexa[7],int face_hexa[24])264 int find_to_hash(int elemID1, int faceID1, int hashID, Hash_table *h_table,
265                  struct hecmwST_local_mesh *mesh, int tmp_connect[2],
266                  int index_face_tetra[5], int face_tetra[12],
267                  int index_face_prism[6], int face_prism[18],
268                  int index_face_hexa[7], int face_hexa[24]) {
269   int elemID2, faceID2;
270   Hash_table *p1;
271   int flag;
272 
273   flag = 0;
274   if (h_table[hashID].elemID > 0) {
275     p1 = h_table[hashID].next_elem;
276     while (p1 != NULL) {
277       if (p1->elemID != elemID1) {
278         elemID2 = p1->elemID;
279         faceID2 = p1->faceID;
280         if (is_connect(elemID1, faceID1, elemID2, faceID2, mesh,
281                        index_face_tetra, face_tetra, index_face_prism,
282                        face_prism, index_face_hexa, face_hexa)) {
283           tmp_connect[0] = elemID2;
284           tmp_connect[1] = faceID2;
285           flag           = 1;
286           return 1;
287         }
288       }
289       p1 = p1->next_elem;
290     }
291   }
292 
293   return (flag);
294 }
295 
h_free(Hash_table * h_table,int maxadd)296 void h_free(Hash_table *h_table, int maxadd) {
297   int i;
298   Hash_table *p1, *p2;
299   for (i = 0; i < maxadd; i++) {
300     p1 = h_table[i].next_elem;
301     while (p1 != NULL) {
302       p2 = p1;
303       p1 = p1->next_elem;
304       HECMW_free(p2);
305     }
306   }
307   HECMW_free(h_table);
308   return;
309 }
310 
free_b_patch(Boundary_patch * b_patch)311 void free_b_patch(Boundary_patch *b_patch) {
312   Boundary_patch *p1, *p2;
313   p1 = b_patch->next_patch;
314   while (p1 != NULL) {
315     p2 = p1->next_patch;
316     HECMW_free(p1);
317     p1 = p2;
318   }
319   HECMW_free(b_patch);
320   return;
321 }
322 
build_connectivity(struct hecmwST_local_mesh * mesh,Hash_table * h_table,int * index_connect,int * connect,int index_face_tetra[5],int face_tetra[12],int index_face_prism[6],int face_prism[18],int index_face_hexa[7],int face_hexa[24])323 void build_connectivity(struct hecmwST_local_mesh *mesh, Hash_table *h_table,
324                         int *index_connect, int *connect,
325                         int index_face_tetra[5], int face_tetra[12],
326                         int index_face_prism[6], int face_prism[18],
327                         int index_face_hexa[7], int face_hexa[24]) {
328   int i, j, k, n_elem, node[MAX_N_NODE], nodesum, tmp_connect[2];
329   int flag;
330 
331   n_elem = mesh->n_elem;
332   for (i = 0; i < n_elem; i++) {
333     for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
334          j++)
335       node[j] = mesh->elem_node_item[mesh->elem_node_index[i] + j];
336     if ((mesh->elem_type[i] == 341) || (mesh->elem_type[i] == 342) ||
337         (mesh->elem_type[i] == 3414)) {
338       for (j = 0; j < 4; j++) {
339         nodesum = 0;
340         for (k = index_face_tetra[j]; k < index_face_tetra[j + 1]; k++) {
341           nodesum += node[face_tetra[k]];
342         }
343         flag = find_to_hash(i, j, nodesum, h_table, mesh, tmp_connect,
344                             index_face_tetra, face_tetra, index_face_prism,
345                             face_prism, index_face_hexa, face_hexa);
346         if (flag == 1) {
347           connect[(index_connect[i] + j) * 3]     = tmp_connect[0];
348           connect[(index_connect[i] + j) * 3 + 1] = tmp_connect[1];
349         } else {
350           connect[(index_connect[i] + j) * 3]     = -1;
351           connect[(index_connect[i] + j) * 3 + 1] = -1;
352         }
353       }
354     } else if ((mesh->elem_type[i] == 351) || (mesh->elem_type[i] == 352)) {
355       for (j = 0; j < 5; j++) {
356         nodesum = 0;
357         for (k = index_face_prism[j]; k < index_face_prism[j + 1]; k++) {
358           nodesum += node[face_prism[k]];
359         }
360         flag = find_to_hash(i, j, nodesum, h_table, mesh, tmp_connect,
361                             index_face_tetra, face_tetra, index_face_prism,
362                             face_prism, index_face_hexa, face_hexa);
363         if (flag == 1) {
364           connect[(index_connect[i] + j) * 3]     = tmp_connect[0];
365           connect[(index_connect[i] + j) * 3 + 1] = tmp_connect[1];
366         } else {
367           connect[(index_connect[i] + j) * 3]     = -1;
368           connect[(index_connect[i] + j) * 3 + 1] = -1;
369         }
370       }
371     } else if ((mesh->elem_type[i] == 361) || (mesh->elem_type[i] == 541) ||
372                (mesh->elem_type[i] == 362) || (mesh->elem_type[i] == 542)) {
373       for (j = 0; j < 6; j++) {
374         nodesum = 0;
375         for (k = index_face_hexa[j]; k < index_face_hexa[j + 1]; k++) {
376           nodesum += node[face_hexa[k]];
377         }
378         flag = find_to_hash(i, j, nodesum, h_table, mesh, tmp_connect,
379                             index_face_tetra, face_tetra, index_face_prism,
380                             face_prism, index_face_hexa, face_hexa);
381         if (flag == 1) {
382           connect[(index_connect[i] + j) * 3]     = tmp_connect[0];
383           connect[(index_connect[i] + j) * 3 + 1] = tmp_connect[1];
384         } else {
385           connect[(index_connect[i] + j) * 3]     = -1;
386           connect[(index_connect[i] + j) * 3 + 1] = -1;
387         }
388       }
389     }
390   }
391 
392   return;
393 }
394 
add_one_patch(Surface * sff,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int * node_hit,Tetra_point * b_point,Head_patch_tetra * head_b_patch,int node[4],int c_base,int d_base,int tn_component)395 void add_one_patch(Surface *sff, struct hecmwST_local_mesh *mesh,
396                    struct hecmwST_result_data *data, int *node_hit,
397                    Tetra_point *b_point, Head_patch_tetra *head_b_patch,
398                    int node[4], int c_base, int d_base, int tn_component) {
399   int k, m;
400   int index_patch, patch[3];
401   int flag;
402   double c_data;
403   Tetra_point *p1, *p2;
404   Patch_tetra *t1, *t2;
405   double tmp;
406 
407   index_patch = -1;
408   flag        = 1;
409   if ((node[0] == node[1]) || (node[0] == node[2]) || (node[1] == node[2]))
410     flag = 0;
411   if (flag == 1) {
412     for (k = 0; k < 3; k++) {
413       if (node_hit[node[k] - 1] == -1) {
414         p1 = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
415         if (p1 == NULL) HECMW_vis_memory_exit("p1");
416         b_point->ident++;
417         p2                 = b_point->nextpoint;
418         b_point->nextpoint = p1;
419         if ((data->nn_dof[sff->color_comp] > 1) && (sff->color_subcomp == 0)) {
420           c_data = 0.0;
421           for (m = 0; m < data->nn_dof[sff->color_comp]; m++) {
422             tmp =
423                 data->node_val_item[(node[k] - 1) * tn_component + c_base + m];
424             c_data += tmp * tmp;
425           }
426           c_data = sqrt(c_data);
427         } else if ((data->nn_dof[sff->color_comp] > 1) &&
428                    (sff->color_subcomp != 0))
429           c_data = data->node_val_item[(node[k] - 1) * tn_component + c_base +
430                                        (sff->color_subcomp - 1)];
431         else if (data->nn_dof[sff->color_comp] == 1)
432           c_data = data->node_val_item[(node[k] - 1) * tn_component + c_base];
433 
434         p1->cdata = c_data;
435         if (sff->deform_display_on == 1) {
436           for (m = 0; m < 3; m++)
437             p1->disp[m] =
438                 data->node_val_item[(node[k] - 1) * tn_component + d_base + m];
439         }
440         p1->ident     = b_point->ident - 1;
441         p1->geom[0]   = mesh->node[(node[k] - 1) * 3];
442         p1->geom[1]   = mesh->node[(node[k] - 1) * 3 + 1];
443         p1->geom[2]   = mesh->node[(node[k] - 1) * 3 + 2];
444         p1->nextpoint = p2;
445         index_patch++;
446         patch[index_patch]    = b_point->ident;
447         node_hit[node[k] - 1] = b_point->ident;
448       } else {
449         index_patch++;
450         patch[index_patch] = node_hit[node[k] - 1];
451       }
452     }
453     t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
454     if (t1 == NULL) HECMW_vis_memory_exit("t1");
455     t2 = head_b_patch->patch_link;
456     head_b_patch->num_patch++;
457     head_b_patch->patch_link = t1;
458     t1->next_patch           = t2;
459     t1->patch[0]             = patch[0];
460     t1->patch[1]             = patch[1];
461     t1->patch[2]             = patch[2];
462   }
463   return;
464 }
465 
add_two_patch(Surface * sff,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int * node_hit,Tetra_point * b_point,Head_patch_tetra * head_b_patch,int node[4],int c_base,int d_base,int tn_component)466 void add_two_patch(Surface *sff, struct hecmwST_local_mesh *mesh,
467                    struct hecmwST_result_data *data, int *node_hit,
468                    Tetra_point *b_point, Head_patch_tetra *head_b_patch,
469                    int node[4], int c_base, int d_base, int tn_component) {
470   int k, m;
471   int index_patch, patch[4];
472   int flag;
473   double c_data;
474   Tetra_point *p1, *p2;
475   Patch_tetra *t1, *t2;
476   double tmp;
477 
478   index_patch = -1;
479   for (k = 0; k < 4; k++) {
480     if (node_hit[node[k] - 1] == -1) {
481       p1 = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
482       if (p1 == NULL) HECMW_vis_memory_exit("p1");
483       b_point->ident++;
484       p2                 = b_point->nextpoint;
485       b_point->nextpoint = p1;
486       if ((data->nn_dof[sff->color_comp] > 1) && (sff->color_subcomp == 0)) {
487         c_data = 0.0;
488         for (m = 0; m < data->nn_dof[sff->color_comp]; m++) {
489           tmp = data->node_val_item[(node[k] - 1) * tn_component + c_base + m];
490           c_data += tmp * tmp;
491         }
492         c_data = sqrt(c_data);
493       } else if ((data->nn_dof[sff->color_comp] > 1) &&
494                  (sff->color_subcomp != 0))
495         c_data = data->node_val_item[(node[k] - 1) * tn_component + c_base +
496                                      (sff->color_subcomp - 1)];
497       else if (data->nn_dof[sff->color_comp] == 1)
498         c_data = data->node_val_item[(node[k] - 1) * tn_component + c_base];
499 
500       p1->cdata = c_data;
501       if (sff->deform_display_on == 1) {
502         for (m = 0; m < 3; m++)
503           p1->disp[m] =
504               data->node_val_item[(node[k] - 1) * tn_component + d_base + m];
505       }
506 
507       p1->ident     = b_point->ident - 1;
508       p1->geom[0]   = mesh->node[(node[k] - 1) * 3];
509       p1->geom[1]   = mesh->node[(node[k] - 1) * 3 + 1];
510       p1->geom[2]   = mesh->node[(node[k] - 1) * 3 + 2];
511       p1->nextpoint = p2;
512       index_patch++;
513       patch[index_patch]    = b_point->ident;
514       node_hit[node[k] - 1] = b_point->ident;
515     } else {
516       index_patch++;
517       patch[index_patch] = node_hit[node[k] - 1];
518     }
519   }
520   flag = 1;
521   if ((node[0] == node[1]) || (node[0] == node[2]) || (node[1] == node[2]))
522     flag = 0;
523   if (flag == 1) {
524     t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
525     if (t1 == NULL) HECMW_vis_memory_exit("t1");
526     t2 = head_b_patch->patch_link;
527     head_b_patch->num_patch++;
528     head_b_patch->patch_link = t1;
529     t1->next_patch           = t2;
530     t1->patch[0]             = patch[0];
531     t1->patch[1]             = patch[1];
532     t1->patch[2]             = patch[2];
533   }
534   flag = 1;
535   if ((node[0] == node[2]) || (node[0] == node[3]) || (node[2] == node[3]))
536     flag = 0;
537   if (flag == 1) {
538     t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
539     if (t1 == NULL) HECMW_vis_memory_exit("t1");
540     t2 = head_b_patch->patch_link;
541     head_b_patch->num_patch++;
542     head_b_patch->patch_link = t1;
543     t1->next_patch           = t2;
544     t1->patch[0]             = patch[0];
545     t1->patch[1]             = patch[2];
546     t1->patch[2]             = patch[3];
547   }
548 
549   return;
550 }
551 
HECMW_vis_find_boundary_surface(Surface * sff,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int * tvertex,int * tpatch,double * minc,double * maxc,Result * result,int sf_i,HECMW_Comm VIS_COMM,int init_flag,Connect_inf * global_connect)552 void HECMW_vis_find_boundary_surface(Surface *sff,
553                                      struct hecmwST_local_mesh *mesh,
554                                      struct hecmwST_result_data *data,
555                                      int *tvertex, int *tpatch, double *minc,
556                                      double *maxc, Result *result, int sf_i,
557                                      HECMW_Comm VIS_COMM, int init_flag,
558                                      Connect_inf *global_connect)
559 
560 {
561   HECMW_Status stat;
562   int pesize, mynode;
563   Hash_table *h_table;
564   int i, j, k, m;
565   int *export_no_element, *import_no_element, *export_elem, *flag,
566       *index_import_element;
567   int *global_id, *among_connect, **g_id, **a_connect, **index_a;
568   int flag1, nnode[MAX_N_NODE], pe_no, startnode;
569   /*  Boundary_patch *b_patch, *p1, *p2;
570    */
571   int node[4];
572   int c_base;
573   int index_face_tetra[5], face_tetra[12], index_face_prism[6], face_prism[18],
574       index_face_hexa[7], face_hexa[24];
575   int *index_a_connect;
576   int n_elem;
577   int *node_hit;
578   Tetra_point *b_point;
579   Head_patch_tetra *head_b_patch;
580   Tetra_point *p1, *p2;
581   Patch_tetra *t1, *t2;
582   int *nelem_dist;
583   int tmp_sum, tmp_nelem;
584   int *connect, *index_connect;
585   int tn_component;
586   int d_base;
587 
588   HECMW_Comm_size(VIS_COMM, &pesize);
589   HECMW_Comm_rank(VIS_COMM, &mynode);
590   n_elem       = mesh->n_elem;
591   tn_component = 0;
592   for (i = 0; i < data->nn_component; i++) tn_component += data->nn_dof[i];
593   generate_face(index_face_tetra, face_tetra, index_face_prism, face_prism,
594                 index_face_hexa, face_hexa);
595   if (init_flag == 1) {
596     nelem_dist = (int *)HECMW_calloc(pesize + 1, sizeof(int));
597     if (mynode == 0) {
598       nelem_dist[0] = 0;
599       nelem_dist[1] = mesh->ne_internal;
600       tmp_sum       = mesh->ne_internal;
601       for (i = 1; i < pesize; i++) {
602         HECMW_Recv(&tmp_nelem, 1, HECMW_INT, i, HECMW_ANY_TAG, VIS_COMM, &stat);
603         tmp_sum += tmp_nelem;
604         nelem_dist[i + 1] = tmp_sum;
605       }
606       for (i = 1; i < pesize; i++)
607         HECMW_Send(nelem_dist, pesize + 1, HECMW_INT, i, 0, VIS_COMM);
608     } else {
609       HECMW_Send(&mesh->ne_internal, 1, HECMW_INT, 0, 0, VIS_COMM);
610       HECMW_Recv(nelem_dist, pesize + 1, HECMW_INT, 0, HECMW_ANY_TAG, VIS_COMM,
611                  &stat);
612     }
613     /*	if(mynode==0) {
614     for(i=0;i<pesize+1;i++)
615             fprintf(stderr, "nelem_dist=%d ", nelem_dist[i]);
616     fprintf(stderr, "\n");
617 }
618      */
619     mesh->global_elem_ID = (int *)HECMW_calloc(mesh->n_elem, sizeof(int));
620     if (mesh->global_elem_ID == NULL) HECMW_vis_memory_exit("Global_elem_ID");
621     for (i = 0; i < mesh->n_elem; i++)
622       mesh->global_elem_ID[i] =
623           mesh->elem_ID[i * 2] + nelem_dist[mesh->elem_ID[i * 2 + 1]];
624 
625     HECMW_Barrier(VIS_COMM);
626     if ((h_table = (Hash_table *)HECMW_calloc(mesh->n_node * 4 - 3,
627                                               sizeof(Hash_table))) == NULL)
628       HECMW_vis_memory_exit("h_table");
629     for (i = 0; i < mesh->n_node * 4 - 3; i++) {
630       h_table[i].elemID    = 0;
631       h_table[i].next_elem = NULL;
632     }
633     index_connect = (int *)HECMW_calloc(mesh->n_elem + 1, sizeof(int));
634     if (index_connect == NULL) HECMW_vis_memory_exit("index_connect");
635     find_index_connectivity(mesh, index_connect);
636 
637     /*   if((connect=(int *)HECMW_calloc(index_connect[mesh->n_elem]*3,
638      * sizeof(int)))==NULL) */
639     if ((connect = (int *)HECMW_calloc(index_connect[mesh->n_elem] * 3 + 20,
640                                        sizeof(int))) ==
641         NULL) /* Modyfied S. Ito at 10/3/2007 */
642       HECMW_vis_memory_exit("connect");
643     /*  i*3: elem_id; i*3+1: face_id   i*3+2: pe_id */
644     for (i = 0; i < index_connect[mesh->n_elem] * 3; i++) connect[i] = -1;
645     build_hash_table(mesh, index_connect, h_table, index_face_tetra, face_tetra,
646                      index_face_prism, face_prism, index_face_hexa, face_hexa);
647 
648     build_connectivity(mesh, h_table, index_connect, connect, index_face_tetra,
649                        face_tetra, index_face_prism, face_prism,
650                        index_face_hexa, face_hexa);
651     h_free(h_table, mesh->n_node * 4 - 3);
652 
653     HECMW_Barrier(VIS_COMM);
654     for (i = 0; i < index_connect[mesh->n_elem]; i++) {
655       if (connect[i * 3] != -1) connect[i * 3 + 2] = mynode;
656     }
657 
658     if (mesh->n_neighbor_pe > 0) {
659       if ((export_no_element =
660                (int *)HECMW_calloc(mesh->n_neighbor_pe, sizeof(int))) == NULL)
661         HECMW_vis_memory_exit("export_no_element");
662       if ((import_no_element =
663                (int *)HECMW_calloc(mesh->n_neighbor_pe, sizeof(int))) == NULL)
664         HECMW_vis_memory_exit("import_no_element");
665       if ((index_import_element =
666                (int *)HECMW_calloc(mesh->n_neighbor_pe, sizeof(int))) == NULL)
667         HECMW_vis_memory_exit("import_no_element");
668 
669       for (i = 0; i < mesh->n_neighbor_pe; i++) export_no_element[i] = 0;
670       if ((export_elem = (int *)HECMW_calloc(mesh->n_neighbor_pe * mesh->n_elem,
671                                              sizeof(int))) == NULL)
672         HECMW_vis_memory_exit("export_elem");
673       if ((flag = (int *)HECMW_calloc(mesh->n_neighbor_pe, sizeof(int))) ==
674           NULL)
675         HECMW_vis_memory_exit("flag_connect");
676 
677       for (i           = 0; i < mesh->n_neighbor_pe * mesh->n_elem; i++)
678         export_elem[i] = -1;
679       for (i = 0; i < mesh->n_elem; i++) {
680         for (j = 0; j < mesh->n_neighbor_pe; j++) flag[j] = 0;
681         for (j = 0; j < mesh->elem_node_index[i + 1] - mesh->elem_node_index[i];
682              j++) {
683           nnode[j] = mesh->elem_node_item[mesh->elem_node_index[i] + j];
684           if (nnode[j] > mesh->nn_internal) {
685             for (k = 0; k < mesh->n_neighbor_pe; k++) {
686               if (flag[k] == 0) {
687                 /*			if(k==0) startnode=0;
688 else startnode=v->mesh->export_index[k-1];
689                  */
690                 startnode = mesh->import_index[k];
691                 for (m = startnode; m < mesh->import_index[k + 1]; m++) {
692                   if (mesh->import_item[m] == nnode[j]) {
693                     flag[k] = 1;
694                     export_no_element[k]++;
695                     export_elem[k * mesh->n_elem + export_no_element[k] - 1] =
696                         i;
697                     break;
698                   }
699                 }
700               }
701             }
702           }
703         }
704       }
705     } /* end of if n_neighbour_pe >0 */
706 
707     for (pe_no = 0; pe_no < pesize; pe_no++) {
708       if (mynode != pe_no) {
709         for (i = 0; i < mesh->n_neighbor_pe; i++) {
710           if (mesh->neighbor_pe[i] == pe_no) {
711             /* send export data */
712             if (export_no_element[i] > 0) {
713               /*                   fprintf(stderr, "the export number for PE %d
714                * in pe %d is %d\n", pe_no, mynode, export_no_element[i]);
715                */
716               if ((global_id = (int *)HECMW_calloc(export_no_element[i],
717                                                    sizeof(int))) == NULL)
718                 HECMW_vis_memory_exit("global_id");
719               index_a_connect =
720                   (int *)HECMW_calloc(export_no_element[i] + 1, sizeof(int));
721               if (index_a_connect == NULL)
722                 HECMW_vis_memory_exit("index_a_connect");
723               find_index_a_connect(mesh, export_no_element[i], i, export_elem,
724                                    index_a_connect);
725               if ((among_connect = (int *)HECMW_calloc(
726                        index_a_connect[export_no_element[i]] * 3,
727                        sizeof(int))) == NULL)
728                 HECMW_vis_memory_exit("among_connect");
729               for (j = 0; j < export_no_element[i]; j++) {
730                 global_id[j] =
731                     *(export_elem[i * mesh->n_elem + j] + mesh->global_elem_ID);
732                 for (k = 0; k < index_a_connect[j + 1] - index_a_connect[j];
733                      k++) {
734                   among_connect[(index_a_connect[j] + k) * 3] =
735                       connect[(index_connect[export_elem[i * n_elem + j]] + k) *
736                               3];
737                   among_connect[(index_a_connect[j] + k) * 3 + 1] =
738                       connect[(index_connect[export_elem[i * n_elem + j]] + k) *
739                                   3 +
740                               1];
741                   among_connect[(index_a_connect[j] + k) * 3 + 2] =
742                       connect[(index_connect[export_elem[i * n_elem + j]] + k) *
743                                   3 +
744                               2];
745                 }
746               }
747             }
748             HECMW_Send(&export_no_element[i], 1, HECMW_INT, pe_no, 0, VIS_COMM);
749             if (export_no_element[i] > 0) {
750               HECMW_Send(&index_a_connect[export_no_element[i]], 1, HECMW_INT,
751                          pe_no, 0, VIS_COMM);
752               HECMW_Send(global_id, export_no_element[i], HECMW_INT, pe_no, 0,
753                          VIS_COMM);
754               HECMW_Send(index_a_connect, export_no_element[i] + 1, HECMW_INT,
755                          pe_no, 0, VIS_COMM);
756               HECMW_Send(among_connect,
757                          index_a_connect[export_no_element[i]] * 3, HECMW_INT,
758                          pe_no, 0, VIS_COMM);
759 
760               HECMW_free(global_id);
761               HECMW_free(index_a_connect);
762               HECMW_free(among_connect);
763             }
764           } /* end of if ==pe_no */
765         }   /* end of for i */
766       }     /* end of if mynode !=pe_no */
767       /* else if (mynode==pe_no) { */
768       else if (mynode == pe_no &&
769                pesize > 1) { /* Modified by S. Ito at 10/2/2007 */
770         /* receive overlapped element data from other pes */
771         if ((g_id = (int **)HECMW_calloc(mesh->n_neighbor_pe, sizeof(int *))) ==
772             NULL)
773           HECMW_vis_memory_exit("g_id");
774 
775         if ((a_connect = (int **)HECMW_calloc(mesh->n_neighbor_pe,
776                                               sizeof(int *))) == NULL)
777           HECMW_vis_memory_exit("a_connect");
778         if ((index_a = (int **)HECMW_calloc(mesh->n_neighbor_pe,
779                                             sizeof(int *))) == NULL)
780           HECMW_vis_memory_exit("index_a");
781 
782         for (i = 0; i < mesh->n_neighbor_pe; i++) {
783           j = mesh->neighbor_pe[i];
784           HECMW_Recv(&import_no_element[i], 1, HECMW_INT, j, HECMW_ANY_TAG,
785                      VIS_COMM, &stat);
786           if (import_no_element[i] > 0) {
787             HECMW_Recv(&index_import_element[i], 1, HECMW_INT, j, HECMW_ANY_TAG,
788                        VIS_COMM, &stat);
789             if ((g_id[i] = (int *)HECMW_calloc(import_no_element[i],
790                                                sizeof(int))) == NULL)
791               HECMW_vis_memory_exit("g_id");
792 
793             HECMW_Recv(g_id[i], import_no_element[i], HECMW_INT, j,
794                        HECMW_ANY_TAG, VIS_COMM, &stat);
795             if ((index_a[i] = (int *)HECMW_calloc(import_no_element[i] + 1,
796                                                   sizeof(int))) == NULL)
797               HECMW_vis_memory_exit("index_a");
798 
799             if ((a_connect[i] = (int *)HECMW_calloc(index_import_element[i] * 3,
800                                                     sizeof(int))) == NULL)
801               HECMW_vis_memory_exit("a_connect");
802             HECMW_Recv(index_a[i], import_no_element[i] + 1, HECMW_INT, j,
803                        HECMW_ANY_TAG, VIS_COMM, &stat);
804 
805             HECMW_Recv(a_connect[i], index_import_element[i] * 3, HECMW_INT, j,
806                        HECMW_ANY_TAG, VIS_COMM, &stat);
807           }
808         }
809         /* recompute boundary element */
810         for (i = 0; i < mesh->n_elem; i++) {
811           for (j = 0; j < index_connect[i + 1] - index_connect[i]; j++) {
812             if (connect[(index_connect[i] + j) * 3] == -1) {
813               flag1 = 1;
814               for (k = 0; k < mesh->n_neighbor_pe; k++) {
815                 if (flag1 == 1) {
816                   for (m = 0; m < import_no_element[k]; m++) {
817                     if ((*(g_id[k] + m) == mesh->global_elem_ID[i]) &&
818                         (*(a_connect[k] + (*(index_a[k] + m) + j) * 3) != -1)) {
819                       connect[(index_connect[i] + j) * 3] =
820                           *(a_connect[k] + (*(index_a[k] + m) + j) * 3);
821                       connect[(index_connect[i] + j) * 3 + 1] =
822                           *(a_connect[k] + (*(index_a[k] + m) + j) * 3 + 1);
823                       connect[(index_connect[i] + j) * 3 + 2] =
824                           *(a_connect[k] + (*(index_a[k] + m) + j) * 3 + 2);
825 
826                       flag1 = 0;
827                       break;
828                     }
829                   }
830                 }
831               }
832             }
833           }
834         }
835         HECMW_free(a_connect);
836         HECMW_free(g_id);
837         HECMW_free(index_a);
838       } /* end of if mynode==pe_no */
839       /* finish in one pe_no */
840       HECMW_Barrier(VIS_COMM);
841     } /* end of for pe_no */
842       /*	  fprintf(stderr,"Finish the connectivity build among PES\n");
843        */
844   }
845   /*  b_patch=(Boundary_patch *)HECMW_malloc(sizeof(Boundary_patch));
846 if(b_patch==NULL)
847     HECMW_vis_memory_exit("b_patch");
848 p1=b_patch;
849    */
850   if (init_flag != 1) {
851     index_connect = global_connect->index_connect;
852     connect       = global_connect->connect;
853   }
854 
855   b_point      = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
856   head_b_patch = (Head_patch_tetra *)HECMW_malloc(sizeof(Head_patch_tetra));
857   if ((b_point == NULL) || (head_b_patch == NULL))
858     HECMW_vis_memory_exit("boundary: initialization");
859   b_point->ident          = 0;
860   head_b_patch->num_patch = 0;
861 
862   node_hit = (int *)HECMW_calloc(mesh->n_node, sizeof(int));
863   if (node_hit == NULL) HECMW_vis_memory_exit("node_hit");
864   for (i = 0; i < mesh->n_node; i++) node_hit[i] = -1;
865   /*  if(sff->deform_display_on==1) {
866     flag_deform=0;
867     for(i=0;i<data->nn_component;i++) {
868   if(strncmp("DISPLACEMENT", data->node_label[i], 10)==0) {
869                   flag_deform=1;
870                   sff->disp_comp=i;
871                   break;
872           }
873     }
874     if(flag_deform==0)
875             HECMW_vis_print_exit("No displament data for displaying deform");
876 }
877    */
878 
879   n_elem = mesh->n_elem;
880   c_base = 0;
881   for (i = 0; i < sff->color_comp; i++) {
882     c_base += data->nn_dof[i];
883   }
884   d_base = 0;
885   if (sff->deform_display_on == 1) {
886     for (i = 0; i < sff->disp_comp; i++) d_base += data->nn_dof[i];
887   }
888   *minc = 1.0E17;
889   *maxc = -1.0E17;
890   for (i = 0; i < n_elem; i++)
891     if (mesh->elem_ID[i * 2 + 1] == mynode) {
892       for (j = 0; j < index_connect[i + 1] - index_connect[i]; j++) {
893         if (connect[(index_connect[i] + j) * 3] == -1) {
894           if ((mesh->elem_type[i] == 341) || (mesh->elem_type[i] == 342) ||
895               (mesh->elem_type[i] == 3414)) {
896             for (k = 0; k < index_face_tetra[j + 1] - index_face_tetra[j]; k++)
897               node[k] =
898                   mesh->elem_node_item[mesh->elem_node_index[i] +
899                                        face_tetra[index_face_tetra[j] + k]];
900             add_one_patch(sff, mesh, data, node_hit, b_point, head_b_patch,
901                           node, c_base, d_base, tn_component);
902           } else if ((mesh->elem_type[i] == 351) ||
903                      (mesh->elem_type[i] == 352)) {
904             for (k = 0; k < index_face_prism[j + 1] - index_face_prism[j]; k++)
905               node[k] =
906                   mesh->elem_node_item[mesh->elem_node_index[i] +
907                                        face_prism[index_face_prism[j] + k]];
908             if (j < 3)
909               add_two_patch(sff, mesh, data, node_hit, b_point, head_b_patch,
910                             node, c_base, d_base, tn_component);
911             else if (j >= 3)
912               add_one_patch(sff, mesh, data, node_hit, b_point, head_b_patch,
913                             node, c_base, d_base, tn_component);
914           } else if ((mesh->elem_type[i] == 361) ||
915                      (mesh->elem_type[i] == 362)) {
916             for (k    = 0; k < index_face_hexa[j + 1] - index_face_hexa[j]; k++)
917               node[k] = mesh->elem_node_item[mesh->elem_node_index[i] +
918                                              face_hexa[index_face_hexa[j] + k]];
919             add_two_patch(sff, mesh, data, node_hit, b_point, head_b_patch,
920                           node, c_base, d_base, tn_component);
921           }
922         }
923       }
924     }
925   HECMW_free(node_hit);
926   if (head_b_patch->num_patch > 0) {
927     result[sf_i].n_vertex = b_point->ident;
928     result[sf_i].n_patch  = head_b_patch->num_patch;
929     result[sf_i].vertex =
930         (double *)HECMW_calloc(result[sf_i].n_vertex * 3, sizeof(double));
931     result[sf_i].color =
932         (double *)HECMW_calloc(result[sf_i].n_vertex, sizeof(double));
933     result[sf_i].patch =
934         (int *)HECMW_calloc(result[sf_i].n_patch * 3, sizeof(int));
935     if ((result[sf_i].vertex == NULL) || (result[sf_i].patch == NULL) ||
936         (result[sf_i].color == NULL))
937       HECMW_vis_memory_exit("Boundary: vertex, patch and color");
938     if (sff->deform_display_on == 1) {
939       result[sf_i].disp =
940           (double *)HECMW_calloc(result[sf_i].n_vertex * 3, sizeof(double));
941       if (result[sf_i].disp == NULL) HECMW_vis_memory_exit("Boundary: disp");
942     }
943   }
944   if (b_point->ident > 0) {
945     p1 = b_point->nextpoint;
946     for (i = 0; i < b_point->ident; i++) {
947       for (j                                     = 0; j < 3; j++)
948         result[sf_i].vertex[(p1->ident) * 3 + j] = p1->geom[j];
949       result[sf_i].color[p1->ident]              = p1->cdata;
950       if (p1->cdata < *minc) *minc               = p1->cdata;
951       if (p1->cdata > *maxc) *maxc               = p1->cdata;
952       if (sff->deform_display_on == 1) {
953         for (j                                   = 0; j < 3; j++)
954           result[sf_i].disp[(p1->ident) * 3 + j] = p1->disp[j];
955       }
956 
957       p2 = p1;
958       p1 = p1->nextpoint;
959       HECMW_free(p2);
960     }
961     HECMW_free(b_point);
962   }
963   if (head_b_patch->num_patch > 0) {
964     t1 = head_b_patch->patch_link;
965     for (i = 0; i < head_b_patch->num_patch; i++) {
966       for (j                          = 0; j < 3; j++)
967         result[sf_i].patch[i * 3 + j] = *tvertex + t1->patch[j];
968       t2                              = t1;
969       t1                              = t1->next_patch;
970       HECMW_free(t2);
971     }
972     HECMW_free(head_b_patch);
973   }
974 
975   *tvertex += result[sf_i].n_vertex;
976   *tpatch += result[sf_i].n_patch;
977   /*          HECMW_free(index_connect);
978     HECMW_free(connect);
979    */
980   if (init_flag == 1) {
981     global_connect->index_connect = index_connect;
982     global_connect->connect       = connect;
983   }
984   return;
985 }
986