FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 
6 #include "hecmw_vis_calc_attr.h"
7 
8 #include <math.h>
9 
10 /*----------------------------------------------------------------------
11 # Subroutines in this file on isosurface generation by Marching Cubes is
12 based
13  on the revision of Dr. Yuriko Takeshima's codes when she was working
14 part time in RIST
15 #---------------------------------------------------------------------- */
16 
17 /* calculate the geometry at the point in cube */
18 int 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 */
33 void 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 
41 void 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 */
146 void 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 */
207 double 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 */
218 double 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 
227 double facial_average(double f00, double f10, double f11, double f01) {
228  return ((f00 + f11 + f10 + f01) / 4);
229 }
facial_average
double facial_average(double f00, double f10, double f11, double f01)
Definition: hecmw_vis_calc_attr.c:227
_fgeom_struct::z
double z
Definition: hecmw_vis_SF_geom.h:182
calc_cross_field
double calc_cross_field(double f00, double f10, double f11, double f01)
Definition: hecmw_vis_calc_attr.c:218
_cell_struct::s_data
double s_data[8]
Definition: hecmw_vis_SF_geom.h:254
get_insidepoint
void get_insidepoint(int inside_index, Cell *cell, double fvalue, Fgeom *vert_geom, double *cdata, int disamb_flag)
Definition: hecmw_vis_calc_attr.c:146
_fgeom_struct
Definition: hecmw_vis_SF_geom.h:179
get_edgepoint
void get_edgepoint(int edge_index, Cell *cell, double fvalue, Fgeom *vert_geom, double *cdata)
Definition: hecmw_vis_calc_attr.c:41
_cell_struct::axis
double axis[3 *8]
Definition: hecmw_vis_SF_geom.h:253
_cell_struct::c_data
double c_data[8]
Definition: hecmw_vis_SF_geom.h:255
_cell_struct
Definition: hecmw_vis_SF_geom.h:252
linear_interpolate
double linear_interpolate(double left, double right, double fvalue)
Definition: hecmw_vis_calc_attr.c:207
hecmw_vis_calc_attr.h
get_gridpoint
void get_gridpoint(int voxel_index, Cell *cell, Fgeom *vert_geom, double *cdata)
Definition: hecmw_vis_calc_attr.c:33
get_point_geom
int get_point_geom(int point_index, Cell *cell, double fvalue, Fgeom *point_geom, double *cdata, int disamb_flag)
Definition: hecmw_vis_calc_attr.c:18
EPSILON
#define EPSILON
Definition: hecmw_vis_color_mapping.c:10
_fgeom_struct::y
double y
Definition: hecmw_vis_SF_geom.h:181
_fgeom_struct::x
double x
Definition: hecmw_vis_SF_geom.h:180