1 /*****************************************************************************
2  * Copyright (c) 2019 FrontISTR Commons
3  * This software is released under the MIT License, see LICENSE.txt
4  *****************************************************************************/
5 
6 /**
7  * @brief 並列計算された結果を読込み処理するためのユーティリティ
8  */
9 
10 #include "fstr_rmerge_util.h"
11 
12 extern FILE* log_fp;
13 extern int nrank;
14 extern int strid;
15 extern int endid;
16 extern int intid;
17 
out_log(const char * fmt,...)18 static void out_log(const char* fmt, ...) {
19   va_list arg;
20   va_start(arg, fmt);
21   vfprintf(log_fp, fmt, arg);
22   va_end(arg);
23 }
24 
25 /**
26  * @brief 全分散メッシュの読込み
27  */
28 
get_dist_fname(char * name_ID,char * fheader,int * fg_single,int * refine,int irank)29 static int get_dist_fname(char* name_ID, char* fheader, int* fg_single,
30                           int* refine, int irank) {
31   struct hecmw_ctrl_meshfiles* files;
32 
33   files = HECMW_ctrl_get_meshfiles_header_sub(name_ID, nrank, irank);
34   if (!files) return -1;
35 
36   if (files->n_mesh == 1) {
37     strcpy(fheader, files->meshfiles[0].filename);
38     if (files->meshfiles[0].type == HECMW_CTRL_FTYPE_HECMW_DIST) {
39       *fg_single = 0;
40     } else {
41       *fg_single = 1;
42     }
43     *refine = files->meshfiles[0].refine;
44     out_log("refine number is %d\n", *refine);
45   } else {
46     HECMW_ctrl_free_meshfiles(files);
47     return -1;
48   }
49 
50   HECMW_ctrl_free_meshfiles(files);
51   return 0;
52 }
53 
get_area_n(char * fheader)54 static int get_area_n(char* fheader) {
55   FILE* fp;
56   char buff[HECMW_FILENAME_LEN + 1];
57   int area = 0;
58 
59   while (1) {
60     sprintf(buff, "%s.%d", fheader, area);
61     out_log("try open : %s  ... ", buff);
62     fp = fopen(buff, "r");
63     if (!fp) {
64       out_log("fail\n");
65       out_log("area number is %d\n", area);
66       return area;
67     } else {
68       out_log("success\n");
69       fclose(fp);
70     }
71     area++;
72   }
73 }
74 
75 /**
76  * @brief 全分散メッシュの読込み
77  */
78 
fstr_get_all_local_mesh(char * name_ID,int * area_number,int * refine)79 struct hecmwST_local_mesh** fstr_get_all_local_mesh(char* name_ID,
80                                                     int* area_number,
81                                                     int* refine) {
82   int i;
83   int area_n;
84   int fg_single;
85   char fheader[HECMW_HEADER_LEN + 1];
86   char fname[HECMW_FILENAME_LEN + 1];
87   struct hecmwST_local_mesh** mesh;
88 
89   if (get_dist_fname(name_ID, fheader, &fg_single, refine, 0)) return NULL;
90 
91   if (fg_single) {
92     out_log("mesh file type is NOT HECMW_DIST.\n");
93     area_n = 1;
94     out_log("area number is %d\n", area_n);
95     mesh    = HECMW_malloc(area_n * sizeof(struct hecmwST_local_mesh*));
96     mesh[0] = HECMW_get_mesh(name_ID);
97     if (!mesh[0]) return NULL;
98   } else {
99     out_log("mesh file type is HECMW_DIST.\n");
100     if (nrank == 0) {
101       area_n = get_area_n(fheader);
102     } else {
103       area_n = nrank;
104     }
105     if (area_n == 0) return NULL;
106     mesh = HECMW_malloc(area_n * sizeof(struct hecmwST_local_mesh*));
107     for (i = 0; i < area_n; i++) {
108       if (nrank == 0) {
109         sprintf(fname, "%s.%d", fheader, i);
110       } else {
111         get_dist_fname(name_ID, fheader, &fg_single, refine, i);
112         sprintf(fname, "%s.%d", fheader, i);
113       }
114       out_log("loading dist mesh from %s\n", fname);
115       mesh[i] = HECMW_get_dist_mesh(fname);
116       if (!mesh[i]) return NULL;
117     }
118   }
119   *area_number = area_n;
120 
121   return mesh;
122 }
123 
124 /**
125  * @brief メッシュの削除
126  */
127 
fstr_free_mesh(struct hecmwST_local_mesh ** mesh,int area_n)128 void fstr_free_mesh(struct hecmwST_local_mesh** mesh, int area_n) {
129   int i;
130 
131   if (!mesh) return;
132 
133   for (i = 0; i < area_n; i++) {
134     HECMW_dist_free(mesh[i]);
135   }
136   HECMW_free(mesh);
137   return;
138 }
139 
140 /**
141  * @brief ステップ数を調べる(ファイルの存在を調べる)
142  */
143 
fstr_get_step_n(char * name_ID)144 int fstr_get_step_n(char* name_ID) {
145   FILE* fp;
146   int step, fg_text;
147   char* fheader;
148   char fname[HECMW_FILENAME_LEN + 1];
149 
150   if (endid > -1) return endid;
151 
152   if (nrank == 0) {
153     if ((fheader = HECMW_ctrl_get_result_fileheader(name_ID, 1, &fg_text)) ==
154         NULL)
155       return 0;
156   } else {
157     if ((fheader = HECMW_ctrl_get_result_fileheader_sub(name_ID, 1, nrank, 0,
158                                                         &fg_text)) == NULL)
159       return 0;
160   }
161 
162   step = 1;
163   while (1) {
164     sprintf(fname, "%s.0.%d", fheader, step);
165     out_log("try open : %s  ... ", fname);
166     fp = fopen(fname, "r");
167     if (!fp) {
168       out_log("fail\n");
169       out_log("step number is %d\n", step - 1);
170       return step - 1;
171     } else {
172       out_log("success\n");
173       fclose(fp);
174     }
175     step++;
176   }
177 }
178 
179 /**
180  * @brief ステップの全領域データの読み込み
181  */
182 
fstr_get_all_result(char * name_ID,int step,int area_n,int refine)183 fstr_res_info** fstr_get_all_result(char* name_ID, int step, int area_n,
184                                     int refine) {
185   char* fheader;
186   char fname[HECMW_FILENAME_LEN + 1];
187   fstr_res_info** res;
188   struct hecmwST_result_data* data;
189   int i, j, k, count, num, nnode, nelem, flag, fg_text;
190   int *node_gid, *elem_gid;
191   int refine_nnode = 0;
192 
193   if (nrank == 0) {
194     if ((fheader = HECMW_ctrl_get_result_fileheader(name_ID, step,
195                                                     &fg_text)) == NULL)
196       return 0;
197   }
198 
199   res = HECMW_malloc(area_n * sizeof(fstr_res_info*));
200   if (!res) return NULL;
201 
202   for (i = 0; i < area_n; i++) {
203     res[i] = HECMW_malloc(sizeof(fstr_res_info));
204     if (!res[i]) return NULL;
205 
206     if (nrank != 0) {
207       if ((fheader = HECMW_ctrl_get_result_fileheader_sub(
208                name_ID, step, nrank, i, &fg_text)) == NULL)
209         return 0;
210     }
211     sprintf(fname, "%s.%d.%d", fheader, i, step);
212     data = HECMW_result_read_by_fname(fname);
213     if (!data) return NULL;
214     nnode    = HECMW_result_get_nnode();
215     node_gid = HECMW_malloc(nnode * sizeof(int));
216     if (!node_gid) return NULL;
217     HECMW_result_get_nodeID(node_gid);
218     nelem    = HECMW_result_get_nelem();
219     elem_gid = HECMW_malloc(nelem * sizeof(int));
220     if (!elem_gid) return NULL;
221     if (data->ne_component) HECMW_result_get_elemID(elem_gid);
222 
223     if (refine) {
224       res[i]->result = HECMW_malloc(sizeof(struct hecmwST_result_data));
225       if (!res[i]->result) return NULL;
226 
227       res[i]->result->ng_component = data->ng_component;
228       res[i]->result->ng_dof = HECMW_malloc(data->ng_component * sizeof(int));
229       if (!res[i]->result->ng_dof) return NULL;
230       res[i]->result->global_label =
231           HECMW_malloc(data->ng_component * sizeof(char*));
232       if (!res[i]->result->global_label) return NULL;
233       for (j = 0; j < data->ng_component; j++) {
234         res[i]->result->ng_dof[j]     = data->ng_dof[j];
235         res[i]->result->global_label[j] = HECMW_strdup(data->global_label[j]);
236       }
237 
238       res[i]->result->nn_component = data->nn_component;
239       res[i]->result->nn_dof = HECMW_malloc(data->nn_component * sizeof(int));
240       if (!res[i]->result->nn_dof) return NULL;
241       res[i]->result->node_label =
242           HECMW_malloc(data->nn_component * sizeof(char*));
243       if (!res[i]->result->node_label) return NULL;
244       num = 0;
245       for (j = 0; j < data->nn_component; j++) {
246         res[i]->result->nn_dof[j]     = data->nn_dof[j];
247         res[i]->result->node_label[j] = HECMW_strdup(data->node_label[j]);
248         num += data->nn_dof[j];
249       }
250 
251       count = 1;
252       flag  = 0;
253       for (j = 1; j < nnode; j++) {
254         if (flag == refine) break;
255         count++;
256         if (node_gid[j] > 0 && node_gid[j - 1] < 0) flag++;
257         if (node_gid[j] < 0 && node_gid[j] > node_gid[j - 1]) flag++;
258       }
259       count--;
260       out_log("\narea:%d -- refined_nn_internal:%d", i, count);
261       refine_nnode += count;
262 
263       count = 0;
264       for (j = 0; j < nnode; j++)
265         if (node_gid[j] > 0) count++;
266       res[i]->nnode_gid = count;
267       res[i]->result->node_val_item =
268           HECMW_malloc(num * count * sizeof(double));
269       if (!res[i]->result->node_val_item) return NULL;
270       count = 0;
271       for (j = 0; j < nnode; j++) {
272         if (node_gid[j] > 0) {
273           for (k = 0; k < num; k++) {
274             res[i]->result->node_val_item[count++] =
275                 data->node_val_item[num * j + k];
276           }
277         }
278       }
279 
280       res[i]->result->ne_component = data->ne_component;
281       res[i]->result->ne_dof = HECMW_malloc(data->ne_component * sizeof(int));
282       if (!res[i]->result->ne_dof) return NULL;
283       res[i]->result->elem_label =
284           HECMW_malloc(data->ne_component * sizeof(char*));
285       if (!res[i]->result->elem_label) return NULL;
286       num = 0;
287       for (j = 0; j < data->ne_component; j++) {
288         res[i]->result->ne_dof[j]     = data->ne_dof[j];
289         res[i]->result->elem_label[j] = HECMW_strdup(data->elem_label[j]);
290         num += data->ne_dof[j];
291       }
292 
293       count = 0;
294       for (j = 0; j < nelem; j++)
295         if (elem_gid[j] > 0) count++;
296       out_log("\narea:%d -- ne_original from result:%d", i, count);
297       res[i]->nelem_gid = count;
298       res[i]->result->elem_val_item =
299           HECMW_malloc(num * count * sizeof(double));
300       if (!res[i]->result->elem_val_item) return NULL;
301       count = 0;
302       for (j = 0; j < nelem; j++) {
303         if (elem_gid[j] > 0) {
304           for (k = 0; k < num; k++) {
305             res[i]->result->elem_val_item[count++] =
306                 data->elem_val_item[num * j + k];
307           }
308         }
309       }
310 
311       HECMW_result_free(data);
312       HECMW_result_free_nodeID();
313       HECMW_result_free_elemID();
314 
315       res[i]->node_gid = HECMW_malloc(res[i]->nnode_gid * sizeof(int));
316       if (!res[i]->node_gid) return NULL;
317       count = 0;
318       for (j = 0; j < nnode; j++) {
319         if (node_gid[j] > 0) res[i]->node_gid[count++] = node_gid[j];
320       }
321       free(node_gid);
322 
323       res[i]->elem_gid = HECMW_malloc(res[i]->nelem_gid * sizeof(int));
324       if (!res[i]->elem_gid) return NULL;
325       count = 0;
326       for (j = 0; j < nelem; j++) {
327         if (elem_gid[j] > 0) res[i]->elem_gid[count++] = elem_gid[j];
328       }
329       free(elem_gid);
330     } else {
331       res[i]->result    = data;
332       res[i]->nnode_gid = nnode;
333       res[i]->node_gid  = node_gid;
334       res[i]->nelem_gid = nelem;
335       res[i]->elem_gid  = elem_gid;
336     }
337   }
338 
339   if (refine) out_log("\ntotal refined_nn_internal:%d\n", refine_nnode);
340 
341   return res;
342 }
343 
344 /**
345  * @brief ステップの全領域データの結合
346  */
347 
fstr_all_result(fstr_glt * glt,fstr_res_info ** res,int refine)348 struct hecmwST_result_data* fstr_all_result(fstr_glt* glt, fstr_res_info** res,
349                                             int refine) {
350   int i, j, k, l_id, area;
351   int gitem, nitem, eitem, count, irec;
352   struct hecmwST_result_data* data;
353 
354   gitem = 0;
355   for (i = 0; i < res[0]->result->ng_component; i++)
356     gitem += res[0]->result->ng_dof[i];
357   nitem = 0;
358   for (i = 0; i < res[0]->result->nn_component; i++)
359     nitem += res[0]->result->nn_dof[i];
360   eitem = 0;
361   for (i = 0; i < res[0]->result->ne_component; i++)
362     eitem += res[0]->result->ne_dof[i];
363 
364   data             = HECMW_malloc(sizeof(struct hecmwST_result_data));
365   data->ng_dof     = HECMW_malloc(res[0]->result->ng_component * sizeof(int));
366   data->global_label = HECMW_malloc(res[0]->result->ng_component * sizeof(char*));
367   data->global_val_item = HECMW_malloc(gitem * sizeof(double));
368   data->nn_dof     = HECMW_malloc(res[0]->result->nn_component * sizeof(int));
369   data->node_label = HECMW_malloc(res[0]->result->nn_component * sizeof(char*));
370   data->node_val_item = HECMW_malloc(nitem * glt->node_n * sizeof(double));
371   data->ne_dof     = HECMW_malloc(res[0]->result->ne_component * sizeof(int));
372   data->elem_label = HECMW_malloc(res[0]->result->ne_component * sizeof(char*));
373   data->elem_val_item = HECMW_malloc(eitem * glt->elem_n * sizeof(double));
374 
375   data->ng_component = res[0]->result->ng_component;
376   for (i = 0; i < res[0]->result->ng_component; i++) {
377     data->ng_dof[i]     = res[0]->result->ng_dof[i];
378     data->global_label[i] = HECMW_strdup(res[0]->result->global_label[i]);
379   }
380   for (i = 0; i < gitem; i++) {
381     data->global_val_item[i] = res[0]->result->global_val_item[i];
382   }
383   count = 0;
384   for (i = 0; i < glt->node_n; i++) {
385     l_id = glt->nrec[i].local;
386     area = glt->nrec[i].area;
387     irec = nitem * l_id;
388     for (j = 0; j < res[0]->result->nn_component; j++) {
389       for (k = 0; k < res[0]->result->nn_dof[j]; k++) {
390         data->node_val_item[count++] = res[area]->result->node_val_item[irec++];
391       }
392     }
393   }
394   data->nn_component = res[0]->result->nn_component;
395   for (i = 0; i < res[0]->result->nn_component; i++) {
396     data->nn_dof[i]     = res[0]->result->nn_dof[i];
397     data->node_label[i] = HECMW_strdup(res[0]->result->node_label[i]);
398   }
399 
400   count = 0;
401   for (i = 0; i < glt->elem_n; i++) {
402     l_id = glt->erec[i].local;
403     area = glt->erec[i].area;
404     irec = eitem * l_id;
405     for (j = 0; j < res[0]->result->ne_component; j++) {
406       for (k = 0; k < res[0]->result->ne_dof[j]; k++) {
407         if (refine) {
408           data->elem_val_item[count++] = 0.0;
409         } else {
410           data->elem_val_item[count++] =
411               res[area]->result->elem_val_item[irec++];
412         }
413       }
414     }
415   }
416   data->ne_component = res[0]->result->ne_component;
417   for (i = 0; i < res[0]->result->ne_component; i++) {
418     data->ne_dof[i]     = res[0]->result->ne_dof[i];
419     data->elem_label[i] = HECMW_strdup(res[0]->result->elem_label[i]);
420   }
421 
422   return data;
423 }
424 
425 /**
426  * @brief fstr_res_info の削除
427  */
428 
fstr_free_result(fstr_res_info ** res,int area_n)429 void fstr_free_result(fstr_res_info** res, int area_n) {
430   int i;
431 
432   if (!res) return;
433 
434   for (i = 0; i < area_n; i++) {
435     HECMW_result_free(res[i]->result);
436     HECMW_free(res[i]->node_gid);
437     HECMW_free(res[i]->elem_gid);
438     HECMW_free(res[i]);
439   }
440   HECMW_free(res);
441   return;
442 }
443 
444 /**
445  * @brief グローバルとローカル、所属領域のテーブル fstr_glt の作成
446  */
447 
cmp_global_glt(const fstr_gl_rec * g1,const fstr_gl_rec * g2)448 static int cmp_global_glt(const fstr_gl_rec* g1, const fstr_gl_rec* g2) {
449   return (g1->global - g2->global);
450 }
451 
452 typedef int (*cmp_func)(const void*, const void*);
453 
fstr_create_glt(struct hecmwST_local_mesh ** mesh,int area_n)454 fstr_glt* fstr_create_glt(struct hecmwST_local_mesh** mesh, int area_n) {
455   int i, j, k, eid, count;
456   int all_n, all_e;
457   int area_e;
458   fstr_gl_rec* nrec;
459   fstr_gl_rec* erec;
460   fstr_glt* glt;
461   int* area_etype_list;
462 
463   all_n = 0;
464   for (i = 0; i < area_n; i++) {
465     out_log("area:%d -- nn_internal:%d\n", i, mesh[i]->nn_internal);
466     all_n += mesh[i]->nn_internal;
467   }
468   out_log("total nn_internal:%d\n", all_n);
469 
470   nrec = HECMW_malloc(sizeof(fstr_gl_rec) * all_n);
471   if (!nrec) return NULL;
472 
473   count = 0;
474   for (i = 0; i < area_n; i++) {
475     for (j = 0; j < mesh[i]->nn_internal; j++) {
476       nrec[count].global = mesh[i]->global_node_ID[j];
477       nrec[count].local  = j;
478       nrec[count].area   = i;
479       count++;
480     }
481   }
482 
483   qsort(nrec, all_n, sizeof(fstr_gl_rec), (cmp_func)cmp_global_glt);
484 
485   all_e = 0;
486   for (i = 0; i < area_n; i++) {
487     area_e = 0;
488     area_etype_list = HECMW_malloc(sizeof(int) * mesh[i]->n_elem);
489     for (j = 0; j < mesh[i]->n_elem_type; j++) {
490       for (k = mesh[i]->elem_type_index[j]; k < mesh[i]->elem_type_index[j+1]; k++) {
491          area_etype_list[k] = mesh[i]->elem_type_item[j];
492       }
493     }
494 
495     for (j = 0; j < mesh[i]->ne_internal; j++) {
496       eid = mesh[i]->elem_internal_list[j] - 1;
497       if ( HECMW_is_etype_patch(area_etype_list[eid]) ) continue;
498       area_e++;
499     }
500 
501     HECMW_free(area_etype_list);
502     all_e += area_e;
503     out_log("area:%d -- ne_internal:%d\n", i, area_e);
504   }
505   out_log("total ne_internal:%d\n", all_e);
506 
507   erec = HECMW_malloc(sizeof(fstr_gl_rec) * all_e);
508   if (!erec) return NULL;
509 
510   count = 0;
511   for (i = 0; i < area_n; i++) {
512 
513     area_etype_list = HECMW_malloc(sizeof(int) * mesh[i]->n_elem);
514     for (j = 0; j < mesh[i]->n_elem_type; j++) {
515       for (k = mesh[i]->elem_type_index[j]; k < mesh[i]->elem_type_index[j+1]; k++) {
516          area_etype_list[k] = mesh[i]->elem_type_item[j];
517       }
518     }
519 
520     for (j = 0; j < mesh[i]->ne_internal; j++) {
521       eid = mesh[i]->elem_internal_list[j] - 1;
522       if ( HECMW_is_etype_patch(area_etype_list[eid]) ) continue;
523       if ( HECMW_is_etype_link(area_etype_list[eid]) ) continue;
524 
525       erec[count].global = mesh[i]->global_elem_ID[eid];
526       erec[count].local = eid;
527       erec[count].area  = i;
528       count++;
529     }
530 
531     HECMW_free(area_etype_list);
532   }
533 
534   qsort(erec, all_e, sizeof(fstr_gl_rec), (cmp_func)cmp_global_glt);
535 
536   glt = HECMW_malloc(sizeof(fstr_glt));
537   if (!glt) return NULL;
538   glt->nrec   = nrec;
539   glt->erec   = erec;
540   glt->node_n = all_n;
541   glt->elem_n = all_e;
542 
543   return glt;
544 }
545 
546 /**
547  * @brief fstr_glt の削除
548  */
549 
fstr_free_glt(fstr_glt * glt)550 void fstr_free_glt(fstr_glt* glt) {
551   if (!glt) return;
552 
553   HECMW_free(glt->nrec);
554   HECMW_free(glt->erec);
555   HECMW_free(glt);
556   return;
557 }
558 
559 /**
560  * @brief 単一領域メッシュの作成
561  */
562 
fstr_create_glmesh(fstr_glt * glt)563 struct hecmwST_local_mesh* fstr_create_glmesh(fstr_glt* glt) {
564   struct hecmwST_local_mesh* mesh;
565   int i;
566 
567   mesh                 = HECMW_calloc(1, sizeof(struct hecmwST_local_mesh));
568   mesh->global_node_ID = HECMW_malloc(glt->node_n * sizeof(int));
569   mesh->global_elem_ID = HECMW_malloc(glt->elem_n * sizeof(int));
570 
571   for (i = 0; i < glt->node_n; i++) {
572     mesh->global_node_ID[i] = glt->nrec[i].global;
573   }
574   mesh->n_node = glt->node_n;
575 
576   for (i = 0; i < glt->elem_n; i++) {
577     mesh->global_elem_ID[i] = glt->erec[i].global;
578   }
579   mesh->n_elem = glt->elem_n;
580 
581   return mesh;
582 }
583 
584 /**
585  * @brief 単一領域メッシュの削除
586  */
587 
fstr_free_glmesh(struct hecmwST_local_mesh * mesh)588 void fstr_free_glmesh(struct hecmwST_local_mesh* mesh) {
589   if (!mesh) return;
590 
591   HECMW_free(mesh->global_node_ID);
592   HECMW_free(mesh->global_elem_ID);
593   HECMW_free(mesh);
594   return;
595 }
596