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_tetra_intersect.h"
7 
8 #include <math.h>
9 #include "hecmw_vis_mem_util.h"
10 #include "hecmw_vis_case_table.h"
11 #include "hecmw_malloc.h"
12 
get_tetra_data(Surface * sff,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int elemID,Tetra * tetra,int tn_component)13 int get_tetra_data(Surface *sff, struct hecmwST_local_mesh *mesh,
14                    struct hecmwST_result_data *data, int elemID, Tetra *tetra,
15                    int tn_component) {
16   int i, j;
17   int nodeID;
18   int c_base, s_base;
19   double tmp;
20   c_base = 0;
21 
22   if (sff->display_way != 4) {
23     for (i = 0; i < sff->color_comp; i++) {
24       c_base += data->nn_dof[i];
25     }
26   }
27   if (sff->surface_style == 2) {
28     s_base = 0;
29     for (i = 0; i < sff->data_comp; i++) s_base += data->nn_dof[i];
30   }
31 
32   /*  set field data of voxel in cube  */
33   for (i = 0; i < 4; i++) {
34     nodeID = mesh->elem_node_item[mesh->elem_node_index[elemID] + i];
35     tetra->local_vid[i] = nodeID - 1;
36     if (sff->display_way != 4) {
37       if ((data->nn_dof[sff->color_comp] > 1) && (sff->color_subcomp == 0)) {
38         tetra->c_data[i] = 0.0;
39         for (j = 0; j < data->nn_dof[sff->color_comp]; j++) {
40           tmp = data->node_val_item[(nodeID - 1) * tn_component + c_base + j];
41           tetra->c_data[i] += tmp * tmp;
42         }
43         tetra->c_data[i] = sqrt(tetra->c_data[i]);
44       } else if ((data->nn_dof[sff->color_comp] > 1) &&
45                  (sff->color_subcomp != 0))
46         tetra->c_data[i] =
47             data->node_val_item[(nodeID - 1) * tn_component + c_base +
48                                 (sff->color_subcomp - 1)];
49       else if (data->nn_dof[sff->color_comp] == 1)
50         tetra->c_data[i] =
51             data->node_val_item[(nodeID - 1) * tn_component + c_base];
52 
53     } else if (sff->display_way == 4)
54       tetra->c_data[i] = sff->specified_color;
55 
56     tetra->axis[i * 3]     = mesh->node[(nodeID - 1) * 3];
57     tetra->axis[i * 3 + 1] = mesh->node[(nodeID - 1) * 3 + 1];
58     tetra->axis[i * 3 + 2] = mesh->node[(nodeID - 1) * 3 + 2];
59 
60     if (sff->surface_style == 2) {
61       if ((data->nn_dof[sff->data_comp] > 1) && (sff->data_subcomp == 0)) {
62         tetra->s_data[i] = 0.0;
63         for (j = 0; j < data->nn_dof[sff->data_comp]; j++) {
64           tmp = data->node_val_item[(nodeID - 1) * tn_component + s_base + j];
65           tetra->s_data[i] += tmp * tmp;
66         }
67         tetra->s_data[i] = sqrt(tetra->s_data[i]);
68       } else if ((data->nn_dof[sff->data_comp] > 1) && (sff->data_subcomp != 0))
69         tetra->s_data[i] =
70             data->node_val_item[(nodeID - 1) * tn_component + s_base +
71                                 (sff->data_subcomp - 1)];
72       else if (data->nn_dof[sff->data_comp] == 1)
73         tetra->s_data[i] =
74             data->node_val_item[(nodeID - 1) * tn_component + s_base];
75     } else if (sff->surface_style == 3)
76 
77       tetra->s_data[i] =
78           get_value_equ(sff->cont_equ, sff->cross_type, tetra->axis[i * 3],
79                         tetra->axis[i * 3 + 1], tetra->axis[i * 3 + 2]);
80   }
81 
82   tetra->elem_id[0] = mesh->elem_ID[elemID * 2];
83   tetra->elem_id[1] = mesh->elem_ID[elemID * 2 + 1];
84 
85   return 1;
86 }
87 
intersect_line(int v0,int v1,double isovalue,Tetra * tetra,double point[3])88 static double intersect_line(int v0, int v1, double isovalue, Tetra *tetra,
89                              double point[3]) {
90   int j;
91   double rate, color;
92 
93   if (fabs(tetra->s_data[v1] - tetra->s_data[v0]) < EPSILON)
94     HECMW_vis_print_exit("There is something wrong in data precison");
95   else {
96     rate = (isovalue - tetra->s_data[v0]) /
97            (tetra->s_data[v1] - tetra->s_data[v0]);
98     for (j     = 0; j < 3; j++)
99       point[j] = rate * (tetra->axis[v1 * 3 + j] - tetra->axis[v0 * 3 + j]) +
100                  tetra->axis[v0 * 3 + j];
101     color = rate * (tetra->c_data[v1] - tetra->c_data[v0]) + tetra->c_data[v0];
102   }
103   return (color);
104 }
105 
find_intersection_tetra(Tetra * tetra,double isovalue,Tetra_point * tetra_point,Head_patch_tetra * head_patch_tetra,Hash_vertex * vertex_hash_table)106 void find_intersection_tetra(Tetra *tetra, double isovalue,
107                              Tetra_point *tetra_point,
108                              Head_patch_tetra *head_patch_tetra,
109                              Hash_vertex *vertex_hash_table) {
110   int i, j, k;
111   int num_gt_0, v0_id;
112   int tmp_int;
113   double point[3], color;
114   Hash_vertex *h1, *h2;
115   Tetra_point *p1, *p2;
116   Patch_tetra *t1, *t2;
117   int flag_existing, index_patch, patch[4];
118   int v1, v2, v[4], count_minus, count_plus;
119   double fp[4][3], n_f[3], n_norm, f_cen_point[3], sign_f, n_c, c_c[3];
120   num_gt_0 = 0;
121   for (i = 0; i < 4; i++) {
122     if (tetra->s_data[i] > isovalue) num_gt_0++;
123   }
124 
125   if ((num_gt_0 != 0) &&
126       (num_gt_0 != 4)) { /* There are patches inside the tetra */
127 
128     if ((num_gt_0 == 1) ||
129         (num_gt_0 == 3)) { /* There is one patch inside the tetra */
130       t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
131       if (t1 == NULL) HECMW_vis_memory_exit("t1");
132       t2 = head_patch_tetra->patch_link;
133       head_patch_tetra->num_patch++;
134       head_patch_tetra->patch_link = t1;
135       t1->next_patch               = t2;
136 
137       if (num_gt_0 == 1) {
138         v0_id = -1;
139         for (i = 0; i < 4; i++) {
140           if (tetra->s_data[i] > isovalue) v0_id = i;
141         }
142       } else if (num_gt_0 == 3) {
143         v0_id = -1;
144         for (i = 0; i < 4; i++) {
145           if (tetra->s_data[i] <= isovalue) v0_id = i;
146         }
147       }
148       /* find intersection */
149       index_patch = -1;
150       for (i = 0; i < 4; i++) {
151         if (i != v0_id) {
152           color   = intersect_line(v0_id, i, isovalue, tetra, point);
153           tmp_int = tetra->local_vid[v0_id] + tetra->local_vid[i];
154           /*					if(vertex_hash_table[tmp_int].ident==0)
155           {
156                   vertex_hash_table[tmp_int].ident++;
157                   h1=(Hash_vertex *)HECMW_malloc(sizeof(Hash_vertex));
158                   if(h1==NULL)
159                           HECMW_vis_memory_exit("h1");
160               vertex_hash_table[tmp_int].next_vertex=h1;
161                   h1->next_vertex=NULL;
162                   h1->ident=tetra_point->ident;
163                   for(j=0;j<3;j++)
164                           h1->geom[j]=point[j];
165                   p1=(Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
166                   if(p1==NULL)
167                           HECMW_vis_memory_exit("p1");
168                   tetra_point->ident++;
169                   p2=tetra_point->nextpoint;
170                   tetra_point->nextpoint=p1;
171                   p1->cdata=color;
172                   p1->ident=tetra_point->ident-1;
173                   for(j=0;j<3;j++)
174                           p1->geom[j]=point[j];
175                   p1->nextpoint=p2;
176                   index_patch++;
177                   t1->patch[index_patch]=p1->ident;
178           }
179           else if(vertex_hash_table[tmp_int].ident>0) { */
180           flag_existing = 0;
181           h1            = vertex_hash_table[tmp_int].next_vertex;
182           for (j = 0; j < vertex_hash_table[tmp_int].ident; j++) {
183             if ((fabs(h1->geom[0] - point[0]) < EPSILON) &&
184                 (fabs(h1->geom[1] - point[1]) < EPSILON) &&
185                 (fabs(h1->geom[2] - point[2]) < EPSILON)) {
186               flag_existing = 1;
187               index_patch++;
188               patch[index_patch] = h1->ident;
189               for (k = 0; k < 3; k++) fp[index_patch][k] = point[k];
190               break;
191             }
192             h1 = h1->next_vertex;
193           }
194           if (flag_existing == 0) { /*adding new vertex */
195             vertex_hash_table[tmp_int].ident++;
196             h1 = (Hash_vertex *)HECMW_malloc(sizeof(Hash_vertex));
197             if (h1 == NULL) HECMW_vis_memory_exit("h1");
198             h2 = vertex_hash_table[tmp_int].next_vertex;
199             vertex_hash_table[tmp_int].next_vertex = h1;
200             h1->next_vertex                        = h2;
201             h1->ident                              = tetra_point->ident;
202             for (j = 0; j < 3; j++) h1->geom[j] = point[j];
203             p1 = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
204             if (p1 == NULL) HECMW_vis_memory_exit("p1");
205             tetra_point->ident++;
206             p2                     = tetra_point->nextpoint;
207             tetra_point->nextpoint = p1;
208             p1->cdata              = color;
209             p1->ident              = tetra_point->ident - 1;
210             for (j = 0; j < 3; j++) p1->geom[j] = point[j];
211             p1->nextpoint                       = p2;
212             index_patch++;
213             patch[index_patch] = p1->ident;
214             for (k = 0; k < 3; k++) fp[index_patch][k] = point[k];
215           }
216                                         /*					}
217 					 */				}
218 
219       } /* end for 0, 4 */
220       /*  judge whether the rotation direction pointed out to isosurface */
221 
222       n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
223                (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
224       n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
225                (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
226       n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
227                (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
228       for (j = 0; j < 3; j++) fp[3][j] = tetra->axis[v0_id * 3 + j];
229       n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
230       if (fabs(n_norm) > EPSILON) {
231         for (j = 0; j < 3; j++) n_f[j] /= n_norm;
232       }
233       /*selce the direction point to inside the element */
234       for (j           = 0; j < 3; j++)
235         f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j]) / 3.0;
236       for (j = 0; j < 3; j++) c_c[j] = fp[3][j] - f_cen_point[j];
237       n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
238       if (fabs(n_c) > EPSILON) {
239         for (j = 0; j < 3; j++) c_c[j] /= n_c;
240       }
241       sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
242       if (((tetra->s_data[v0_id] > isovalue) && (sign_f < -EPSILON)) ||
243           ((tetra->s_data[v0_id] <= isovalue) && (sign_f > EPSILON))) {
244         tmp_int  = patch[1];
245         patch[1] = patch[2];
246         patch[2] = tmp_int;
247       }
248       for (j = 0; j < 3; j++) t1->patch[j] = patch[j];
249 
250     } /* end if (num_gt_0==1 or 3) */
251 
252     else if (num_gt_0 == 2) { /* Two patches inside the tetra */
253       index_patch = -1;
254       count_minus = 0;
255       count_plus  = 0;
256       for (i = 0; i < 4; i++) {
257         if (tetra->s_data[i] <= isovalue) {
258           v[count_minus] = i;
259           count_minus++;
260         } else {
261           v[2 + count_plus] = i;
262           count_plus++;
263         }
264       }
265 #ifdef old
266       /*   judge the tetrahedra whether in the right rotation side */
267       for (i = 0; i < 4; i++)
268         for (j = 0; j < 3; j++) fp[i][j] = tetra->axis[v[i] * 3 + j];
269       n_f[0] = (fp[2][1] - fp[1][1]) * (fp[1][2] - fp[0][2]) -
270                (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[1][2]);
271       n_f[1] = -(fp[2][0] - fp[1][0]) * (fp[1][2] - fp[0][2]) +
272                (fp[1][0] - fp[0][0]) * (fp[2][2] - fp[1][2]);
273       n_f[2] = (fp[2][0] - fp[1][0]) * (fp[1][1] - fp[0][1]) -
274                (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[1][1]);
275       n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
276       if (fabs(n_norm) > EPSILON) {
277         for (j = 0; j < 3; j++) n_f[j] /= n_norm;
278       }
279       /*selce the direction point to inside the element */
280       for (j           = 0; j < 3; j++)
281         f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j]) / 3.0;
282       for (j = 0; j < 3; j++) c_c[j] = fp[3][j] - f_cen_point[j];
283       n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
284       if (fabs(n_c) > EPSILON) {
285         for (j = 0; j < 3; j++) c_c[j] /= n_c;
286       }
287       sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
288       if (sign_f < -EPSILON) {
289         tmp_int = v[2];
290         v[2]    = v[3];
291         v[3]    = tmp_int;
292       }
293 #endif
294 
295       for (i = 0; i < 4; i++) {
296         if (i == 0) {
297           v1 = v[0];
298           v2 = v[2];
299         } else if (i == 1) {
300           v1 = v[1];
301           v2 = v[2];
302         } else if (i == 2) {
303           v1 = v[1];
304           v2 = v[3];
305         } else if (i == 3) {
306           v1 = v[0];
307           v2 = v[3];
308         }
309 
310         color   = intersect_line(v1, v2, isovalue, tetra, point);
311         tmp_int = tetra->local_vid[v1] + tetra->local_vid[v2];
312         /*					if(vertex_hash_table[tmp_int].ident==0)
313            {
314                         vertex_hash_table[tmp_int].ident++;
315                         h1=(Hash_vertex *)HECMW_malloc(sizeof(Hash_vertex));
316                         if(h1==NULL)
317                                 HECMW_vis_memory_exit("h1");
318                     vertex_hash_table[tmp_int].next_vertex=h1;
319                         h1->next_vertex=NULL;
320                         h1->ident=tetra_point->ident;
321                         for(j=0;j<3;j++)
322                                 h1->geom[j]=point[j];
323                         p1=(Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
324                         if(p1==NULL)
325                                 HECMW_vis_memory_exit("p1");
326                         tetra_point->ident++;
327                         p2=tetra_point->nextpoint;
328                         tetra_point->nextpoint=p1;
329                         p1->cdata=color;
330                         p1->ident=tetra_point->ident-1;
331                         for(j=0;j<3;j++)
332                                 p1->geom[j]=point[j];
333                         p1->nextpoint=p2;
334                         index_patch++;
335                         patch[index_patch]=p1->ident;
336                 }
337                 else if(vertex_hash_table[tmp_int].ident>0) {
338          */
339         flag_existing = 0;
340         h1            = vertex_hash_table[tmp_int].next_vertex;
341         for (j = 0; j < vertex_hash_table[tmp_int].ident; j++) {
342           if ((fabs(h1->geom[0] - point[0]) < EPSILON) &&
343               (fabs(h1->geom[1] - point[1]) < EPSILON) &&
344               (fabs(h1->geom[2] - point[2]) < EPSILON)) {
345             flag_existing = 1;
346             index_patch++;
347             patch[index_patch] = h1->ident;
348             for (k = 0; k < 3; k++) fp[index_patch][k] = point[k];
349             break;
350           }
351           h1 = h1->next_vertex;
352         }
353         if (flag_existing == 0) { /*adding new vertex */
354           vertex_hash_table[tmp_int].ident++;
355           h1 = (Hash_vertex *)HECMW_malloc(sizeof(Hash_vertex));
356           if (h1 == NULL) HECMW_vis_memory_exit("h1");
357           h2 = vertex_hash_table[tmp_int].next_vertex;
358           vertex_hash_table[tmp_int].next_vertex = h1;
359           h1->next_vertex                        = h2;
360           h1->ident                              = tetra_point->ident;
361           for (j = 0; j < 3; j++) h1->geom[j] = point[j];
362           p1 = (Tetra_point *)HECMW_malloc(sizeof(Tetra_point));
363           if (p1 == NULL) HECMW_vis_memory_exit("p1");
364           tetra_point->ident++;
365           p2                     = tetra_point->nextpoint;
366           tetra_point->nextpoint = p1;
367           p1->cdata              = color;
368           p1->ident              = tetra_point->ident - 1;
369           for (j = 0; j < 3; j++) p1->geom[j] = point[j];
370           p1->nextpoint                       = p2;
371           index_patch++;
372           patch[index_patch] = p1->ident;
373           for (k = 0; k < 3; k++) fp[index_patch][k] = point[k];
374           /*						}
375           }
376            */
377         }
378       } /* end for i*/
379 
380       /*  judge whether the rotation direction pointed out to isosurface */
381 
382       n_f[0] = (fp[1][1] - fp[0][1]) * (fp[2][2] - fp[0][2]) -
383                (fp[2][1] - fp[0][1]) * (fp[1][2] - fp[0][2]);
384       n_f[1] = -(fp[1][0] - fp[0][0]) * (fp[2][2] - fp[0][2]) +
385                (fp[2][0] - fp[0][0]) * (fp[1][2] - fp[0][2]);
386       n_f[2] = (fp[1][0] - fp[0][0]) * (fp[2][1] - fp[0][1]) -
387                (fp[2][0] - fp[0][0]) * (fp[1][1] - fp[0][1]);
388       for (j = 0; j < 3; j++) fp[3][j] = tetra->axis[0 * 3 + j];
389       n_norm = sqrt(n_f[0] * n_f[0] + n_f[1] * n_f[1] + n_f[2] * n_f[2]);
390       if (fabs(n_norm) > EPSILON) {
391         for (j = 0; j < 3; j++) n_f[j] /= n_norm;
392       }
393       /*selce the direction point to inside the element */
394       for (j           = 0; j < 3; j++)
395         f_cen_point[j] = (fp[0][j] + fp[1][j] + fp[2][j]) / 3.0;
396       for (j = 0; j < 3; j++) c_c[j] = fp[3][j] - f_cen_point[j];
397       n_c = sqrt(c_c[0] * c_c[0] + c_c[1] * c_c[1] + c_c[2] * c_c[2]);
398       if (fabs(n_c) > EPSILON) {
399         for (j = 0; j < 3; j++) c_c[j] /= n_c;
400       }
401       sign_f = n_f[0] * c_c[0] + n_f[1] * c_c[1] + n_f[2] * c_c[2];
402       if (((tetra->s_data[0] > isovalue) && (sign_f < -EPSILON)) ||
403           ((tetra->s_data[0] <= isovalue) && (sign_f > EPSILON))) {
404         tmp_int  = patch[1];
405         patch[1] = patch[3];
406         patch[3] = tmp_int;
407       }
408 
409       t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
410       if (t1 == NULL) HECMW_vis_memory_exit("t1");
411       t2 = head_patch_tetra->patch_link;
412       head_patch_tetra->num_patch++;
413       head_patch_tetra->patch_link = t1;
414       t1->next_patch               = t2;
415       for (j = 0; j < 3; j++) t1->patch[j] = patch[j];
416       t1 = (Patch_tetra *)HECMW_malloc(sizeof(Patch_tetra));
417       if (t1 == NULL) HECMW_vis_memory_exit("t1");
418       t2 = head_patch_tetra->patch_link;
419       head_patch_tetra->num_patch++;
420       head_patch_tetra->patch_link = t1;
421       t1->next_patch               = t2;
422       t1->patch[0]                 = patch[0];
423       t1->patch[1]                 = patch[2];
424       t1->patch[2]                 = patch[3];
425 
426     } /* end if (num_gt_0==2) */
427   }
428   return;
429 }
430