FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_vis_calc_attr.c
Go to the documentation of this file.
1/*****************************************************************************
2 * Copyright (c) 2019 FrontISTR Commons
3 * This software is released under the MIT License, see LICENSE.txt
4 *****************************************************************************/
5
7
8#include <math.h>
9
10/*----------------------------------------------------------------------
11# Subroutines in this file on isosurface generation by Marching Cubes is
12based
13 on the revision of Dr. Yuriko Takeshima's codes when she was working
14part time in RIST
15#---------------------------------------------------------------------- */
16
17/* calculate the geometry at the point in cube */
18int get_point_geom(int point_index, Cell *cell, double fvalue,
19 Fgeom *point_geom, double *cdata, int disamb_flag) {
20 if (point_index < 100) {
21 get_edgepoint(point_index, cell, fvalue, point_geom, cdata);
22 } else if (point_index < 200) {
23 get_insidepoint((point_index - 100), cell, fvalue, point_geom, cdata,
24 disamb_flag);
25 } else {
26 get_gridpoint((point_index - 200), cell, point_geom, cdata);
27 }
28
29 return 1;
30}
31
32/* return the geometry at the vertex of cube */
33void get_gridpoint(int voxel_index, Cell *cell, Fgeom *vert_geom,
34 double *cdata) {
35 vert_geom->x = cell->axis[voxel_index * 3];
36 vert_geom->y = cell->axis[voxel_index * 3 + 1];
37 vert_geom->z = cell->axis[voxel_index * 3 + 2];
38 *cdata = cell->c_data[voxel_index];
39}
40
41void get_edgepoint(int edge_index, Cell *cell, double fvalue, Fgeom *vert_geom,
42 double *cdata) {
43 double flip_ratio, ratio;
44
45 switch (edge_index) {
46 case 0:
47 ratio = linear_interpolate(cell->s_data[0], cell->s_data[1], fvalue);
48 flip_ratio = 1 - ratio;
49 vert_geom->x = flip_ratio * cell->axis[0] + ratio * cell->axis[3];
50 vert_geom->y = flip_ratio * cell->axis[1] + ratio * cell->axis[4];
51 vert_geom->z = flip_ratio * cell->axis[2] + ratio * cell->axis[5];
52 *cdata = flip_ratio * cell->c_data[0] + ratio * cell->c_data[1];
53 break;
54 case 1:
55 ratio = linear_interpolate(cell->s_data[1], cell->s_data[2], fvalue);
56 flip_ratio = 1 - ratio;
57 vert_geom->x = flip_ratio * cell->axis[3] + ratio * cell->axis[6];
58 vert_geom->y = flip_ratio * cell->axis[4] + ratio * cell->axis[7];
59 vert_geom->z = flip_ratio * cell->axis[5] + ratio * cell->axis[8];
60 *cdata = flip_ratio * cell->c_data[1] + ratio * cell->c_data[2];
61 break;
62 case 2:
63 ratio = linear_interpolate(cell->s_data[3], cell->s_data[2], fvalue);
64 flip_ratio = 1 - ratio;
65 vert_geom->x = flip_ratio * cell->axis[9] + ratio * cell->axis[6];
66 vert_geom->y = flip_ratio * cell->axis[10] + ratio * cell->axis[7];
67 vert_geom->z = flip_ratio * cell->axis[11] + ratio * cell->axis[8];
68 *cdata = flip_ratio * cell->c_data[3] + ratio * cell->c_data[2];
69 break;
70 case 3:
71 ratio = linear_interpolate(cell->s_data[0], cell->s_data[3], fvalue);
72 flip_ratio = 1 - ratio;
73 vert_geom->x = flip_ratio * cell->axis[0] + ratio * cell->axis[9];
74 vert_geom->y = flip_ratio * cell->axis[1] + ratio * cell->axis[10];
75 vert_geom->z = flip_ratio * cell->axis[2] + ratio * cell->axis[11];
76 *cdata = flip_ratio * cell->c_data[0] + ratio * cell->c_data[3];
77 break;
78 case 4:
79 ratio = linear_interpolate(cell->s_data[4], cell->s_data[5], fvalue);
80 flip_ratio = 1 - ratio;
81 vert_geom->x = flip_ratio * cell->axis[12] + ratio * cell->axis[15];
82 vert_geom->y = flip_ratio * cell->axis[13] + ratio * cell->axis[16];
83 vert_geom->z = flip_ratio * cell->axis[14] + ratio * cell->axis[17];
84 *cdata = flip_ratio * cell->c_data[4] + ratio * cell->c_data[5];
85 break;
86 case 5:
87 ratio = linear_interpolate(cell->s_data[5], cell->s_data[6], fvalue);
88 flip_ratio = 1 - ratio;
89 vert_geom->x = flip_ratio * cell->axis[15] + ratio * cell->axis[18];
90 vert_geom->y = flip_ratio * cell->axis[16] + ratio * cell->axis[19];
91 vert_geom->z = flip_ratio * cell->axis[17] + ratio * cell->axis[20];
92 *cdata = flip_ratio * cell->c_data[5] + ratio * cell->c_data[6];
93 break;
94 case 6:
95 ratio = linear_interpolate(cell->s_data[7], cell->s_data[6], fvalue);
96 flip_ratio = 1 - ratio;
97 vert_geom->x = flip_ratio * cell->axis[21] + ratio * cell->axis[18];
98 vert_geom->y = flip_ratio * cell->axis[22] + ratio * cell->axis[19];
99 vert_geom->z = flip_ratio * cell->axis[23] + ratio * cell->axis[20];
100 *cdata = flip_ratio * cell->c_data[7] + ratio * cell->c_data[6];
101 break;
102 case 7:
103 ratio = linear_interpolate(cell->s_data[4], cell->s_data[7], fvalue);
104 flip_ratio = 1 - ratio;
105 vert_geom->x = flip_ratio * cell->axis[12] + ratio * cell->axis[21];
106 vert_geom->y = flip_ratio * cell->axis[13] + ratio * cell->axis[22];
107 vert_geom->z = flip_ratio * cell->axis[14] + ratio * cell->axis[23];
108 *cdata = flip_ratio * cell->c_data[4] + ratio * cell->c_data[7];
109 break;
110 case 8:
111 ratio = linear_interpolate(cell->s_data[0], cell->s_data[4], fvalue);
112 flip_ratio = 1 - ratio;
113 vert_geom->x = flip_ratio * cell->axis[0] + ratio * cell->axis[12];
114 vert_geom->y = flip_ratio * cell->axis[1] + ratio * cell->axis[13];
115 vert_geom->z = flip_ratio * cell->axis[2] + ratio * cell->axis[14];
116 *cdata = flip_ratio * cell->c_data[0] + ratio * cell->c_data[4];
117 break;
118 case 9:
119 ratio = linear_interpolate(cell->s_data[1], cell->s_data[5], fvalue);
120 flip_ratio = 1 - ratio;
121 vert_geom->x = flip_ratio * cell->axis[3] + ratio * cell->axis[15];
122 vert_geom->y = flip_ratio * cell->axis[4] + ratio * cell->axis[16];
123 vert_geom->z = flip_ratio * cell->axis[5] + ratio * cell->axis[17];
124 *cdata = flip_ratio * cell->c_data[1] + ratio * cell->c_data[5];
125 break;
126 case 10:
127 ratio = linear_interpolate(cell->s_data[3], cell->s_data[7], fvalue);
128 flip_ratio = 1 - ratio;
129 vert_geom->x = flip_ratio * cell->axis[9] + ratio * cell->axis[21];
130 vert_geom->y = flip_ratio * cell->axis[10] + ratio * cell->axis[22];
131 vert_geom->z = flip_ratio * cell->axis[11] + ratio * cell->axis[23];
132 *cdata = flip_ratio * cell->c_data[3] + ratio * cell->c_data[7];
133 break;
134 case 11:
135 ratio = linear_interpolate(cell->s_data[2], cell->s_data[6], fvalue);
136 flip_ratio = 1 - ratio;
137 vert_geom->x = flip_ratio * cell->axis[6] + ratio * cell->axis[18];
138 vert_geom->y = flip_ratio * cell->axis[7] + ratio * cell->axis[19];
139 vert_geom->z = flip_ratio * cell->axis[8] + ratio * cell->axis[20];
140 *cdata = flip_ratio * cell->c_data[2] + ratio * cell->c_data[6];
141 break;
142 }
143}
144
145/* return the geometry at the interpolate point in cube */
146void get_insidepoint(int inside_index, Cell *cell, double fvalue,
147 Fgeom *vert_geom, double *cdata, int disamb_flag) {
148 double ratio, flip_ratio;
149
150 switch (inside_index) {
151 case 0:
152 ratio = linear_interpolate(cell->s_data[0], cell->s_data[6], fvalue);
153 flip_ratio = 1 - ratio;
154 vert_geom->x = flip_ratio * cell->axis[0] + ratio * cell->axis[18];
155 vert_geom->y = flip_ratio * cell->axis[1] + ratio * cell->axis[19];
156 vert_geom->z = flip_ratio * cell->axis[2] + ratio * cell->axis[20];
157 *cdata = flip_ratio * cell->c_data[0] + ratio * cell->c_data[6];
158 break;
159 case 1:
160 ratio = linear_interpolate(cell->s_data[1], cell->s_data[7], fvalue);
161 flip_ratio = 1 - ratio;
162 vert_geom->x = flip_ratio * cell->axis[3] + ratio * cell->axis[21];
163 vert_geom->y = flip_ratio * cell->axis[4] + ratio * cell->axis[22];
164 vert_geom->z = flip_ratio * cell->axis[5] + ratio * cell->axis[23];
165 *cdata = flip_ratio * cell->c_data[1] + ratio * cell->c_data[7];
166 break;
167 case 2:
168 ratio = linear_interpolate(cell->s_data[2], cell->s_data[4], fvalue);
169 flip_ratio = 1 - ratio;
170 vert_geom->x = flip_ratio * cell->axis[6] + ratio * cell->axis[12];
171 vert_geom->y = flip_ratio * cell->axis[7] + ratio * cell->axis[13];
172 vert_geom->z = flip_ratio * cell->axis[8] + ratio * cell->axis[14];
173 *cdata = flip_ratio * cell->c_data[2] + ratio * cell->c_data[4];
174 break;
175 case 3:
176 ratio = linear_interpolate(cell->s_data[3], cell->s_data[5], fvalue);
177 flip_ratio = 1 - ratio;
178 vert_geom->x = flip_ratio * cell->axis[9] + ratio * cell->axis[15];
179 vert_geom->y = flip_ratio * cell->axis[10] + ratio * cell->axis[16];
180 vert_geom->z = flip_ratio * cell->axis[11] + ratio * cell->axis[17];
181 *cdata = flip_ratio * cell->c_data[3] + ratio * cell->c_data[5];
182 break;
183 case 4:
184 case 5:
185 case 6:
186 vert_geom->x =
187 (cell->axis[0] + cell->axis[3] + cell->axis[6] + cell->axis[9] +
188 cell->axis[12] + cell->axis[15] + cell->axis[18] + cell->axis[21]) /
189 8.0;
190 vert_geom->y =
191 (cell->axis[1] + cell->axis[4] + cell->axis[7] + cell->axis[10] +
192 cell->axis[13] + cell->axis[16] + cell->axis[19] + cell->axis[22]) /
193 8.0;
194 vert_geom->z =
195 (cell->axis[2] + cell->axis[5] + cell->axis[8] + cell->axis[11] +
196 cell->axis[14] + cell->axis[17] + cell->axis[20] + cell->axis[23]) /
197 8.0;
198 *cdata = (cell->c_data[0] + cell->c_data[1] + cell->c_data[2] +
199 cell->c_data[3] + cell->c_data[4] + cell->c_data[5] +
200 cell->c_data[6] + cell->c_data[7]) /
201 8.0;
202 break;
203 }
204}
205
206/* calculate linear interpolate value */
207double linear_interpolate(double left, double right, double fvalue) {
208 double ratio;
209
210 if (fabs(left - right) < EPSILON)
211 ratio = 0.0;
212 else
213 ratio = ((fvalue - left) / (right - left));
214 return ratio;
215}
216
217/* calculate the field value at the intersection point of the asymptotes */
218double calc_cross_field(double f00, double f10, double f11, double f01) {
219 double value;
220 if (fabs(f00 + f11 - f01 - f10) < EPSILON)
221 value = (f00 + f11 + f01 + f10) / 4.0;
222 else
223 value = ((f00 * f11 - f10 * f01) / (f00 + f11 - f01 - f10));
224 return value;
225}
226
227double facial_average(double f00, double f10, double f11, double f01) {
228 return ((f00 + f11 + f10 + f01) / 4);
229}
double calc_cross_field(double f00, double f10, double f11, double f01)
double linear_interpolate(double left, double right, double fvalue)
double facial_average(double f00, double f10, double f11, double f01)
void get_gridpoint(int voxel_index, Cell *cell, Fgeom *vert_geom, double *cdata)
void get_insidepoint(int inside_index, Cell *cell, double fvalue, Fgeom *vert_geom, double *cdata, int disamb_flag)
void get_edgepoint(int edge_index, Cell *cell, double fvalue, Fgeom *vert_geom, double *cdata)
int get_point_geom(int point_index, Cell *cell, double fvalue, Fgeom *point_geom, double *cdata, int disamb_flag)
#define EPSILON
double axis[3 *8]