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_generate_histogram_sf.h"
7 
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "hecmw_vis_bmp.h"
12 #include "hecmw_vis_font_texture.h"
13 #include "hecmw_vis_mem_util.h"
14 #include "hecmw_malloc.h"
15 
generate_histogram_graph_sf(struct surface_module * sf,int * color_list,struct hecmwST_result_data * data,double * mivalue,double * mavalue,Result * result,int mynode,int pesize,HECMW_Comm VIS_COMM,int color_system_type)16 void generate_histogram_graph_sf(struct surface_module *sf, int *color_list,
17                                  struct hecmwST_result_data *data,
18                                  double *mivalue, double *mavalue,
19                                  Result *result, int mynode, int pesize,
20                                  HECMW_Comm VIS_COMM, int color_system_type) {
21   int i, j, k, m, ii;
22   double delta, value, color[3];
23   int count[500], t_count[500], max_number, max_length, start_x, end_x, start_y,
24       end_y;
25   FILE *fp;
26   double *graph;
27   BITMAPFILEHEADER header; /* File header */
28 
29   unsigned char r, g, b;
30   int ri, gi, bi;
31   BITMAPINFOHEADER info;
32 
33   int start_xs, start_ys;
34   char buf[128];
35   int output7[7][7];
36   int color_id;
37 
38   for (i = 0; i < data->nn_component; i++) {
39     if (color_list[i] == 1) {
40       color_id = i;
41       delta    = (mavalue[i] - mivalue[i]) / 500.0;
42     }
43   }
44 
45   for (i = 0; i < 500; i++) {
46     count[i]   = 0;
47     t_count[i] = 0;
48   }
49 
50   for (ii = 1; ii <= sf[0].surface_style; ii++) {
51     if (sf[ii].display_method != 4) {
52       for (i = 0; i < result[ii - 1].n_vertex; i++) {
53         j = (int)((result[ii - 1].color[i] - mivalue[color_id]) / delta);
54         if (j < 0) j   = 0;
55         if (j > 499) j = 499;
56         count[j]++;
57       }
58     }
59   }
60   if (pesize > 1)
61     HECMW_Allreduce(count, t_count, 500, HECMW_INT, HECMW_SUM, VIS_COMM);
62   else {
63     for (i = 0; i < 500; i++) t_count[i] = count[i];
64   }
65   /*   fprintf(stderr, "count[2]=%d  t_count[2]=%d\n", count[2], t_count[2]);
66    */
67   if (mynode == 0) {
68     fp = fopen("histogram.bmp", "wb");
69     if (fp == NULL)
70       HECMW_vis_print_exit("Cannot generate the histogram output file");
71     graph = (double *)HECMW_calloc(400 * 530 * 3, sizeof(double));
72     if (graph == NULL) HECMW_vis_memory_exit("graph");
73     for (i = 0; i < 400 * 530 * 3; i++) graph[i] = 0.0;
74     max_number                                   = 0;
75     for (i = 0; i < 500; i++) {
76       if (t_count[i] > max_number) max_number = t_count[i];
77     }
78     if (max_number == 0)
79       HECMW_vis_print_exit(
80           "ERROR: HEC-MW-VIS-E2003:Cannot generate histogram graph, the number "
81           "of voxels is 0");
82     max_length = (int)(400 - 30 - 5 - 45 * 1.5);
83     start_x    = (int)(5 + 45 * 1.5 + 15);
84     start_y    = 15;
85     for (j = 0; j < 500; j++) {
86       end_x =
87           (int)((double)t_count[j] * (double)max_length / (double)max_number +
88                 start_x) +
89           2;
90       value = j / 500.0;
91       value2_to_rgb(value, color, color_system_type);
92       for (i = start_x; i < end_x; i++) {
93         graph[((j + 15) * 400 + i) * 3]     = color[0];
94         graph[((j + 15) * 400 + i) * 3 + 1] = color[1];
95         graph[((j + 15) * 400 + i) * 3 + 2] = color[2];
96       }
97     }
98     /*start mark scales */
99     start_y = 15;
100     end_y   = 515;
101     for (k = 0; k < 11; k++) {
102       value = mivalue[color_id] +
103               (mavalue[color_id] - mivalue[color_id]) / 10.0 * k;
104       start_ys = start_y + (int)((double)500.0 / 10 * k) - (int)7 / 2;
105       start_xs = 15;
106       sprintf(buf, "%3.2E", value);
107       for (m = 0; m < 9; m++) {
108         font7_generate(buf[8 - m], output7);
109         for (j = 0; j < 7; j++)
110           for (i = 0; i < 7; i++) {
111             graph[((start_ys + j) * 400 + start_xs - i) * 3] =
112                 (double)output7[6 - j][i];
113             graph[((start_ys + j) * 400 + start_xs - i) * 3 + 1] =
114                 (double)output7[6 - j][i];
115             graph[((start_ys + j) * 400 + start_xs - i) * 3 + 2] =
116                 (double)output7[6 - j][i];
117           }
118         start_xs += 7;
119         if ((value >= 0) && (m == 0)) start_xs -= 7;
120       }
121       if ((k != 0) && (k != 10)) {
122         start_ys += (int)7 / 2;
123 
124         for (i = start_x; i < start_x + 5; i++)
125           for (j = 0; j < 3; j++) graph[(start_ys * 400 + i) * 3 + j] = 1.0;
126       }
127     }
128     header.bfSize = 54 + 3 * 400 * 530;
129 #ifdef CONVERSE_ORDER
130     header.bfSize = change_unsigned_int_order(54 + 3 * 400 * 530);
131 #endif
132     header.bfReserved1 = 0;
133 #ifdef CONVERSE_ORDER
134     header.bfReserved1 = change_short_int_order(0);
135 #endif
136 
137     header.bfReserved2 = 0;
138 #ifdef CONVERSE_ORDER
139     header.bfReserved2 = change_short_int_order(0);
140 #endif
141 
142     header.bfOffBits = 54;
143 #ifdef CONVERSE_ORDER
144     header.bfOffBits = change_unsigned_int_order(54);
145 #endif
146 
147     info.biBitCount = 24;
148 #ifdef CONVERSE_ORDER
149     info.biBitCount = change_short_int_order(24);
150 #endif
151 
152     info.biSize = 40;
153 #ifdef CONVERSE_ORDER
154     info.biSize = change_unsigned_int_order(40);
155 #endif
156 
157     info.biWidth = 400;
158 #ifdef CONVERSE_ORDER
159     info.biWidth = change_int_order(400);
160 #endif
161 
162     info.biHeight = 530;
163 #ifdef CONVERSE_ORDER
164     info.biHeight = change_int_order(530);
165 #endif
166 
167     info.biSizeImage = 3 * 400 * 530;
168 #ifdef CONVERSE_ORDER
169     info.biSizeImage = change_unsigned_int_order(3 * 400 * 530);
170 #endif
171 
172     info.biClrImportant = 0;
173 #ifdef CONVERSE_ORDER
174     info.biClrImportant = change_unsigned_int_order(0);
175 #endif
176 
177     info.biClrUsed = 0;
178 #ifdef CONVERSE_ORDER
179     info.biClrUsed = change_unsigned_int_order(0);
180 #endif
181 
182     info.biCompression = 0;
183 #ifdef CONVERSE_ORDER
184     info.biCompression = change_unsigned_int_order(0);
185 #endif
186 
187     info.biPlanes = 1;
188 #ifdef CONVERSE_ORDER
189     info.biPlanes = change_short_int_order(1);
190 #endif
191 
192     info.biXPelsPerMeter = 3780;
193 #ifdef CONVERSE_ORDER
194     info.biXPelsPerMeter = change_int_order(3780);
195 #endif
196 
197     info.biYPelsPerMeter = 3780;
198 #ifdef CONVERSE_ORDER
199     info.biYPelsPerMeter = change_int_order(3780);
200 #endif
201 
202     putc('B', fp);
203     putc('M', fp);
204     fwrite(&(header.bfSize), sizeof(unsigned int), 1, fp);
205     fwrite(&header.bfReserved1, sizeof(unsigned short int), 1, fp);
206     fwrite(&header.bfReserved2, sizeof(unsigned short int), 1, fp);
207     fwrite(&header.bfOffBits, sizeof(unsigned int), 1, fp);
208     fwrite(&info, 40, 1, fp);
209     for (j = 0; j < 530; j++)
210       for (i = 400 - 1; i >= 0; i--) {
211         ri               = (int)(graph[(j * 400 + i) * 3] * 256);
212         gi               = (int)(graph[(j * 400 + i) * 3 + 1] * 256);
213         bi               = (int)(graph[(j * 400 + i) * 3 + 2] * 256);
214         if (ri < 0) ri   = 0;
215         if (ri > 255) ri = 255;
216         if (gi < 0) gi   = 0;
217         if (gi > 255) gi = 255;
218         if (bi < 0) bi   = 0;
219         if (bi > 255) bi = 255;
220         r                = ri;
221         g                = gi;
222         b                = bi;
223         putc(b, fp);
224         putc(g, fp);
225         putc(r, fp);
226       }
227 
228     fclose(fp);
229     HECMW_free(graph);
230   }
231   return;
232 }
233 
generate_interval_point_sf(struct surface_module * sf,int * color_list,struct hecmwST_result_data * data,double * mivalue,double * mavalue,Result * result,int mynode,int pesize,HECMW_Comm VIS_COMM,double * interval_point)234 void generate_interval_point_sf(struct surface_module *sf, int *color_list,
235                                 struct hecmwST_result_data *data,
236                                 double *mivalue, double *mavalue,
237                                 Result *result, int mynode, int pesize,
238                                 HECMW_Comm VIS_COMM, double *interval_point) {
239   int i, j, ii;
240   double delta;
241   int count[500], t_count[500], tmp_count[500], sum_count, interv_count, sum,
242       current_j;
243   int color_id;
244 
245   for (i = 0; i < data->nn_component; i++) {
246     if (color_list[i] == 1) {
247       color_id = i;
248       delta    = (mavalue[i] - mivalue[i]) / 500.0;
249     }
250   }
251 
252   for (i = 0; i < 500; i++) {
253     count[i]   = 0;
254     t_count[i] = 0;
255   }
256 
257   for (ii = 1; ii <= sf[0].surface_style; ii++) {
258     if (sf[ii].display_method != 4) {
259       for (i = 0; i < result[ii - 1].n_vertex; i++) {
260         j = (int)((result[ii - 1].color[i] - mivalue[color_id]) / delta);
261         if (j < 0) j   = 0;
262         if (j > 499) j = 499;
263         count[j]++;
264       }
265     }
266   }
267   if (pesize > 1)
268     HECMW_Allreduce(count, t_count, 500, HECMW_INT, HECMW_SUM, VIS_COMM);
269   else {
270     for (i = 0; i < 500; i++) t_count[i] = count[i];
271   }
272   sum_count = 0;
273   for (i = 0; i < 500; i++) {
274     sum_count += t_count[i];
275     tmp_count[i] = t_count[i];
276   }
277   interv_count      = (int)sum_count / 10;
278   current_j         = 0;
279   interval_point[0] = mivalue[color_id];
280   interval_point[1] = 0.0;
281   for (i = 1; i < 10; i++) {
282     interval_point[i * 2 + 1] = (double)i / 10.0;
283     sum                       = 0;
284     j                         = current_j;
285     while ((j < 500) && (sum < interv_count)) {
286       sum += t_count[j];
287       j++;
288     }
289     j--;
290     interval_point[i * 2] =
291         ((double)j / 500.0 +
292          1.0 / 500.0 *
293              (1.0 - (double)(sum - interv_count) / (double)tmp_count[j])) *
294             (mavalue[color_id] - mivalue[color_id]) +
295         mivalue[color_id];
296     t_count[j] = sum - interv_count;
297     current_j  = j;
298   }
299   interval_point[20] = mavalue[color_id];
300   interval_point[21] = 1.0;
301   fprintf(stderr, "The automatic color mapping set is :\n");
302   for (i = 0; i < 11; i++)
303     fprintf(stderr, "%lf    %lf   \n", interval_point[i * 2],
304             interval_point[i * 2 + 1]);
305   return;
306 }
307 
output_histogram_sf(struct surface_module * sf,int * color_list,struct hecmwST_result_data * data,double * mivalue,double * mavalue,Result * result,int mynode,int pesize,HECMW_Comm VIS_COMM)308 void output_histogram_sf(struct surface_module *sf, int *color_list,
309                          struct hecmwST_result_data *data, double *mivalue,
310                          double *mavalue, Result *result, int mynode,
311                          int pesize, HECMW_Comm VIS_COMM) {
312   int i, j, ii;
313   double delta;
314   int count[100], t_count[100];
315   FILE *fp;
316   int color_id;
317   for (i = 0; i < data->nn_component; i++) {
318     if (color_list[i] == 1) {
319       color_id = i;
320       delta    = (mavalue[i] - mivalue[i]) / 100.0;
321     }
322   }
323 
324   for (i = 0; i < 100; i++) {
325     count[i]   = 0;
326     t_count[i] = 0;
327   }
328 
329   for (ii = 1; ii <= sf[0].surface_style; ii++) {
330     if (sf[ii].display_method != 4) {
331       for (i = 0; i < result[ii - 1].n_vertex; i++) {
332         j = (int)((result[ii - 1].color[i] - mivalue[color_id]) / delta);
333         if (j < 0) j  = 0;
334         if (j > 99) j = 99;
335         count[j]++;
336       }
337     }
338   }
339   if (pesize > 1)
340     HECMW_Allreduce(count, t_count, 100, HECMW_INT, HECMW_SUM, VIS_COMM);
341   else {
342     for (i = 0; i < 100; i++) t_count[i] = count[i];
343   }
344 
345   if (mynode == 0) {
346     fp = fopen("histogram.file", "w");
347     if (fp == NULL)
348       HECMW_vis_print_exit("Cannot generate the histogram output file");
349     for (i = 0; i < 100; i++)
350       fprintf(fp, "%d   %d   -----(%lf --- %lf)\n", i, t_count[i],
351               mivalue[color_id] + i * delta,
352               mivalue[color_id] + (i + 1) * delta);
353     fclose(fp);
354   }
355   return;
356 }
357 
find_minmax_sf(struct hecmwST_local_mesh * mesh,int mynode,double range[6])358 void find_minmax_sf(struct hecmwST_local_mesh *mesh, int mynode,
359                     double range[6]) {
360   int i;
361   for (i = 0; i < 3; i++) {
362     range[i * 2]     = 1.0E17;
363     range[i * 2 + 1] = -1.0E17;
364   }
365   for (i = 0; i < mesh->n_node; i++) {
366     if (mesh->node[i * 3] < range[0]) range[0]     = mesh->node[i * 3];
367     if (mesh->node[i * 3] > range[1]) range[1]     = mesh->node[i * 3];
368     if (mesh->node[i * 3 + 1] < range[2]) range[2] = mesh->node[i * 3 + 1];
369     if (mesh->node[i * 3 + 1] > range[3]) range[3] = mesh->node[i * 3 + 1];
370     if (mesh->node[i * 3 + 2] < range[4]) range[4] = mesh->node[i * 3 + 2];
371     if (mesh->node[i * 3 + 2] > range[5]) range[5] = mesh->node[i * 3 + 2];
372   }
373 
374   return;
375 }
376 
find_patch_minmax_sf(Result * result,struct surface_module * sf,double range[6])377 void find_patch_minmax_sf(Result *result, struct surface_module *sf,
378                           double range[6]) {
379   int i, ii, j;
380 
381   for (i = 0; i < 3; i++) {
382     range[i * 2]     = 1.0E17;
383     range[i * 2 + 1] = -1.0E17;
384   }
385 
386   for (ii = 0; ii < sf[0].surface_style; ii++) {
387     for (i = 0; i < result[ii].n_vertex; i++) {
388       for (j = 0; j < 3; j++) {
389         if (result[ii].vertex[i * 3 + j] < range[j * 2])
390           range[j * 2] = result[ii].vertex[i * 3 + j];
391         if (result[ii].vertex[i * 3 + j] > range[j * 2 + 1])
392           range[j * 2 + 1] = result[ii].vertex[i * 3 + j];
393       }
394     }
395   }
396   return;
397 }
398 
find_new_patch_minmax_sf(Result * result,struct surface_module * sf,double range[6])399 void find_new_patch_minmax_sf(Result *result, struct surface_module *sf,
400                               double range[6]) {
401   int i, ii, j;
402   double new_v;
403 
404   for (ii = 0; ii < sf[0].surface_style; ii++) {
405     if (sf[ii + 1].deform_display_on == 1) {
406       for (i = 0; i < result[ii].n_vertex; i++) {
407         for (j = 0; j < 3; j++) {
408           new_v = result[ii].vertex[i * 3 + j] +
409                   result[ii].disp[i * 3 + j] * sf[ii + 1].disp_scale;
410           if (new_v < range[j * 2]) range[j * 2]         = new_v;
411           if (new_v > range[j * 2 + 1]) range[j * 2 + 1] = new_v;
412         }
413       }
414     }
415   }
416   return;
417 }
418 
get_deform_scale(struct surface_module * sf,int ii,double range_x,double range_y,double range_z,struct hecmwST_local_mesh * mesh,struct hecmwST_result_data * data,int pesize,HECMW_Comm VIS_COMM)419 void get_deform_scale(struct surface_module *sf, int ii, double range_x,
420                       double range_y, double range_z,
421                       struct hecmwST_local_mesh *mesh,
422                       struct hecmwST_result_data *data, int pesize,
423                       HECMW_Comm VIS_COMM) {
424   double max_disp, t_max_disp, tmp[3], disp, max_range, s_scale;
425   int i, j;
426   int tn_component, d_base;
427   int mynode;
428 
429   tn_component = 0;
430   for (i = 0; i < data->nn_component; i++) tn_component += data->nn_dof[i];
431   d_base = 0;
432   for (i = 0; i < sf[ii].disp_comp; i++) d_base += data->nn_dof[i];
433 
434   max_disp = 0.0;
435 
436   for (i = 0; i < mesh->n_node; i++) {
437     for (j   = 0; j < 3; j++)
438       tmp[j] = data->node_val_item[i * tn_component + d_base + j];
439     disp     = sqrt(tmp[0] * tmp[0] + tmp[1] * tmp[1] + tmp[2] * tmp[2]);
440     if (disp > max_disp) max_disp = disp;
441   }
442   HECMW_Comm_rank(VIS_COMM, &mynode);
443   if (pesize > 1)
444     HECMW_Allreduce(&max_disp, &t_max_disp, 1, HECMW_DOUBLE, HECMW_MAX,
445                     VIS_COMM);
446   else
447     t_max_disp = max_disp;
448 
449   max_range = sqrt(range_x * range_x + range_y * range_y + range_z * range_z);
450   if (fabs(t_max_disp) < EPSILON * 0.0001)
451     HECMW_vis_print_exit("No deformation in this dataset");
452   s_scale = max_range * 0.1 / t_max_disp;
453   if (sf[ii].disp_scale < 0.0)
454     sf[ii].disp_scale = s_scale;
455   else
456     sf[ii].disp_scale *= s_scale;
457   if (mynode == 0)
458     fprintf(stderr, "The real deformation scale = %lf\n", sf[ii].disp_scale);
459   return;
460 }
461