1 /*****************************************************************************
2  * Copyright (c) 2019 FrontISTR Commons
3  * This software is released under the MIT License, see LICENSE.txt
4  *****************************************************************************/
5 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include "hecmw_struct.h"
9 #include "hecmw_util.h"
10 #include "hecmw_common_define.h"
11 #include "hecmw_etype.h"
12 #include "hecmw_reorder.h"
13 
14 #define MASK_BIT(map, bit) ((map) |= (bit))
15 #define EVAL_BIT(map, bit) ((map) & (bit))
16 #define INV_BIT(map, bit) ((map) ^= (bit))
17 #define CLEAR_BIT(map, bit) \
18   ((map) |= (bit));         \
19   ((map) ^= (bit))
20 
21 #define BIT_DOF_TWO 1
22 #define BIT_DOF_THREE 2
23 #define BIT_DOF_SIX 4
24 #define BIT_DOF_FOUR 8
25 #define BIT_DOF_ALL (BIT_DOF_TWO | BIT_DOF_THREE | BIT_DOF_SIX | BIT_DOF_FOUR)
26 
27 #define HECMW_COMMON_EQUATION_BLOCK_NAME "EQUATION_BLOCK"
28 #define MY_RANK 1
29 #define NEIGHBOR_RANK 2
30 #define BOTH_RANK 3
31 #define MPC_BLOCK 4
32 #define CANDIDATE 8
33 #define ALL_BIT 255
34 
35 #ifdef DEBUG
36 #define dw_node_flag(ptr, nmemb) dw_node_flag_(ptr, nmemb, __FILE__, __LINE__)
37 #else
38 #define dw_node_flag(ptr, nmemb) ((void)0)
39 #endif
40 
41 struct equation_block {
42   int n_eqn_block;
43   int n_mpc_block;
44   int *eqn_block_index;
45 };
46 
47 /*----------------------------------------------------------------------------*/
48 #if 0
49 static int
50 get_eqn_block_idx( struct hecmwST_local_mesh *local_mesh )
51 {
52   int i;
53 
54   for( i=0; i<local_mesh->node_group->n_grp; i++ )
55     if( !strcmp( local_mesh->node_group->grp_name[i], HECMW_COMMON_EQUATION_BLOCK_NAME ) ) return i;
56 
57   return -1;
58 }
59 #endif
60 
61 /*                                                                            */
62 /*  convert node id from old to new                                           */
63 /*                                                                            */
64 
65 /*----------------------------------------------------------------------------*/
66 /*  nodal coordinates < node >                                                */
67 /*----------------------------------------------------------------------------*/
old2new_node(struct hecmwST_local_mesh * local_mesh,int * node_new2old)68 static int old2new_node(struct hecmwST_local_mesh *local_mesh,
69                         int *node_new2old) {
70   double *new, *old;
71   int i;
72 
73   new = (double *)HECMW_malloc(sizeof(double) * local_mesh->n_node * 3);
74   if (new == NULL) {
75     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
76     return -1;
77   }
78 
79   old = local_mesh->node;
80   HECMW_assert(old);
81 
82   for (i = 0; i < local_mesh->n_node; i++) {
83     new[3 * i]     = old[3 * (node_new2old[i] - 1)];
84     new[3 * i + 1] = old[3 * (node_new2old[i] - 1) + 1];
85     new[3 * i + 2] = old[3 * (node_new2old[i] - 1) + 2];
86   }
87 
88   local_mesh->node = new;
89 
90   HECMW_free(old);
91 
92   return 0;
93 }
94 
95 /*----------------------------------------------------------------------------*/
96 /*  local node id & belonging domain of node < node_ID >                      */
97 /*----------------------------------------------------------------------------*/
old2new_node_ID(struct hecmwST_local_mesh * local_mesh,int * node_new2old)98 static int old2new_node_ID(struct hecmwST_local_mesh *local_mesh,
99                            int *node_new2old) {
100   int *new, *old;
101   int i;
102 
103   new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node * 2);
104   if (new == NULL) {
105     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
106     return -1;
107   }
108 
109   old = local_mesh->node_ID;
110   HECMW_assert(old);
111 
112   for (i = 0; i < local_mesh->n_node; i++) {
113     new[2 * i]     = old[2 * (node_new2old[i] - 1)];
114     new[2 * i + 1] = old[2 * (node_new2old[i] - 1) + 1];
115   }
116 
117   local_mesh->node_ID = new;
118 
119   HECMW_free(old);
120 
121   return 0;
122 }
123 
124 /*----------------------------------------------------------------------------*/
125 /*  global node id                                                            */
126 /*----------------------------------------------------------------------------*/
old2new_global_node_ID(struct hecmwST_local_mesh * local_mesh,int * node_new2old)127 static int old2new_global_node_ID(struct hecmwST_local_mesh *local_mesh,
128                                   int *node_new2old) {
129   int *new, *old;
130   int i;
131 
132   new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node);
133   if (new == NULL) {
134     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
135     return -1;
136   }
137 
138   old = local_mesh->global_node_ID;
139   HECMW_assert(old);
140 
141   for (i = 0; i < local_mesh->n_node; i++) {
142     new[i] = old[node_new2old[i] - 1];
143   }
144 
145   local_mesh->global_node_ID = new;
146 
147   HECMW_free(old);
148 
149   return 0;
150 }
151 
152 /*----------------------------------------------------------------------------*/
153 /*  initial conditions of node                                                */
154 /*----------------------------------------------------------------------------*/
old2new_node_init_val(struct hecmwST_local_mesh * local_mesh,int * node_new2old)155 static int old2new_node_init_val(struct hecmwST_local_mesh *local_mesh,
156                                  int *node_new2old) {
157   int *new_index, *old_index;
158   double *new_item, *old_item;
159   int old_id;
160   int counter;
161   int i, j;
162 
163   new_index = (int *)HECMW_calloc(local_mesh->n_node + 1, sizeof(int));
164   if (new_index == NULL) {
165     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
166     return -1;
167   }
168   new_item = (double *)HECMW_malloc(
169       sizeof(double) * local_mesh->node_init_val_index[local_mesh->n_node]);
170   if (new_item == NULL) {
171     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
172     return -1;
173   }
174 
175   old_index = local_mesh->node_init_val_index;
176   old_item  = local_mesh->node_init_val_item;
177   HECMW_assert(old_index);
178   HECMW_assert(old_item);
179 
180   for (counter = 0, i = 0; i < local_mesh->n_node; i++) {
181     old_id = node_new2old[i];
182 
183     for (j = old_index[old_id - 1]; j < old_index[old_id]; j++) {
184       new_item[counter++] = old_item[j];
185     }
186     new_index[i + 1] = counter;
187   }
188 
189   local_mesh->node_init_val_index = new_index;
190   local_mesh->node_init_val_item  = new_item;
191 
192   HECMW_free(old_item);
193   HECMW_free(old_index);
194 
195   return 0;
196 }
197 
198 /*----------------------------------------------------------------------------*/
199 /*  component nodes of element < elem_node_item >                             */
200 /*----------------------------------------------------------------------------*/
old2new_elem_node_item(struct hecmwST_local_mesh * local_mesh,int * node_old2new)201 static int old2new_elem_node_item(struct hecmwST_local_mesh *local_mesh,
202                                   int *node_old2new) {
203   int i;
204 
205   for (i = 0; i < local_mesh->elem_node_index[local_mesh->n_elem]; i++) {
206     local_mesh->elem_node_item[i] =
207         node_old2new[local_mesh->elem_node_item[i] - 1];
208   }
209 
210   return 0;
211 }
212 
213 /*----------------------------------------------------------------------------*/
214 /*  MPC group < local_mesh->mpc >                                             */
215 /*----------------------------------------------------------------------------*/
old2new_mpc_item(struct hecmwST_local_mesh * local_mesh,int * node_old2new)216 static int old2new_mpc_item(struct hecmwST_local_mesh *local_mesh,
217                             int *node_old2new) {
218   int i;
219 
220   if (!local_mesh->mpc->n_mpc) return 0;
221 
222   for (i = 0; i < local_mesh->mpc->mpc_index[local_mesh->mpc->n_mpc]; i++) {
223     local_mesh->mpc->mpc_item[i] =
224         node_old2new[local_mesh->mpc->mpc_item[i] - 1];
225   }
226 
227   return 0;
228 }
229 
230 /*----------------------------------------------------------------------------*/
231 /*  node group < local_mesh->node_group->grp_item >                           */
232 /*----------------------------------------------------------------------------*/
old2new_node_grp_item(struct hecmwST_local_mesh * local_mesh,int * node_old2new)233 static int old2new_node_grp_item(struct hecmwST_local_mesh *local_mesh,
234                                  int *node_old2new) {
235   int i, j;
236 
237   if (!local_mesh->node_group->n_grp) return 0;
238 
239   for (i = 0; i < local_mesh->node_group->n_grp; i++) {
240     if (strcmp(local_mesh->node_group->grp_name[i],
241                HECMW_COMMON_EQUATION_BLOCK_NAME)) {
242       for (j = local_mesh->node_group->grp_index[i];
243            j < local_mesh->node_group->grp_index[i + 1]; j++) {
244         local_mesh->node_group->grp_item[j] =
245             node_old2new[local_mesh->node_group->grp_item[j] - 1];
246       }
247     }
248   }
249 
250   return 0;
251 }
252 
253 /*============================================================================*/
254 /*  convert node id from old to new                                           */
255 /*============================================================================*/
old2new_node_info(struct hecmwST_local_mesh * local_mesh,int * node_new2old,int * node_old2new)256 static int old2new_node_info(struct hecmwST_local_mesh *local_mesh,
257                              int *node_new2old, int *node_old2new) {
258   /* nodal coordinates */
259   if (old2new_node(local_mesh, node_new2old)) return -1;
260 
261   /* local node id & belonging domain of node */
262   if (old2new_node_ID(local_mesh, node_new2old)) return -1;
263 
264   /* global node id */
265   if (old2new_global_node_ID(local_mesh, node_new2old)) return -1;
266 
267   /* initial conditions of node */
268   if (local_mesh->hecmw_flag_initcon) {
269     if (old2new_node_init_val(local_mesh, node_new2old)) return -1;
270   }
271 
272   /* component nodes of element */
273   if (old2new_elem_node_item(local_mesh, node_old2new)) return -1;
274 
275   /* MPC group */
276   if (local_mesh->mpc->n_mpc) {
277     if (old2new_mpc_item(local_mesh, node_old2new)) return -1;
278   }
279 
280   /* node group */
281   if (local_mesh->node_group->n_grp) {
282     if (old2new_node_grp_item(local_mesh, node_old2new)) return -1;
283   }
284 
285   return 0;
286 }
287 
288 /*                                                                            */
289 /*  convert element id from old to new                                        */
290 /*                                                                            */
291 
292 /*----------------------------------------------------------------------------*/
293 /*  finite element type < elem_type >                                         */
294 /*----------------------------------------------------------------------------*/
old2new_elem_type(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)295 static int old2new_elem_type(struct hecmwST_local_mesh *local_mesh,
296                              int *elem_new2old) {
297   int *new, *old;
298   int i;
299 
300   new = HECMW_malloc(sizeof(int) * local_mesh->n_elem);
301   if (new == NULL) {
302     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
303     return -1;
304   }
305 
306   old = local_mesh->elem_type;
307   HECMW_assert(old);
308 
309   for (i = 0; i < local_mesh->n_elem; i++) {
310     new[i] = local_mesh->elem_type[elem_new2old[i] - 1];
311   }
312 
313   local_mesh->elem_type = new;
314 
315   HECMW_free(old);
316 
317   return 0;
318 }
319 
320 /*----------------------------------------------------------------------------*/
321 /*  component nodes of element < elem_node_index, elem_node_item >            */
322 /*----------------------------------------------------------------------------*/
old2new_elem_node(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)323 static int old2new_elem_node(struct hecmwST_local_mesh *local_mesh,
324                              int *elem_new2old) {
325   int *new_index, *old_index;
326   int *new_item, *old_item;
327   int old_id;
328   int counter;
329   int i, j;
330 
331   new_index = HECMW_calloc(local_mesh->n_elem + 1, sizeof(int));
332   if (new_index == NULL) {
333     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
334     return -1;
335   }
336   new_item = HECMW_malloc(sizeof(int) *
337                           local_mesh->elem_node_index[local_mesh->n_elem]);
338   if (new_item == NULL) {
339     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
340     return -1;
341   }
342 
343   old_index = local_mesh->elem_node_index;
344   old_item  = local_mesh->elem_node_item;
345   HECMW_assert(old_index);
346   HECMW_assert(old_item);
347 
348   for (counter = 0, i = 0; i < local_mesh->n_elem; i++) {
349     old_id = elem_new2old[i];
350 
351     for (j = old_index[old_id - 1]; j < old_index[old_id]; j++) {
352       new_item[counter++] = old_item[j];
353     }
354     new_index[i + 1] = counter;
355   }
356   HECMW_assert(counter > local_mesh->n_elem);
357 
358   local_mesh->elem_node_index = new_index;
359   local_mesh->elem_node_item  = new_item;
360 
361   HECMW_free(old_item);
362   HECMW_free(old_index);
363 
364   return 0;
365 }
366 
367 /*----------------------------------------------------------------------------*/
368 /*  local element ID & belonging domain of element < elem_ID >                */
369 /*----------------------------------------------------------------------------*/
old2new_elem_ID(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)370 static int old2new_elem_ID(struct hecmwST_local_mesh *local_mesh,
371                            int *elem_new2old) {
372   int *new, *old;
373   int i;
374 
375   new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_elem * 2);
376   if (new == NULL) {
377     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
378     return -1;
379   }
380 
381   old = local_mesh->elem_ID;
382   HECMW_assert(old);
383 
384   for (i = 0; i < local_mesh->n_elem; i++) {
385     new[2 * i]     = old[2 * (elem_new2old[i] - 1)];
386     new[2 * i + 1] = old[2 * (elem_new2old[i] - 1) + 1];
387   }
388 
389   local_mesh->elem_ID = new;
390 
391   HECMW_free(old);
392 
393   return 0;
394 }
395 
396 /*----------------------------------------------------------------------------*/
397 /*  global element id < global_elem_ID >                                      */
398 /*----------------------------------------------------------------------------*/
old2new_global_elem_ID(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)399 static int old2new_global_elem_ID(struct hecmwST_local_mesh *local_mesh,
400                                   int *elem_new2old) {
401   int *new, *old;
402   int i;
403 
404   new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_elem);
405   if (new == NULL) {
406     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
407     return -1;
408   }
409 
410   old = local_mesh->global_elem_ID;
411   HECMW_assert(old);
412 
413   for (i = 0; i < local_mesh->n_elem; i++) {
414     new[i] = old[elem_new2old[i] - 1];
415   }
416 
417   local_mesh->global_elem_ID = new;
418 
419   HECMW_free(old);
420 
421   return 0;
422 }
423 
424 /*----------------------------------------------------------------------------*/
425 /*  section id < section_ID >                                                 */
426 /*----------------------------------------------------------------------------*/
old2new_section_ID(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)427 static int old2new_section_ID(struct hecmwST_local_mesh *local_mesh,
428                               int *elem_new2old) {
429   int *new, *old;
430   int i;
431 
432   new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_elem);
433   if (new == NULL) {
434     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
435     return -1;
436   }
437 
438   old = local_mesh->section_ID;
439   HECMW_assert(old);
440 
441   for (i = 0; i < local_mesh->n_elem; i++) {
442     new[i] = local_mesh->section_ID[elem_new2old[i] - 1];
443   }
444 
445   local_mesh->section_ID = new;
446 
447   HECMW_free(old);
448 
449   return 0;
450 }
451 
452 /*----------------------------------------------------------------------------*/
453 /*  material id < elem_mat_ID_index, elem_mat_ID_item >                       */
454 /*----------------------------------------------------------------------------*/
old2new_mat_ID(struct hecmwST_local_mesh * local_mesh,int * elem_new2old)455 static int old2new_mat_ID(struct hecmwST_local_mesh *local_mesh,
456                           int *elem_new2old) {
457   int *new_index, *old_index;
458   int *new_item, *old_item;
459   int old_id;
460   int counter;
461   int i, j;
462 
463   new_index = (int *)HECMW_calloc(local_mesh->n_elem + 1, sizeof(int));
464   if (new_index == NULL) {
465     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
466     return -1;
467   }
468   new_item = (int *)HECMW_malloc(
469       sizeof(int) * local_mesh->elem_mat_ID_index[local_mesh->n_elem]);
470   if (new_item == NULL) {
471     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
472     return -1;
473   }
474 
475   old_index = local_mesh->elem_mat_ID_index;
476   old_item  = local_mesh->elem_mat_ID_item;
477   HECMW_assert(old_index);
478   HECMW_assert(old_item);
479 
480   for (counter = 0, i = 0; i < local_mesh->n_elem; i++) {
481     old_id = elem_new2old[i];
482 
483     for (j = old_index[old_id - 1]; j < old_index[old_id]; j++) {
484       new_item[counter++] = old_item[j];
485     }
486     new_index[i + 1] = counter;
487   }
488   HECMW_assert(counter >= local_mesh->n_elem);
489 
490   local_mesh->elem_mat_ID_index = new_index;
491   local_mesh->elem_mat_ID_item  = new_item;
492 
493   HECMW_free(old_item);
494   HECMW_free(old_index);
495 
496   return 0;
497 }
498 
499 /*----------------------------------------------------------------------------*/
500 /*  list of internal element < elem_internal_list >                           */
501 /*----------------------------------------------------------------------------*/
old2new_elem_internal_list(struct hecmwST_local_mesh * local_mesh,int * elem_old2new)502 static int old2new_elem_internal_list(struct hecmwST_local_mesh *local_mesh,
503                                       int *elem_old2new) {
504   int i;
505 
506   for (i = 0; i < local_mesh->ne_internal; i++) {
507     local_mesh->elem_internal_list[i] =
508         elem_old2new[local_mesh->elem_internal_list[i] - 1];
509   }
510 
511   return 0;
512 }
513 
514 /*----------------------------------------------------------------------------*/
515 /*  element group < elem_group->grp_item >                                    */
516 /*----------------------------------------------------------------------------*/
old2new_elem_grp_item(struct hecmwST_local_mesh * local_mesh,int * elem_old2new)517 static int old2new_elem_grp_item(struct hecmwST_local_mesh *local_mesh,
518                                  int *elem_old2new) {
519   int i;
520 
521   for (i = 0;
522        i < local_mesh->elem_group->grp_index[local_mesh->elem_group->n_grp];
523        i++) {
524     local_mesh->elem_group->grp_item[i] =
525         elem_old2new[local_mesh->elem_group->grp_item[i] - 1];
526   }
527 
528   return 0;
529 }
530 
531 /*----------------------------------------------------------------------------*/
532 /*  surface group < surf_group->grp_item >                                    */
533 /*----------------------------------------------------------------------------*/
old2new_surf_grp_item(struct hecmwST_local_mesh * local_mesh,int * elem_old2new)534 static int old2new_surf_grp_item(struct hecmwST_local_mesh *local_mesh,
535                                  int *elem_old2new) {
536   int i;
537 
538   for (i = 0;
539        i < local_mesh->surf_group->grp_index[local_mesh->surf_group->n_grp];
540        i++) {
541     local_mesh->surf_group->grp_item[2 * i] =
542         elem_old2new[local_mesh->surf_group->grp_item[2 * i] - 1];
543   }
544 
545   return 0;
546 }
547 
548 /*============================================================================*/
549 /*  convert element id from old to new                                        */
550 /*============================================================================*/
old2new_elem_info(struct hecmwST_local_mesh * local_mesh,int * elem_new2old,int * elem_old2new)551 static int old2new_elem_info(struct hecmwST_local_mesh *local_mesh,
552                              int *elem_new2old, int *elem_old2new) {
553   /* finite element type */
554   if (old2new_elem_type(local_mesh, elem_new2old)) return -1;
555 
556   /* component nodes of element */
557   if (old2new_elem_node(local_mesh, elem_new2old)) return -1;
558 
559   /* local element id & belonging domain of element */
560   if (old2new_elem_ID(local_mesh, elem_new2old)) return -1;
561 
562   /* global element id */
563   if (old2new_global_elem_ID(local_mesh, elem_new2old)) return -1;
564 
565   /* section id */
566   if (old2new_section_ID(local_mesh, elem_new2old)) return -1;
567 
568   /* material id */
569   if (old2new_mat_ID(local_mesh, elem_new2old)) return -1;
570 
571   /* list of internal element */
572   if (old2new_elem_internal_list(local_mesh, elem_old2new)) return -1;
573 
574   /* element group */
575   if (local_mesh->elem_group->n_grp) {
576     if (old2new_elem_grp_item(local_mesh, elem_old2new)) return -1;
577   }
578 
579   /* surface group */
580   if (local_mesh->surf_group->n_grp) {
581     if (old2new_surf_grp_item(local_mesh, elem_old2new)) return -1;
582   }
583 
584   return 0;
585 }
586 
587 /*                                                                            */
588 /*  reorder member of EQUATION_BLOCK                                          */
589 /*                                                                            */
590 #if 0
591 /*----------------------------------------------------------------------------*/
592 /*  node group ( only EQUATION_BLOCK ) < node_group->grp_item >               */
593 /*----------------------------------------------------------------------------*/
594 static int
595 old2new_eqn_block( struct hecmwST_local_mesh *local_mesh, int *eqn_block_old2new )
596 {
597   struct hecmwST_node_grp *grp=local_mesh->node_group;
598   int n_eqn_block, eqn_block_idx;
599   int *new_item;
600   int counter;
601   int i, js;
602 
603   /* index of EQUATION_BLOCK */
604   eqn_block_idx = get_eqn_block_idx( local_mesh );
605   HECMW_assert( eqn_block_idx >= 0 );
606 
607   /* number of EQUATION_BLOCKs */
608   n_eqn_block = grp->grp_index[eqn_block_idx+1] - grp->grp_index[eqn_block_idx];
609   HECMW_assert( n_eqn_block > 0 );
610 
611   /* order EQUATION_BLOCK */
612   new_item = (int *)HECMW_calloc( n_eqn_block, sizeof(int) );
613   if( new_item == NULL ) {
614     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
615     return -1;
616   }
617 
618   for( js=0, i=grp->grp_index[eqn_block_idx]; i<grp->grp_index[eqn_block_idx+1]; i++ ) {
619     new_item[eqn_block_old2new[i-grp->grp_index[eqn_block_idx]]] = grp->grp_item[i] - js;
620     js = grp->grp_item[i];
621   }
622 
623   for( counter=0, i=0; i<n_eqn_block; i++ ) {
624     counter += new_item[i];
625     grp->grp_item[grp->grp_index[eqn_block_idx]+i] = counter;
626   }
627 
628   HECMW_free( new_item );
629 
630   return 0;
631 }
632 
633 /*============================================================================*/
634 /*  reorder member of EQUATION_BLOCK                                          */
635 /*============================================================================*/
636 static int
637 old2new_eqn_block_info( struct hecmwST_local_mesh *local_mesh, int *eqn_block_old2new )
638 {
639   /* node group ( only EQUATION_BLOCK ) */
640   if(old2new_eqn_block( local_mesh, eqn_block_old2new )) return -1;
641 
642   return 0;
643 }
644 #endif
645 /*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$*/
646 
647 /*                                                                            */
648 /*  reorder element according to finite element type                          */
649 /*                                                                            */
650 
651 /*----------------------------------------------------------------------------*/
652 /*  count elements of each finite element type                                */
653 /*----------------------------------------------------------------------------*/
count_each_elem_type(struct hecmwST_local_mesh * local_mesh,int * counter)654 static int count_each_elem_type(struct hecmwST_local_mesh *local_mesh,
655                                 int *counter) {
656   int etype;
657   int i;
658 
659   /* count elements of each finite element type */
660   for (i = 0; i < local_mesh->n_elem; i++) {
661     etype = HECMW_get_etype_HECMW2UTIL(local_mesh->elem_type[i]);
662     counter[etype]++;
663   }
664 
665   return 0;
666 }
667 
668 /*----------------------------------------------------------------------------*/
669 /*  count element types in mesh                                               */
670 /*----------------------------------------------------------------------------*/
set_n_elem_type(struct hecmwST_local_mesh * local_mesh,int * counter)671 static int set_n_elem_type(struct hecmwST_local_mesh *local_mesh,
672                            int *counter) {
673   int i;
674 
675   /* count element types */
676   local_mesh->n_elem_type = 0;
677   for (i = 0; i < HECMW_MESH_ETYPE_MAX + 1; i++) {
678     if (counter[i]) {
679       (local_mesh->n_elem_type)++;
680     }
681   }
682 
683   /* check data */
684   if (local_mesh->n_elem_type <= 0) {
685     HECMW_set_error(HECMW_COMMON_E_OUT_OF_RANGE, "");
686     return -1;
687   }
688 
689   return 0;
690 }
691 
692 /*----------------------------------------------------------------------------*/
693 /*  create item and index for finite element type & set conversion table      */
694 /*  between old and new element id                                            */
695 /*----------------------------------------------------------------------------*/
set_elem_type_index(struct hecmwST_local_mesh * local_mesh,int * counter,int * elem_new2old,int * elem_old2new)696 static int set_elem_type_index(struct hecmwST_local_mesh *local_mesh,
697                                int *counter, int *elem_new2old,
698                                int *elem_old2new) {
699   int etype;
700   int types, elems;
701   int i, j;
702 
703   /* allocation */
704   if (local_mesh->elem_type_index) {
705     HECMW_free(local_mesh->elem_type_index);
706   }
707   local_mesh->elem_type_index =
708       (int *)HECMW_calloc(local_mesh->n_elem_type + 1, sizeof(int));
709   if (local_mesh->elem_type_index == NULL) {
710     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
711     return -1;
712   }
713 
714   if (local_mesh->elem_type_item) {
715     HECMW_free(local_mesh->elem_type_item);
716   }
717   local_mesh->elem_type_item =
718       (int *)HECMW_malloc(sizeof(int) * local_mesh->n_elem_type);
719   if (local_mesh->elem_type_item == NULL) {
720     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
721     return -1;
722   }
723 
724   /* create index for element type */
725   for (elems = 0, types = 0, i = 0; i < HECMW_MESH_ETYPE_MAX + 1; i++) {
726     if (counter[i]) {
727       etype = HECMW_get_etype_UTIL2HECMW(i);
728 
729       for (j = 0; j < local_mesh->n_elem; j++) {
730         if (local_mesh->elem_type[j] == etype) {
731           elem_new2old[elems] = j + 1;
732           elem_old2new[j]     = elems + 1;
733           elems++;
734         }
735       }
736       local_mesh->elem_type_index[types + 1] = elems;
737       local_mesh->elem_type_item[types]      = etype;
738       types++;
739     }
740   }
741 
742   return 0;
743 }
744 
745 /*============================================================================*/
746 /*  reoreder element accordint to finite element type                         */
747 /*============================================================================*/
HECMW_reorder_elem_type(struct hecmwST_local_mesh * local_mesh)748 extern int HECMW_reorder_elem_type(struct hecmwST_local_mesh *local_mesh) {
749   int *counter; /* counter of elements */
750   int *elem_new2old,
751       *elem_old2new; /* conversion table between old and new element id */
752 
753   HECMW_assert(local_mesh);
754 
755   /* allocation */
756   counter = (int *)HECMW_calloc(HECMW_MESH_ETYPE_MAX + 1, sizeof(int));
757   if (counter == NULL) {
758     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
759     return -1;
760   }
761 
762   elem_new2old = (int *)HECMW_calloc(local_mesh->n_elem, sizeof(int));
763   if (elem_new2old == NULL) {
764     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
765     return -1;
766   }
767 
768   elem_old2new = (int *)HECMW_calloc(local_mesh->n_elem, sizeof(int));
769   if (elem_old2new == NULL) {
770     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
771     return -1;
772   }
773 
774   /* count elements of each element type */
775   if (count_each_elem_type(local_mesh, counter)) {
776     return -1;
777   }
778 
779   /* count finite element types in data */
780   if (set_n_elem_type(local_mesh, counter)) {
781     return -1;
782   }
783 
784   /* create index for finite element type */
785   if (set_elem_type_index(local_mesh, counter, elem_new2old, elem_old2new)) {
786     return -1;
787   }
788 
789   /* reorder relevant arrays */
790   if (old2new_elem_info(local_mesh, elem_new2old, elem_old2new)) {
791     return -1;
792   }
793 
794   HECMW_free(counter);
795   HECMW_free(elem_new2old);
796   HECMW_free(elem_old2new);
797 
798   return 0;
799 }
800 
801 /*                                                                            */
802 /*  reorder node accordint to DOF                                             */
803 /*                                                                            */
804 
805 /*----------------------------------------------------------------------------*/
806 /*  mask node according to DOF                                                */
807 /*----------------------------------------------------------------------------*/
mask_node_dof_inner(struct hecmwST_local_mesh * local_mesh,char * node_flag,const int is,const int ie,const int n_comp,const int n_dof)808 static int mask_node_dof_inner(struct hecmwST_local_mesh *local_mesh,
809                                char *node_flag, const int is, const int ie,
810                                const int n_comp, const int n_dof) {
811   int nidx, node;
812   int i, j;
813 
814   for (i = is; i < ie; i++) {
815     nidx = local_mesh->elem_node_index[i];
816 
817     for (j = 0; j < n_comp; j++) {
818       node = local_mesh->elem_node_item[nidx + j];
819       MASK_BIT(node_flag[node - 1], n_dof);
820     }
821   }
822 
823   return 0;
824 }
825 
826 /*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
mask_node_dof(struct hecmwST_local_mesh * local_mesh,char * node_flag)827 static int mask_node_dof(struct hecmwST_local_mesh *local_mesh,
828                          char *node_flag) {
829   int n_comp, n_dof;
830   int is, ie;
831   int i;
832   n_dof = -1;
833 
834   for (i = 0; i < local_mesh->n_elem_type; i++) {
835     is     = local_mesh->elem_type_index[i];
836     ie     = local_mesh->elem_type_index[i + 1];
837     n_comp = HECMW_get_max_node(local_mesh->elem_type_item[i]);
838 
839     switch (local_mesh->elem_type_item[i]) {
840       /* line element */
841       case HECMW_ETYPE_ROD1:
842       case HECMW_ETYPE_ROD2:
843         n_dof = BIT_DOF_TWO;
844         break;
845 
846       /* surface element */
847       case HECMW_ETYPE_TRI1:
848       case HECMW_ETYPE_TRI2:
849       case HECMW_ETYPE_TRI22:
850       case HECMW_ETYPE_QUA1:
851       case HECMW_ETYPE_QUA2:
852         n_dof = BIT_DOF_TWO;
853         break;
854 
855       /* solid element */
856       case HECMW_ETYPE_TET1:
857       case HECMW_ETYPE_TET2:
858       case HECMW_ETYPE_TET22:
859       case HECMW_ETYPE_PRI1:
860       case HECMW_ETYPE_PRI2:
861       case HECMW_ETYPE_HEX1:
862       case HECMW_ETYPE_HEX2:
863       case HECMW_ETYPE_PYR1:
864       case HECMW_ETYPE_PYR2:
865       case HECMW_ETYPE_ROD31:
866         n_dof = BIT_DOF_THREE;
867         break;
868       case HECMW_ETYPE_TET1_4:
869       case HECMW_ETYPE_HEX1_4:
870         n_dof = BIT_DOF_FOUR;
871         break;
872 
873       /* master-slave type element */
874       case HECMW_ETYPE_MST1:
875       case HECMW_ETYPE_MST2:
876       case HECMW_ETYPE_MSQ1:
877       case HECMW_ETYPE_MSQ2:
878         n_dof = BIT_DOF_THREE;
879         break;
880 
881       /* interface element */
882       case HECMW_ETYPE_JTB1:
883       case HECMW_ETYPE_JTT1:
884       case HECMW_ETYPE_JTT2:
885       case HECMW_ETYPE_JTQ1:
886       case HECMW_ETYPE_JTQ2:
887         n_dof = BIT_DOF_THREE;
888         break;
889 
890       /* beam element */
891       case HECMW_ETYPE_BEM1:
892       case HECMW_ETYPE_BEM2:
893         n_dof = BIT_DOF_SIX;
894         break;
895 
896       /* surface shell element */
897       case HECMW_ETYPE_SHT1:
898       case HECMW_ETYPE_SHT2:
899       case HECMW_ETYPE_SHQ1:
900       case HECMW_ETYPE_SHQ2:
901       case HECMW_ETYPE_SHQ3:
902         n_dof = BIT_DOF_SIX;
903         break;
904 
905       /* surface shell element (351 361) */
906       case HECMW_ETYPE_SHT6:
907       case HECMW_ETYPE_SHQ8:
908       case HECMW_ETYPE_BEM3:
909         n_dof = BIT_DOF_THREE;
910         break;
911 
912       /* patch element */
913       case HECMW_ETYPE_PTT1:
914       case HECMW_ETYPE_PTT2:
915       case HECMW_ETYPE_PTQ1:
916       case HECMW_ETYPE_PTQ2:
917         n_dof = BIT_DOF_THREE;
918         break;
919 
920       /* link element for MPC */
921       case HECMW_ETYPE_LN11:
922       case HECMW_ETYPE_LN12:
923       case HECMW_ETYPE_LN13:
924       case HECMW_ETYPE_LN14:
925       case HECMW_ETYPE_LN15:
926       case HECMW_ETYPE_LN16:
927       case HECMW_ETYPE_LN21:
928       case HECMW_ETYPE_LN22:
929       case HECMW_ETYPE_LN23:
930       case HECMW_ETYPE_LN24:
931       case HECMW_ETYPE_LN25:
932       case HECMW_ETYPE_LN26:
933       case HECMW_ETYPE_LN31:
934       case HECMW_ETYPE_LN32:
935       case HECMW_ETYPE_LN33:
936       case HECMW_ETYPE_LN34:
937       case HECMW_ETYPE_LN35:
938       case HECMW_ETYPE_LN36:
939       case HECMW_ETYPE_LN41:
940       case HECMW_ETYPE_LN42:
941       case HECMW_ETYPE_LN43:
942       case HECMW_ETYPE_LN44:
943       case HECMW_ETYPE_LN45:
944       case HECMW_ETYPE_LN46:
945       case HECMW_ETYPE_LN51:
946       case HECMW_ETYPE_LN52:
947       case HECMW_ETYPE_LN53:
948       case HECMW_ETYPE_LN54:
949       case HECMW_ETYPE_LN55:
950       case HECMW_ETYPE_LN56:
951       case HECMW_ETYPE_LN61:
952       case HECMW_ETYPE_LN62:
953       case HECMW_ETYPE_LN63:
954       case HECMW_ETYPE_LN64:
955       case HECMW_ETYPE_LN65:
956       case HECMW_ETYPE_LN66:
957         n_dof = BIT_DOF_TWO;
958         break;
959 
960       default:
961         HECMW_assert(0);
962     }
963 
964     if (mask_node_dof_inner(local_mesh, node_flag, is, ie, n_comp, n_dof))
965       return -1;
966   }
967 
968   return 0;
969 }
970 
971 /*----------------------------------------------------------------------------*/
972 /*  reorder node according to DOF & create conversion table                   */
973 /*----------------------------------------------------------------------------*/
reorder_node_dof(struct hecmwST_local_mesh * local_mesh,char * node_flag,int * node_new2old,int * node_old2new,char * dof_flag,int * n_dof_tot)974 static int reorder_node_dof(struct hecmwST_local_mesh *local_mesh,
975                             char *node_flag, int *node_new2old,
976                             int *node_old2new, char *dof_flag, int *n_dof_tot) {
977   int counter = 0;
978   int i;
979 
980   /* six DOF */
981   for (i = 0; i < local_mesh->n_node; i++) {
982     if (EVAL_BIT(node_flag[i], BIT_DOF_SIX)) {
983       node_old2new[i]       = counter + 1;
984       node_new2old[counter] = i + 1;
985       counter++;
986 
987       (n_dof_tot[HECMW_MESH_DOF_SIX])++;
988       MASK_BIT(*dof_flag, BIT_DOF_SIX);
989       CLEAR_BIT(node_flag[i], BIT_DOF_ALL);
990     }
991   }
992 
993   /* four DOF */
994   for (i = 0; i < local_mesh->n_node; i++) {
995     if (EVAL_BIT(node_flag[i], BIT_DOF_FOUR)) {
996       node_old2new[i]       = counter + 1;
997       node_new2old[counter] = i + 1;
998       counter++;
999 
1000       (n_dof_tot[HECMW_MESH_DOF_FOUR])++;
1001       MASK_BIT(*dof_flag, BIT_DOF_FOUR);
1002       CLEAR_BIT(node_flag[i], BIT_DOF_ALL);
1003     }
1004   }
1005 
1006   /* three DOF */
1007   for (i = 0; i < local_mesh->n_node; i++) {
1008     if (EVAL_BIT(node_flag[i], BIT_DOF_THREE)) {
1009       node_old2new[i]       = counter + 1;
1010       node_new2old[counter] = i + 1;
1011       counter++;
1012 
1013       (n_dof_tot[HECMW_MESH_DOF_THREE])++;
1014       MASK_BIT(*dof_flag, BIT_DOF_THREE);
1015       CLEAR_BIT(node_flag[i], BIT_DOF_ALL);
1016     }
1017   }
1018 
1019   /* two DOF */
1020   for (i = 0; i < local_mesh->n_node; i++) {
1021     if (EVAL_BIT(node_flag[i], BIT_DOF_TWO)) {
1022       node_old2new[i]       = counter + 1;
1023       node_new2old[counter] = i + 1;
1024       counter++;
1025 
1026       (n_dof_tot[HECMW_MESH_DOF_TWO])++;
1027       MASK_BIT(*dof_flag, BIT_DOF_TWO);
1028       CLEAR_BIT(node_flag[i], BIT_DOF_ALL);
1029     }
1030   }
1031 
1032   HECMW_assert(counter == local_mesh->n_node);
1033 
1034   return 0;
1035 }
1036 
1037 /*----------------------------------------------------------------------------*/
1038 /*  reorder node according to DOF with EQUATION_BLOCK                         */
1039 /*----------------------------------------------------------------------------*/
1040 #if 0  /* commented out by K.Goto; begin */
1041 static int
1042 mask_eqn_block( struct hecmwST_local_mesh *local_mesh,
1043                 char *node_flag, char *block_flag, const int eqn_block_idx )
1044 {
1045   struct hecmwST_node_grp *grp=local_mesh->node_group;
1046   int dof_max;
1047   int i, j, js;
1048 
1049   for( js=0, i=grp->grp_index[eqn_block_idx]; i<grp->grp_index[eqn_block_idx+1]; i++ ) {
1050     for( dof_max=0, j=js; j<grp->grp_item[i]; j++ ) {
1051       if( EVAL_BIT( node_flag[j], BIT_DOF_TWO ) )   MASK_BIT( block_flag[i-grp->grp_index[eqn_block_idx]], BIT_DOF_TWO );
1052       if( EVAL_BIT( node_flag[j], BIT_DOF_THREE ) ) MASK_BIT( block_flag[i-grp->grp_index[eqn_block_idx]], BIT_DOF_THREE );
1053       if( EVAL_BIT( node_flag[j], BIT_DOF_SIX ) )   MASK_BIT( block_flag[i-grp->grp_index[eqn_block_idx]], BIT_DOF_SIX );
1054       if( EVAL_BIT( node_flag[j], BIT_DOF_FOUR) )   MASK_BIT( block_flag[i-grp->grp_index[eqn_block_idx]], BIT_DOF_FOUR );
1055     }
1056     js = grp->grp_item[i];
1057   }
1058 
1059   return 0;
1060 }
1061 
1062 /*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1063 static int
1064 reorder_node_dof_4mpc_inner( struct hecmwST_local_mesh *local_mesh,
1065                              char *node_flag, char *block_flag,
1066                              int *node_new2old, int *node_old2new, int *block_old2new,
1067                              const int n_eqn_block, const int eqn_block_idx,
1068                              char *dof_flag, int *n_dof_tot )
1069 {
1070   struct hecmwST_node_grp *grp=local_mesh->node_group;
1071   int idx=grp->grp_index[eqn_block_idx];
1072   int counter=0, blocks=0;
1073   int i, j, js;
1074 
1075   /* six DOF */
1076   for( js=0, i=0; i<n_eqn_block; i++ ) {
1077     if( EVAL_BIT( block_flag[i], BIT_DOF_SIX ) ) {
1078       block_old2new[i] = blocks++;
1079 
1080       for( j=js; j<grp->grp_item[idx+i]; j++ ) {
1081         node_old2new[j]       = counter+1;
1082         node_new2old[counter] = j+1;
1083         counter++;
1084 
1085         (n_dof_tot[HECMW_MESH_DOF_SIX])++;
1086       }
1087       MASK_BIT( *dof_flag, BIT_DOF_SIX );
1088       CLEAR_BIT( block_flag[i], BIT_DOF_ALL );
1089     }
1090     js = grp->grp_item[idx+i];
1091   }
1092 
1093   /* four DOF */
1094   for( js=0, i=0; i<n_eqn_block; i++ ) {
1095     if( EVAL_BIT( block_flag[i], BIT_DOF_FOUR ) ) {
1096       block_old2new[i] = blocks++;
1097 
1098       for( j=js; j<grp->grp_item[idx+i]; j++ ) {
1099         node_old2new[j]       = counter+1;
1100         node_new2old[counter] = j+1;
1101         counter++;
1102 
1103         (n_dof_tot[HECMW_MESH_DOF_FOUR])++;
1104       }
1105       MASK_BIT( *dof_flag, BIT_DOF_FOUR );
1106       CLEAR_BIT( block_flag[i], BIT_DOF_ALL );
1107     }
1108     js = grp->grp_item[idx+i];
1109   }
1110 
1111   /* three DOF */
1112   for( js=0, i=0; i<n_eqn_block; i++ ) {
1113     if( EVAL_BIT( block_flag[i], BIT_DOF_THREE ) ) {
1114       block_old2new[i] = blocks++;
1115 
1116       for( j=js; j<grp->grp_item[idx+i]; j++ ) {
1117         node_old2new[j]       = counter+1;
1118         node_new2old[counter] = j+1;
1119         counter++;
1120 
1121         (n_dof_tot[HECMW_MESH_DOF_THREE])++;
1122       }
1123       MASK_BIT( *dof_flag, BIT_DOF_THREE );
1124       CLEAR_BIT( block_flag[i], BIT_DOF_ALL );
1125     }
1126     js = grp->grp_item[idx+i];
1127   }
1128 
1129   /* two DOF */
1130   for( js=0, i=0; i<n_eqn_block; i++ ) {
1131     if( EVAL_BIT( block_flag[i], BIT_DOF_TWO ) ) {
1132       block_old2new[i] = blocks++;
1133 
1134       for( j=js; j<grp->grp_item[idx+i]; j++ ) {
1135         node_old2new[j]       = counter+1;
1136         node_new2old[counter] = j+1;
1137         counter++;
1138 
1139         (n_dof_tot[HECMW_MESH_DOF_TWO])++;
1140       }
1141       MASK_BIT( *dof_flag, BIT_DOF_TWO );
1142       CLEAR_BIT( block_flag[i], BIT_DOF_ALL );
1143     }
1144     js = grp->grp_item[idx+i];
1145   }
1146 
1147   HECMW_assert( counter == local_mesh->n_node );
1148 
1149   return 0;
1150 }
1151 
1152 /*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1153 static int
1154 reorder_node_dof_4mpc( struct hecmwST_local_mesh *local_mesh,
1155                        char *node_flag, int *node_new2old, int *node_old2new,
1156                        char *dof_flag, int *n_dof_tot )
1157 {
1158   struct hecmwST_node_grp *grp=local_mesh->node_group;
1159   char *block_flag;
1160   int *block_old2new;
1161   int n_eqn_block, eqn_block_idx;
1162 
1163   /* group id of EQUATION_BLOCK */
1164   eqn_block_idx = get_eqn_block_idx( local_mesh );
1165   if( eqn_block_idx < 0 ) {
1166     HECMW_print_msg( HECMW_LOG_WARN, HECMW_COMMON_W_NO_EQN_BLOCK, "");
1167     return 1;
1168   }
1169 
1170   /* number of EQUATION_BLOCKs */
1171   n_eqn_block = grp->grp_index[eqn_block_idx+1] - grp->grp_index[eqn_block_idx];
1172 
1173   /* allocation */
1174   block_flag = (char *)HECMW_calloc( n_eqn_block, sizeof(char) );
1175   if( block_flag == NULL ) {
1176     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1177     return -1;
1178   }
1179 
1180   block_old2new = (int *)HECMW_malloc( sizeof(int)*n_eqn_block );
1181   if( block_old2new == NULL ) {
1182     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1183     return -1;
1184   }
1185 
1186   /* mask EQUATION_BLOCK according to DOF */
1187   if(mask_eqn_block( local_mesh, node_flag, block_flag, eqn_block_idx )) return -1;
1188 
1189   /* reorder nodes accordint to DOF with EQUATION_BLOCK */
1190   if(reorder_node_dof_4mpc_inner( local_mesh, node_flag, block_flag, node_new2old, node_old2new, block_old2new, n_eqn_block, eqn_block_idx, dof_flag, n_dof_tot )) return -1;
1191 
1192   /* reorder relevant arrays */
1193   if(old2new_eqn_block_info( local_mesh, block_old2new )) return -1;
1194 
1195   /* free */
1196   HECMW_free( block_flag );
1197   HECMW_free( block_old2new );
1198 
1199   return 0;
1200 }
1201 #endif /* commented out by K.Goto; end */
1202 
1203 /*----------------------------------------------------------------------------*/
1204 /*  count number of DOF groups & set maximal number of DOF                    */
1205 /*----------------------------------------------------------------------------*/
count_n_dof_grp(struct hecmwST_local_mesh * local_mesh,char * dof_flag)1206 static int count_n_dof_grp(struct hecmwST_local_mesh *local_mesh,
1207                            char *dof_flag) {
1208   local_mesh->n_dof_grp = 0;
1209 
1210   /* two DOF */
1211   if (EVAL_BIT(*dof_flag, BIT_DOF_TWO)) {
1212     (local_mesh->n_dof_grp)++;
1213     local_mesh->n_dof = HECMW_MESH_DOF_TWO;
1214   }
1215 
1216   /* three DOF */
1217   if (EVAL_BIT(*dof_flag, BIT_DOF_THREE)) {
1218     (local_mesh->n_dof_grp)++;
1219     local_mesh->n_dof = HECMW_MESH_DOF_THREE;
1220   }
1221 
1222   /* four DOF */
1223   if (EVAL_BIT(*dof_flag, BIT_DOF_FOUR)) {
1224     (local_mesh->n_dof_grp)++;
1225     local_mesh->n_dof = HECMW_MESH_DOF_FOUR;
1226   }
1227 
1228   /* six DOF */
1229   if (EVAL_BIT(*dof_flag, BIT_DOF_SIX)) {
1230     (local_mesh->n_dof_grp)++;
1231     local_mesh->n_dof = HECMW_MESH_DOF_SIX;
1232   }
1233 
1234   HECMW_assert(local_mesh->n_dof_grp > 0 &&
1235                local_mesh->n_dof_grp <= HECMW_MESH_DOF_TOT);
1236 
1237   return 0;
1238 }
1239 
1240 /*----------------------------------------------------------------------------*/
1241 /*  create index for DOF group & its value                                    */
1242 /*----------------------------------------------------------------------------*/
create_node_dof_item(struct hecmwST_local_mesh * local_mesh,char * dof_flag,int * n_dof_tot)1243 static int create_node_dof_item(struct hecmwST_local_mesh *local_mesh,
1244                                 char *dof_flag, int *n_dof_tot) {
1245   int counter = 0;
1246 
1247   /* allocation */
1248   if (local_mesh->node_dof_index) {
1249     HECMW_free(local_mesh->node_dof_index);
1250   }
1251   local_mesh->node_dof_index =
1252       (int *)HECMW_calloc(local_mesh->n_dof_grp + 1, sizeof(int));
1253   if (local_mesh->node_dof_index == NULL) {
1254     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1255     return -1;
1256   }
1257 
1258   if (local_mesh->node_dof_item) {
1259     HECMW_free(local_mesh->node_dof_item);
1260   }
1261   local_mesh->node_dof_item =
1262       (int *)HECMW_malloc(sizeof(int) * local_mesh->n_dof_grp);
1263   if (local_mesh->node_dof_item == NULL) {
1264     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1265     return -1;
1266   }
1267 
1268   /* six DOF */
1269   if (EVAL_BIT(*dof_flag, BIT_DOF_SIX)) {
1270     local_mesh->node_dof_index[counter + 1] =
1271         local_mesh->node_dof_index[counter] + n_dof_tot[HECMW_MESH_DOF_SIX];
1272     local_mesh->node_dof_item[counter] = HECMW_MESH_DOF_SIX;
1273     counter++;
1274   }
1275 
1276   /* four DOF */
1277   if (EVAL_BIT(*dof_flag, BIT_DOF_FOUR)) {
1278     local_mesh->node_dof_index[counter + 1] =
1279         local_mesh->node_dof_index[counter] + n_dof_tot[HECMW_MESH_DOF_FOUR];
1280     local_mesh->node_dof_item[counter] = HECMW_MESH_DOF_FOUR;
1281     counter++;
1282   }
1283 
1284   /* three DOF */
1285   if (EVAL_BIT(*dof_flag, BIT_DOF_THREE)) {
1286     local_mesh->node_dof_index[counter + 1] =
1287         local_mesh->node_dof_index[counter] + n_dof_tot[HECMW_MESH_DOF_THREE];
1288     local_mesh->node_dof_item[counter] = HECMW_MESH_DOF_THREE;
1289     counter++;
1290   }
1291 
1292   /* two DOF */
1293   if (EVAL_BIT(*dof_flag, BIT_DOF_TWO)) {
1294     local_mesh->node_dof_index[counter + 1] =
1295         local_mesh->node_dof_index[counter] + n_dof_tot[HECMW_MESH_DOF_TWO];
1296     local_mesh->node_dof_item[counter] = HECMW_MESH_DOF_TWO;
1297     counter++;
1298   }
1299 
1300   HECMW_assert(counter == local_mesh->n_dof_grp);
1301   HECMW_assert(local_mesh->node_dof_index[local_mesh->n_dof_grp] ==
1302                local_mesh->n_node);
1303 
1304   return 0;
1305 }
1306 
1307 /*============================================================================*/
1308 /*  reorder nodes according to DOF                                            */
1309 /*============================================================================*/
HECMW_reorder_node_dof(struct hecmwST_local_mesh * local_mesh)1310 extern int HECMW_reorder_node_dof(struct hecmwST_local_mesh *local_mesh) {
1311   int *node_new2old,
1312       *node_old2new; /* conversion table between old and new node id */
1313   char *node_flag;   /* for masking */
1314   char dof_flag = '\0';
1315   int n_dof_tot[HECMW_MESH_DOF_MAX + 1];
1316   int i;
1317 
1318   HECMW_assert(local_mesh);
1319 
1320   /* allocation */
1321   node_flag = (char *)HECMW_calloc(local_mesh->n_node, sizeof(char));
1322   if (node_flag == NULL) {
1323     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1324     return -1;
1325   }
1326 
1327   node_new2old = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node);
1328   if (node_new2old == NULL) {
1329     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1330     return -1;
1331   }
1332 
1333   node_old2new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node);
1334   if (node_old2new == NULL) {
1335     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1336     return -1;
1337   }
1338 
1339   for (i = 0; i < HECMW_MESH_DOF_MAX + 1; i++) {
1340     n_dof_tot[i] = 0;
1341   }
1342 
1343   /* mask node */
1344   if (mask_node_dof(local_mesh, node_flag)) {
1345     return -1;
1346   }
1347 
1348   /* reorder node */
1349   /* commented out by K.Goto; begin */
1350   /*
1351   if( local_mesh->mpc->n_mpc ) {
1352     int rtc;
1353     if((rtc = reorder_node_dof_4mpc( local_mesh, node_flag, node_new2old,
1354   node_old2new, &dof_flag, n_dof_tot )) < 0) {
1355       return -1;
1356         }
1357     if( rtc ) {
1358       if(reorder_node_dof( local_mesh, node_flag, node_new2old, node_old2new,
1359   &dof_flag, n_dof_tot )) {
1360         return -1;
1361           }
1362     }
1363   } else {
1364   */
1365   /* commented out by K.Goto; end */
1366   if (reorder_node_dof(local_mesh, node_flag, node_new2old, node_old2new,
1367                        &dof_flag, n_dof_tot)) {
1368     return -1;
1369   }
1370   /* commented out by K.Goto; begin */
1371   /*
1372   }
1373   */
1374   /* commented out by K.Goto; end */
1375 
1376   HECMW_free(node_flag);
1377 
1378   /* create DOF information */
1379   if (count_n_dof_grp(local_mesh, &dof_flag)) {
1380     return -1;
1381   }
1382 
1383   if (create_node_dof_item(local_mesh, &dof_flag, n_dof_tot)) {
1384     return -1;
1385   }
1386 
1387   /* reorder relevant arrays */
1388   if (old2new_node_info(local_mesh, node_new2old, node_old2new)) {
1389     return -1;
1390   }
1391 
1392   /* free */
1393   HECMW_free(node_new2old);
1394   HECMW_free(node_old2new);
1395 
1396   return 0;
1397 }
1398 
1399 /*                                                                            */
1400 /*  reorder node according to MPC group                                       */
1401 /*                                                                            */
1402 
1403 /*----------------------------------------------------------------------------*/
1404 /*  mask MPC group                                                            */
1405 /*----------------------------------------------------------------------------*/
set_mpc_block(struct hecmwST_local_mesh * local_mesh,int * mpc_node_flag,int * mpc_index_flag)1406 static int set_mpc_block(struct hecmwST_local_mesh *local_mesh,
1407                          int *mpc_node_flag, int *mpc_index_flag) {
1408   int node, min_group, group;
1409   int i, j, js, je;
1410 
1411   for (i = 0; i < local_mesh->mpc->n_mpc; i++) {
1412     js = local_mesh->mpc->mpc_index[i];
1413     je = local_mesh->mpc->mpc_index[i + 1];
1414 
1415     /* MPC???롼?פγ???????????륰?롼?פ?mpc_node_flag ?˳??   */
1416     /* ??PC???롼?פν????????å???mpc_index_flag ?˳?? (??? */
1417     min_group = i;
1418     for (j = js; j < je; j++) {
1419       node = local_mesh->mpc->mpc_item[j];
1420 
1421       /* ?????Ƥ?????????????Υ??롼?פ˽?????Ƥ????? ??group >= 0 */
1422       /*                       ¾?Υ??롼?פ˽?????Ƥ??ʤ?????group <  0 */
1423       group                   = mpc_node_flag[node - 1];
1424       mpc_node_flag[node - 1] = (group < 0) ? i : group;
1425 
1426       /* i????PC???롼?פγ????ν?????륰?롼?פκǾ????롼?פ?? */
1427       min_group = (mpc_index_flag[mpc_node_flag[node - 1]] < min_group)
1428                       ? mpc_index_flag[mpc_node_flag[node - 1]]
1429                       : min_group;
1430     }
1431 
1432     /* i????PC???롼?פγ???????????Ƥ??륰?롼?פϺǾ????롼?פ????*/
1433     for (j = js; j < je; j++) {
1434       node  = local_mesh->mpc->mpc_item[j];
1435       group = mpc_node_flag[node - 1];
1436 
1437       mpc_index_flag[group] = min_group;
1438     }
1439     mpc_index_flag[i] = min_group;
1440   }
1441 
1442   /* ???롼?פ?????????*/
1443   for (i = 0; i < local_mesh->mpc->n_mpc; i++) {
1444     group             = mpc_index_flag[i];
1445     mpc_index_flag[i] = mpc_index_flag[group];
1446   }
1447 
1448   return 0;
1449 }
1450 
1451 /*----------------------------------------------------------------------------*/
1452 /*  count MPC blocks                                                          */
1453 /*----------------------------------------------------------------------------*/
count_mpc_block(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block,int * mpc_index_flag,int * mpc_group2block)1454 static int count_mpc_block(struct hecmwST_local_mesh *local_mesh,
1455                            struct equation_block *eqn_block,
1456                            int *mpc_index_flag, int *mpc_group2block) {
1457   int *n_block;
1458   int block, counter;
1459   int i;
1460 
1461   /* allocation */
1462   n_block = (int *)HECMW_calloc(local_mesh->mpc->n_mpc, sizeof(int));
1463   if (n_block == NULL) {
1464     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1465     return -1;
1466   }
1467 
1468   /* count MPC groups in each MPC block */
1469   for (i = 0; i < local_mesh->mpc->n_mpc; i++) {
1470     block = mpc_index_flag[i];
1471     n_block[block]++;
1472   }
1473 
1474   /* count MPC blocks                                              */
1475   /* create conversion table from "MPC group ID" to "MPC block ID" */
1476   for (counter = 0, i = 0; i < local_mesh->mpc->n_mpc; i++) {
1477     if (n_block[i]) mpc_group2block[i] = counter++;
1478   }
1479 
1480   /* number of MPC blocks */
1481   eqn_block->n_mpc_block = counter;
1482   eqn_block->n_eqn_block = counter;
1483 
1484   /* free */
1485   HECMW_free(n_block);
1486 
1487   return 0;
1488 }
1489 
1490 /*----------------------------------------------------------------------------*/
1491 /*  count EQUATION_BLOCKs                                                     */
1492 /*----------------------------------------------------------------------------*/
count_eqn_block(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block,int * mpc_node_flag)1493 static int count_eqn_block(struct hecmwST_local_mesh *local_mesh,
1494                            struct equation_block *eqn_block,
1495                            int *mpc_node_flag) {
1496   int i;
1497 
1498   for (i = 0; i < local_mesh->n_node; i++) {
1499     if (mpc_node_flag[i] < 0) (eqn_block->n_eqn_block)++;
1500   }
1501 
1502   return 0;
1503 }
1504 
1505 /*----------------------------------------------------------------------------*/
1506 /*  set EQUATION_BLOCK id to node                                             */
1507 /*----------------------------------------------------------------------------*/
set_eqn_block_of_node(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block,int * mpc_node_flag,int * mpc_index_flag,int * mpc_group2block)1508 static int set_eqn_block_of_node(struct hecmwST_local_mesh *local_mesh,
1509                                  struct equation_block *eqn_block,
1510                                  int *mpc_node_flag, int *mpc_index_flag,
1511                                  int *mpc_group2block) {
1512   int counter;
1513   int i;
1514 
1515   for (counter = eqn_block->n_mpc_block, i = 0; i < local_mesh->n_node; i++) {
1516     if (mpc_node_flag[i] >= 0) {
1517       mpc_node_flag[i] = mpc_group2block[mpc_index_flag[mpc_node_flag[i]]];
1518     } else {
1519       mpc_node_flag[i] = counter++;
1520     }
1521   }
1522 
1523   HECMW_assert(counter == eqn_block->n_eqn_block);
1524 
1525   return 0;
1526 }
1527 
1528 /*----------------------------------------------------------------------------*/
1529 /*   create EQUATION_BLOCK                                                    */
1530 /*----------------------------------------------------------------------------*/
create_eqn_block_index(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block,int * mpc_node_flag,int * mpc_group2block)1531 static int create_eqn_block_index(struct hecmwST_local_mesh *local_mesh,
1532                                   struct equation_block *eqn_block,
1533                                   int *mpc_node_flag, int *mpc_group2block) {
1534   int *n_block;
1535   int i;
1536 
1537   /* allocation */
1538   n_block = (int *)HECMW_calloc(eqn_block->n_eqn_block, sizeof(int));
1539   if (n_block == NULL) {
1540     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1541     return -1;
1542   }
1543 
1544   eqn_block->eqn_block_index =
1545       (int *)HECMW_calloc(eqn_block->n_eqn_block + 1, sizeof(int));
1546   if (eqn_block->eqn_block_index == NULL) {
1547     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1548     return -1;
1549   }
1550 
1551   /* count nodes in each EQUATION_BLOCK */
1552   for (i = 0; i < local_mesh->n_node; i++) {
1553     (n_block[mpc_node_flag[i]])++;
1554   }
1555 
1556   /* create index for EQUATION_BLOCK */
1557   for (i = 0; i < eqn_block->n_eqn_block; i++) {
1558     eqn_block->eqn_block_index[i + 1] =
1559         eqn_block->eqn_block_index[i] + n_block[i];
1560   }
1561 
1562   /* free */
1563   HECMW_free(n_block);
1564 
1565   HECMW_assert(eqn_block->eqn_block_index[eqn_block->n_eqn_block] ==
1566                local_mesh->n_node);
1567 
1568   return 0;
1569 }
1570 
1571 /*----------------------------------------------------------------------------*/
1572 /*  create conversion table of node                                           */
1573 /*----------------------------------------------------------------------------*/
create_old2new_node(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block,int * mpc_node_flag,int * node_old2new,int * node_new2old)1574 static int create_old2new_node(struct hecmwST_local_mesh *local_mesh,
1575                                struct equation_block *eqn_block,
1576                                int *mpc_node_flag, int *node_old2new,
1577                                int *node_new2old) {
1578   int *n_block;
1579   int new_id;
1580   int i;
1581 
1582   n_block = (int *)HECMW_calloc(eqn_block->n_eqn_block, sizeof(int));
1583   if (n_block == NULL) {
1584     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1585     return -1;
1586   }
1587 
1588   for (i = 0; i < local_mesh->n_node; i++) {
1589     new_id = eqn_block->eqn_block_index[mpc_node_flag[i]] +
1590              n_block[mpc_node_flag[i]];
1591     node_old2new[i]      = new_id + 1;
1592     node_new2old[new_id] = i + 1;
1593     (n_block[mpc_node_flag[i]])++;
1594   }
1595 
1596   HECMW_free(n_block);
1597 
1598   return 0;
1599 }
1600 
1601 /*----------------------------------------------------------------------------*/
1602 /*  reconstruct node group informaion                                         */
1603 /*----------------------------------------------------------------------------*/
reconstruct_node_grp(struct hecmwST_local_mesh * local_mesh,struct equation_block * eqn_block)1604 static int reconstruct_node_grp(struct hecmwST_local_mesh *local_mesh,
1605                                 struct equation_block *eqn_block) {
1606   struct hecmwST_node_grp *grp = local_mesh->node_group;
1607   int i;
1608 
1609   /* number of node groups */
1610   (grp->n_grp)++;
1611 
1612   /* index for node group */
1613   grp->grp_index =
1614       (int *)HECMW_realloc(grp->grp_index, sizeof(int) * (grp->n_grp + 1));
1615   if (grp->grp_index == NULL) {
1616     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1617     return -1;
1618   }
1619 
1620   grp->grp_index[grp->n_grp] =
1621       grp->grp_index[grp->n_grp - 1] + eqn_block->n_eqn_block;
1622 
1623   /* name of node group */
1624   grp->grp_name =
1625       (char **)HECMW_realloc(grp->grp_name, sizeof(char *) * grp->n_grp);
1626   if (grp->grp_name == NULL) {
1627     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1628     return -1;
1629   }
1630   grp->grp_name[grp->n_grp - 1] =
1631       (char *)HECMW_malloc(sizeof(char *) * (HECMW_NAME_LEN + 1));
1632   if (grp->grp_name[grp->n_grp - 1] == NULL) {
1633     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1634     return -1;
1635   }
1636 
1637   strcpy(grp->grp_name[grp->n_grp - 1], HECMW_COMMON_EQUATION_BLOCK_NAME);
1638 
1639   /* member of node group */
1640   grp->grp_item = (int *)HECMW_realloc(
1641       grp->grp_item, sizeof(int) * grp->grp_index[grp->n_grp]);
1642   if (grp->grp_item == NULL) {
1643     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1644     return -1;
1645   }
1646 
1647   for (i = 0; i < eqn_block->n_eqn_block; i++) {
1648     grp->grp_item[grp->grp_index[grp->n_grp - 1] + i] =
1649         eqn_block->eqn_block_index[i + 1];
1650   }
1651 
1652   return 0;
1653 }
1654 
1655 /*============================================================================*/
1656 /*  reorder node according to MPC group                                       */
1657 /*============================================================================*/
HECMW_reorder_node_mpc(struct hecmwST_local_mesh * local_mesh)1658 extern int HECMW_reorder_node_mpc(struct hecmwST_local_mesh *local_mesh) {
1659   int *mpc_node_flag, *mpc_index_flag, *mpc_group2block;
1660   int *node_old2new, *node_new2old;
1661   int i;
1662   struct equation_block *eqn_block;
1663 
1664   if (local_mesh->mpc->n_mpc == 0) return 0;
1665 
1666   /* allocation */
1667   mpc_node_flag = (int *)HECMW_calloc(local_mesh->n_node, sizeof(int));
1668   if (mpc_node_flag == NULL) {
1669     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1670     return -1;
1671   }
1672   for (i = 0; i < local_mesh->n_node; i++) {
1673     mpc_node_flag[i] = -1;
1674   }
1675 
1676   mpc_index_flag = (int *)HECMW_malloc(sizeof(int) * local_mesh->mpc->n_mpc);
1677   if (mpc_index_flag == NULL) {
1678     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1679     return -1;
1680   }
1681   for (i = 0; i < local_mesh->mpc->n_mpc; i++) {
1682     mpc_index_flag[i] = i;
1683   }
1684 
1685   eqn_block =
1686       (struct equation_block *)HECMW_malloc(sizeof(struct equation_block));
1687   if (eqn_block == NULL) {
1688     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1689     return -1;
1690   }
1691 
1692   /* group together */
1693   if (set_mpc_block(local_mesh, mpc_node_flag, mpc_index_flag)) {
1694     return -1;
1695   }
1696 
1697   /* count MPC blocks */
1698   mpc_group2block = (int *)HECMW_malloc(sizeof(int) * local_mesh->mpc->n_mpc);
1699   if (mpc_group2block == NULL) {
1700     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1701     return -1;
1702   }
1703   for (i = 0; i < local_mesh->mpc->n_mpc; i++) mpc_group2block[i] = -1;
1704 
1705   if (count_mpc_block(local_mesh, eqn_block, mpc_index_flag, mpc_group2block)) {
1706     return -1;
1707   }
1708 
1709   /* count EQUATION_BLOCKs */
1710   if (count_eqn_block(local_mesh, eqn_block, mpc_node_flag)) {
1711     return -1;
1712   }
1713 
1714   /* set EQUATION_BLOCK to node */
1715   if (set_eqn_block_of_node(local_mesh, eqn_block, mpc_node_flag,
1716                             mpc_index_flag, mpc_group2block)) {
1717     return -1;
1718   }
1719 
1720   /* create EQUATION_BLOCK */
1721   if (create_eqn_block_index(local_mesh, eqn_block, mpc_node_flag,
1722                              mpc_group2block)) {
1723     return -1;
1724   }
1725 
1726   HECMW_free(mpc_index_flag);
1727   HECMW_free(mpc_group2block);
1728 
1729   /* create conversion table between old and new node id */
1730   node_old2new = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node);
1731   if (node_old2new == NULL) {
1732     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1733     return -1;
1734   }
1735 
1736   node_new2old = (int *)HECMW_malloc(sizeof(int) * local_mesh->n_node);
1737   if (node_new2old == NULL) {
1738     HECMW_set_error(HECMW_COMMON_E_ALLOCATION, "");
1739     return -1;
1740   }
1741 
1742   if (create_old2new_node(local_mesh, eqn_block, mpc_node_flag, node_old2new,
1743                           node_new2old)) {
1744     return -1;
1745   }
1746 
1747   if (old2new_node_info(local_mesh, node_new2old, node_old2new)) {
1748     return -1;
1749   }
1750 
1751   HECMW_free(mpc_node_flag);
1752   HECMW_free(node_old2new);
1753   HECMW_free(node_new2old);
1754 
1755   /* reconstruct node group */
1756   if (reconstruct_node_grp(local_mesh, eqn_block)) {
1757     return -1;
1758   }
1759 
1760   HECMW_free(eqn_block->eqn_block_index);
1761   HECMW_free(eqn_block);
1762 
1763   return 0;
1764 }
1765 
1766 /*                                                                            */
1767 /*  reorder node & element                                                    */
1768 /*                                                                            */
1769 
HECMW_reorder(struct hecmwST_local_mesh * local_mesh)1770 extern int HECMW_reorder(struct hecmwST_local_mesh *local_mesh) {
1771   HECMW_assert(local_mesh);
1772 
1773   /* reorder element accordint to finite element type */
1774   if (HECMW_reorder_elem_type(local_mesh)) return -1;
1775 
1776   /* reorder node according to MPC group */
1777   /* commented out by K.Goto; begin */
1778   /* if(HECMW_reorder_node_mpc( local_mesh )) return -1; */
1779   /* commented out by K.Goto; end */
1780 
1781   /* reorder node accordint to node dof */
1782   if (HECMW_reorder_node_dof(local_mesh)) return -1;
1783 
1784   return 0;
1785 }
1786