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