All Classes Files Functions Variables Pages
partCubitMesh.c
Go to the documentation of this file.
1 
121 
127 // Standard C include
128 #include <stdio.h>
129 #include <stdlib.h>
130 #include <string.h>
131 #include <math.h>
132 #include <assert.h>
133 #include <ctype.h>
134 
135 // my include
136 #include "partCubitMesh.h"
137 #include "mpi.h"
138 
139 #include "exodusII.h"
140 
141 
142 /********************************************************** FUNCTIONS DECLARATION **********************************************************/
143 int count_elt(FILE *unit);
144 
145 void setHexaWeight(mesh_t *mesh);
146 
147 void readCubitMeshAbaqus(char *fileinp, mesh_t *mesh, int *number_of_elemnt_block);
148 
149 void readCubitMeshExodus(char *fileinp, mesh_t *mesh, int *number_of_elemnt_block);
150 
151 void read_part_cubit_mesh_(char *fileinp
152  , int *ig_ncpu
153  , int *ig_myrank
154  , int *size_real_t
155  , int *ig_mesh_nnode
156  , int *ig_ncpu_neighbor
157  , int *ig_nhexa
158  , int *ig_nhexa_outer
159  , int *ig_nquad_parax
160  , int *ig_nquad_fsurf
161  , int *ig_hexa_nnode
162  , int *ig_quad_nnode
163  , float *rg_mesh_xmax
164  , float *rg_mesh_ymax
165  , float *rg_mesh_zmax
166  , float *rg_mesh_xmin
167  , float *rg_mesh_ymin
168  , float *rg_mesh_zmin);
169 
170 int fill_mesh_arrays_ (int *ig_ncpu
171  , int *ig_myrank
172  , int *ig_ngll_total
173  , int *ig_nneighbor_all_kind
174  , int *cpu_neighbor
175  , int *ig_cpu_neighbor_info
176  , int *ig_hexa_gnode_glonum
177  , int *ig_quadp_gnode_glonum
178  , int *ig_quadf_gnode_glonum
179  , int *ig_hexa_gll_glonum
180  , int *ig_quadp_gll_glonum
181  , int *ig_quadf_gll_glonum
182  , int *ig_quadp_neighbor_hexa
183  , int *ig_quadp_neighbor_hexaface
184  , int *ig_hexa_material_number);
185 
186 void fill_efispec_arrays( info_t *info
187  , mesh_t *mesh
188  , int rank
189  , int *ig_ngll_total
190  , int *ig_nneighbor_all_kind
191  , int *cpu_neighbor
192  , int *ig_cpu_neighbor_info
193  , int *ig_hexa_gnode_glonum
194  , int *ig_quadp_gnode_glonum
195  , int *ig_quadf_gnode_glonum
196  , int *ig_hexa_gll_glonum
197  , int *ig_quadp_gll_glonum
198  , int *ig_quadf_gll_glonum
199  , int *ig_quadp_neighbor_hexa
200  , int *ig_quadp_neighbor_hexaface
201  , int *ig_hexa_material_number);
202 
203 info_t *getConnInfo(mesh_t *mesh, int rank, int *ig_mesh_nnode);
204 
205 int getProcConn(int ielt, int iproc, proc_info_t *proc_tab, mesh_t *mesh);
206 
207 void setEltConn(int ielt, int iproc, elt_info_t *elt_tab, mesh_t *mesh, proc_info_t *proc_tab);
208 
209 void setHexaConnDetail(int elt_num, int neigh_num, int nb_conn, int *connex_nodes, elt_info_t *elt_tab, mesh_t *mesh);
210 
211 void setQuadConnDetail(int elt_num, int neigh_num, int nb_conn, int *connex_nodes, elt_info_t *elt_tab, mesh_t *mesh);
212 
213 int findFaceOrientation(int *connex_nodes, int num_conn_source, int num_conn_target);
214 
215 int findEdgeOrientation(int *connex_nodes, int num_conn_source, int num_conn_target);
216 
217 char *ltrim(char*s);
218 char *rtrim(char*s);
219 char *trim (char*s);
220 
221 void printMemUsage(mesh_t *mesh);
222 
223 char* rotateVect(char *ortn, char *ret_or, int cw_quarter_rot);
224 
225 void add_conn(topo_t typecon, elt_info_t *elt_tab, int hexa_src, int numcon_src, idx_t num_neigh, int numcon_neigh, int orientation, elt_t type);
226 
227 int getnbneigh(topo_t typecon, elt_info_t *elt_tab, int numcon_src);
228 
229 void free_all_1();
230 
231 void free_all_2();
232 
233 void get_domain_min_max(mesh_t *mesh, float *rg_mesh_xmax, float *rg_mesh_ymax, float *rg_mesh_zmax, float *rg_mesh_xmin, float *rg_mesh_ymin, float *rg_mesh_zmin);
234 
235 // Fortran routines
236 void mod_init_mesh_mp_init_gll_number_(int*,int*);
237 
238 void mod_init_mesh_mp_propagate_gll_nodes_face_(int*, int*, int*, int*, int*);
239 
240 void mod_init_mesh_mp_propagate_gll_nodes_edge_(int*, int*, int*, int*, int*);
241 
242 void mod_init_mesh_mp_propagate_gll_nodes_corner_(int*, int*, int*, int*);
243 
244 void mod_init_mesh_mp_propagate_gll_nodes_quad_(int*, int*, int*, int*, int*);
245 
246 void pass_address_back_(float *, float *, float *);
247 
248 void allocate_node_coord_arrays_();
249 
250 void flush()
251 {
252  fflush(stdout);
253  fflush(stderr);
254 }
255 
256 /********************************************************** GLOBAL VARIABLES **********************************************************/
257 static mesh_t *mesh_ptr;
258 static info_t *info_ptr;
259 
260 static int ncpu_neighbor, nhexa_local;
261 static int nquad_parax;
262 static int nquad_fsurf;
263 static int num_nodes_hexa;
264 
265 float *xcoord;
266 float *ycoord;
267 float *zcoord;
268 
269 /******************************************************** FUNCTIONS DEFINITION ********************************************************/
270 void pass_address_back_(float *rg_gnode_x, float *rg_gnode_y, float *rg_gnode_z)
271 {
272  xcoord = rg_gnode_x;
273  ycoord = rg_gnode_y;
274  zcoord = rg_gnode_z;
275 }
276 
277 void free_all_1()
278 {
279  free(mesh_ptr->xadj);
280  free(mesh_ptr->adjncy);
281  free(mesh_ptr->xadj_hex);
282  free(mesh_ptr->adjncy_hex);
283  free(mesh_ptr->ncoords_x);
284  free(mesh_ptr->ncoords_y);
285  free(mesh_ptr->ncoords_z);
286 }
287 
288 void free_all_2() // à modifier (cf 896 et avant)
289 {
290  int iproc;
291  for (iproc=0; iproc<mesh_ptr->npart; iproc++) {
292  free(info_ptr->proc[iproc].local_elts);
293  free(info_ptr->proc[iproc].connex);
294  }
295  free(info_ptr->proc);
296 
297  int nbelt = mesh_ptr->nh;
298  for(int i=0; i< nbelt; i++) {
299  elt_info_t elt = info_ptr->elt[i];
300 
301  for (int j=0; j<NEDGE; j++) {
302  // conn_info* ptr = &elt->edges[j];
303  conn_info* ptr = &elt.edges[j];
304  if (ptr->defined) {
305  conn_info* next;
306  int first = 1;
307  while(ptr)
308  { next = ptr->next;
309  if (!first) free(ptr);
310  ptr = next;
311  first = 0;
312  }
313  }
314  }
315  for (int j=0; j<NCORNER; j++) {
316  // conn_info* ptr = &elt->corners[j];
317  conn_info* ptr = &elt.corners[j];
318  if (ptr->defined) {
319  conn_info* next;
320  int first = 1;
321  while(ptr)
322  { next = ptr->next;
323  if (!first) free(ptr);
324  ptr = next;
325  first = 0;
326  }
327  }
328  }
329  }
330 
331  free(info_ptr->elt);
332  free(info_ptr);
333 
334  free(mesh_ptr->eptr);
335  free(mesh_ptr->eind);
336  if (mesh_ptr->hex27) {
337  free(mesh_ptr->eptr_27);
338  free(mesh_ptr->eind_27);
339  }
340  free(mesh_ptr->part);
341  free(mesh_ptr->layer);
342  free(mesh_ptr->vwgt);
343  free(mesh_ptr->types);
344 
345  free(mesh_ptr);
346 }
347 
348 /******************************************************** FUNCTIONS *******************************************************************/
349 void read_part_cubit_mesh_( char * prefix
350  , int * ig_ncpu
351  , int * ig_myrank
352  , int * size_real_t
353  , int * ig_mesh_nnode
354  , int * ig_ncpu_neighbor
355  , int * ig_nhexa
356  , int * ig_nhexa_outer
357  , int * ig_nquad_parax
358  , int * ig_nquad_fsurf
359  , int * ig_hexa_nnode
360  , int * ig_quad_nnode
361  , float* rg_mesh_xmax
362  , float* rg_mesh_ymax
363  , float* rg_mesh_zmax
364  , float* rg_mesh_xmin
365  , float* rg_mesh_ymin
366  , float* rg_mesh_zmin)
367 {
368  int rank, iproc;
369  mesh_t *mesh; // -> the mesh structure
370  idx_t ncommon;
371  idx_t pnumflag;
372  idx_t edgecuts;
373  float r_buf[6];
374  int i_buf[8];
375  int number_of_elemnt_block[N_BLOCK_MAX];
376  MPI_Status status;
377  char fileinp[100];
378 
379 
380  mesh = MALLOC(mesh_t, 1);
381  // input
382  mesh->npart = *ig_ncpu;
383  rank = *ig_myrank;
384 
385  if (rank == 0) {
386 
387  MSG("")
388  // 93 because cg_prefix is a 92 car string + \0
389  strncpy(fileinp, prefix, 93);
390  strcat(fileinp, ".ex2");
391  FILE * check = fopen(fileinp, "r");
392  if (check)
393  {
394  fclose(check);
395  MSG("Found exodusII binary mesh file")
396  readCubitMeshExodus(fileinp, mesh, number_of_elemnt_block);
397  }
398  else
399  {
400  strncpy(fileinp, prefix, 93);
401  strcat(fileinp, ".inp");
402  check = fopen(fileinp, "r");
403  if (check) {
404  fclose(check);
405  MSG("Found abaqus ascii mesh file")
406  readCubitMeshAbaqus(fileinp, mesh, number_of_elemnt_block);
407  } else {
408  printf("\nMesh file not found (%s.ex2 or %s.inp)\naborting ...\n\n",prefix,prefix);
409  MPI_Abort(MPI_COMM_WORLD, 1);
410  }
411  }
412 
413 #ifdef VERBOSE
414  MSG("")
415  MSG("Mesh information :")
416  printf("\t%d nodes\n",mesh->nn);
417  printf("\t%d nodes per hexahedron\n",mesh->num_nodes_hexa);
418  printf("\t%d hexahedron elements\n",mesh->nh);
419  printf("\t%d quadrangle elements as boundary absorption\n",mesh->nq_parax);
420  printf("\t%d quadrangle elements as free surface snapshot\n",mesh->nq_surf);
421  MSG("")
422  MSG("Computation started: see *.lst file for more information")
423 #endif
424  get_domain_min_max( mesh, rg_mesh_xmax, rg_mesh_ymax, rg_mesh_zmax, rg_mesh_xmin, rg_mesh_ymin, rg_mesh_zmin);
425 
426  r_buf[0] = *rg_mesh_xmax;
427  r_buf[1] = *rg_mesh_ymax;
428  r_buf[2] = *rg_mesh_zmax;
429  r_buf[3] = *rg_mesh_xmin;
430  r_buf[4] = *rg_mesh_ymin;
431  r_buf[5] = *rg_mesh_zmin;
432  MPI_Bcast( r_buf, 6, MPI_FLOAT, 0, MPI_COMM_WORLD);
433  } else {
434  MPI_Bcast( r_buf, 6, MPI_FLOAT, 0, MPI_COMM_WORLD);
435  *rg_mesh_xmax = r_buf[0];
436  *rg_mesh_ymax = r_buf[1];
437  *rg_mesh_zmax = r_buf[2];
438  *rg_mesh_xmin = r_buf[3];
439  *rg_mesh_ymin = r_buf[4];
440  *rg_mesh_zmin = r_buf[5];
441  }
442 
443  if (rank == 0) {
444 
445  //MSG("Build complete adjacency graph with 4 nodes connections")
446  ncommon = 4;
447  pnumflag = 0;
448  METIS_MeshToDual(&mesh->ne, &mesh->nn, mesh->eptr, mesh->eind, &ncommon, &pnumflag, &mesh->xadj, &mesh->adjncy);
449 
450  //MSG("Build hexahedron adjacency graph with 1 node connection")
451  ncommon = 1;
452  METIS_MeshToDual(&mesh->nh, &mesh->nn, mesh->eptr, mesh->eind, &ncommon, &pnumflag, &mesh->xadj_hex, &mesh->adjncy_hex);
453 
454  //MSG("Set weight to hexa")
455  setHexaWeight(mesh);
456 
457  //MSG("Part mesh using Metis")
458 
459  if (mesh->npart > 1) {
460 
461  METIS_PartGraphRecursive(&mesh->nh, &mesh->ncon, mesh->xadj_hex, mesh->adjncy_hex, mesh->vwgt, NULL, NULL, &mesh->npart, NULL, NULL, NULL, &edgecuts, mesh->part);
462  int numelem = 0;
463 
464  for (int hh=0; hh<mesh->nh; hh++) {
465  if (mesh->part[hh] == rank) numelem++;
466  }
467 
468  } else {
469 
470  memset(mesh->part, 0, mesh->nh*sizeof(idx_t));
471 
472  }
473 
474  //MSG("Get complete connectivity information")
475 
476  // rg_[xyz]_coord_geom_nodes allocated and filled in getConnInfo()
477  // allocation coord_geom_nodes -> peut-etre redecouper ici !!!
478  info_t* info = getConnInfo(mesh, rank, ig_mesh_nnode);
479 
480  for (iproc = 0; iproc < mesh->npart; iproc++) {
481  if (iproc == 0) {
482  //sizes 4 efispec
483  *size_real_t = sizeof(real_t);
484  *ig_nhexa = info->proc[iproc].nb_elt;
485  *ig_nhexa_outer = info->proc[iproc].nb_ext;
486  *ig_nquad_parax = info->proc[iproc].nb_quad_p;
487  *ig_nquad_fsurf = info->proc[iproc].nb_quad_f;
488  *ig_hexa_nnode = mesh->num_nodes_hexa;
489  *ig_ncpu_neighbor = info->proc[iproc].nb_conn;
490  *ig_quad_nnode = mesh->num_node_per_dim*mesh->num_node_per_dim;
491  } else {
492  i_buf[0] = sizeof(real_t);
493  i_buf[1] = info->proc[iproc].nb_elt;
494  i_buf[2] = info->proc[iproc].nb_ext;
495  i_buf[3] = info->proc[iproc].nb_quad_p;
496  i_buf[4] = info->proc[iproc].nb_quad_f;
497  i_buf[5] = mesh->num_nodes_hexa;
498  i_buf[6] = info->proc[iproc].nb_conn;
499  i_buf[7] = mesh->num_node_per_dim*mesh->num_node_per_dim;
500  MPI_Send(i_buf, 8, MPI_INT, iproc, 5, MPI_COMM_WORLD);
501  }
502  }
503  // DAVID : statics ptr for memory freeing
504  mesh_ptr = mesh;
505  info_ptr = info;
506  } else {
507  getConnInfo(mesh, rank, ig_mesh_nnode);
508  MPI_Recv(i_buf, 8, MPI_INT, 0, 5, MPI_COMM_WORLD, &status);
509  *size_real_t = i_buf[0];
510  *ig_nhexa = i_buf[1];
511  *ig_nhexa_outer = i_buf[2];
512  *ig_nquad_parax = i_buf[3];
513  *ig_nquad_fsurf = i_buf[4];
514  *ig_hexa_nnode = i_buf[5];
515  *ig_ncpu_neighbor = i_buf[6];
516  *ig_quad_nnode = i_buf[7];
517  mesh_ptr = mesh;
518  }
519 
520  // enregistrement pour plus tard
521  ncpu_neighbor = *ig_ncpu_neighbor;
522  nhexa_local = *ig_nhexa;
523  nquad_parax = *ig_nquad_parax;
524  nquad_fsurf = *ig_nquad_fsurf;
525  num_nodes_hexa = *ig_hexa_nnode;
526 
527 
528  // free what's not needed anymore to avoid a memory consumption pick as much as possible
529  if (rank == 0) free_all_1();
530 
531  return;
532 }
533 
534 void get_domain_min_max(mesh_t* mesh, float* rg_mesh_xmax, float* rg_mesh_ymax, float* rg_mesh_zmax, float* rg_mesh_xmin, float* rg_mesh_ymin, float* rg_mesh_zmin)
535 {
536  float xmin, ymin, zmin;
537  float xmax, ymax, zmax;
538  xmin = ymin = zmin = +FLT_MAX;
539  xmax = ymax = zmax = -FLT_MAX;
540 
541  for (int node=0; node<mesh->nn; node++) {
542  if (mesh->ncoords_x[node] < xmin) xmin = mesh->ncoords_x[node];
543  if (mesh->ncoords_y[node] < ymin) ymin = mesh->ncoords_y[node];
544  if (mesh->ncoords_z[node] < zmin) zmin = mesh->ncoords_z[node];
545 
546  if (mesh->ncoords_x[node] > xmax) xmax = mesh->ncoords_x[node];
547  if (mesh->ncoords_y[node] > ymax) ymax = mesh->ncoords_y[node];
548  if (mesh->ncoords_z[node] > zmax) zmax = mesh->ncoords_z[node];
549  }
550 
551  *rg_mesh_xmax = xmax;
552  *rg_mesh_ymax = ymax;
553  *rg_mesh_zmax = zmax;
554 
555  *rg_mesh_xmin = xmin;
556  *rg_mesh_ymin = ymin;
557  *rg_mesh_zmin = zmin;
558 }
559 
560 int fill_mesh_arrays_(int* nb_part_,
561  int* rank_,
562  int* ig_ngll_total,
563  int* ig_nneighbor_all_kind,
564  int* cpu_neighbor,
565  int* ig_cpu_neighbor_info,
566  int* ig_hexa_gnode_glonum,
567  int* ig_quadp_gnode_glonum,
568  int* ig_quadf_gnode_glonum,
569  int* ig_hexa_gll_glonum,
570  int* ig_quadp_gll_glonum,
571  int* ig_quadf_gll_glonum,
572  int* ig_quadp_neighbor_hexa,
573  int* ig_quadp_neighbor_hexaface,
574  int* ig_hexa_material_number)
575 {
576 
577  mesh_t* mesh = mesh_ptr;
578  info_t* info = info_ptr;
579 
580  fill_efispec_arrays( info
581  , mesh
582  ,*rank_
583  , ig_ngll_total
584  , ig_nneighbor_all_kind
585  , cpu_neighbor
586  , ig_cpu_neighbor_info
587  , ig_hexa_gnode_glonum
588  , ig_quadp_gnode_glonum
589  , ig_quadf_gnode_glonum
590  , ig_hexa_gll_glonum
591  , ig_quadp_gll_glonum
592  , ig_quadf_gll_glonum
593  , ig_quadp_neighbor_hexa
594  , ig_quadp_neighbor_hexaface
595  , ig_hexa_material_number);
596 
597  if (*rank_ == 0) free_all_2();
598  else free(mesh_ptr);
599 
600  //if (*rank_ == 0) MSG("Mesh partition done.")
601 
602  return 0;
603 }
604 
605 // write domain files for each proc of efispec3D
606 void fill_efispec_arrays( info_t* info
607  , mesh_t* mesh
608  , int rank
609  , int * ig_ngll_total
610  , int * ig_nneighbor_all_kind
611  , int * cpu_neighbor
612  , int * ig_cpu_neighbor_info //*
613  , int * ig_hexa_gnode_glonum //*
614  , int * ig_quadp_gnode_glonum
615  , int * ig_quadf_gnode_glonum
616  , int * ig_hexa_gll_glonum
617  , int * ig_quadp_gll_glonum
618  , int * ig_quadf_gll_glonum
619  , int * ig_quadp_neighbor_hexa //*
620  , int * ig_quadp_neighbor_hexaface //*
621  , int * ig_hexa_material_number) //*
622 {
623  int iproc, ineigh, ielt, i, j, k;
624  int* i_buf, *ibuf_hexa_material_number, *ibuf_hexa_gnode_glonum;
625  proc_info_t* proc_tab;
626  elt_info_t* elt_tab;
627  int* nneighbor_all_kind;
628  int* cpu_neighbor_info;
629  int* ibuf_quadp_neighbor_hexa;
630  int* ibuf_quadp_neighbor_hexaface;
631  int ibuf2[2];
632  MPI_Status status;
633  int icpu, iparax, ifsurf;
634 
635  idx_t *p_eptr, *p_eind;
636  if (mesh->hex27) {
637  p_eptr = mesh->eptr_27;
638  p_eind = mesh->eind_27;
639  } else {
640  p_eptr = mesh->eptr;
641  p_eind = mesh->eind;
642  }
643 
644  if (rank == 0) {
645  proc_tab = info->proc;
646  elt_tab = info->elt;
647 
648  for(iproc = 0; iproc<mesh->npart; iproc++) {
649  if (iproc == 0) {
650  icpu=0;//nb connected proc
651  for(ineigh = 0; ineigh<mesh->npart; ineigh++) {
652  if (proc_tab[iproc].connex[ineigh]) cpu_neighbor[icpu++] = ineigh; // list of connected proc
653  }
654  } else {
655  // envoi
656  i_buf = MALLOC(int, proc_tab[iproc].nb_conn);
657  icpu=0;//nb connected proc
658  for(ineigh = 0; ineigh<mesh->npart; ineigh++) {
659  if (proc_tab[iproc].connex[ineigh]) i_buf[icpu++] = ineigh; // list of connected proc
660  }
661  MPI_Send(i_buf, proc_tab[iproc].nb_conn, MPI_INT, iproc, 10, MPI_COMM_WORLD);
662  free(i_buf);
663  }
664  }
665  } else {
666  MPI_Recv(cpu_neighbor, ncpu_neighbor, MPI_INT, 0, 10, MPI_COMM_WORLD, &status);
667  }
668 
669  for(iproc = 0; iproc<mesh->npart; iproc++) {
670  int nneighbor_all_kind = 0;
671  int nb_items = 0;
672  if (rank == 0) {
673  if (iproc != 0) {
674  ibuf_hexa_material_number = MALLOC(int, proc_tab[iproc].nb_elt);
675  ibuf_hexa_gnode_glonum = MALLOC(int, proc_tab[iproc].nb_elt*mesh->num_nodes_hexa);
676 
677  ibuf_quadp_neighbor_hexa = MALLOC(int, proc_tab[iproc].nb_quad_p);
678  ibuf_quadp_neighbor_hexaface = MALLOC(int, proc_tab[iproc].nb_quad_p);
679 
680  i_buf = MALLOC(int, 7*26*info->proc[iproc].nb_elt); // 1 type + 5 params
681  cpu_neighbor_info = MALLOC(int, 3*26*info->proc[iproc].nb_ext);
682  }
683 
684  iparax = 0;
685  ifsurf = 0;
686 
687  for (ielt=0; ielt<proc_tab[iproc].nb_elt; ielt++) {// for each hexa of proc iproc
688 
689  int ihexa = ielt+1;
690  int ielt_num, ielt_number, ielt_face, ielt_edge, ielt_corner, ielt_coty, ielt_cpu;
691 
692  elt_info_t* curhex = proc_tab[iproc].local_elts[ielt];
693  int globnum = curhex->globalnum;
694 
695  // materiau
696  if (iproc != 0) ibuf_hexa_material_number[ielt] = mesh->layer[globnum];
697  else ig_hexa_material_number[ielt] = mesh->layer[globnum];
698 
699  // geometry
700  for (i=0; i<mesh->num_nodes_hexa; i++) {// geom coords of hexa geom nodes
701  // use efispec node order
702  // ig_hexa_gnode_glonum(inode,ihexa)
703  if (p_eind[p_eptr[globnum]+corner2efispec[i]-1] == 0) {
704  STOP("num node should not be 0 (fortran arrays start at 1)");
705  }
706  if (iproc != 0) ibuf_hexa_gnode_glonum[ielt * mesh->num_nodes_hexa + i] = p_eind[p_eptr[globnum]+corner2efispec[i]-1]; // numerotation starts from 1
707  else ig_hexa_gnode_glonum[ielt * mesh->num_nodes_hexa + i] = p_eind[p_eptr[globnum]+corner2efispec[i]-1]; // numerotation starts from 1
708  }
709  if (iproc == 0) mod_init_mesh_mp_init_gll_number_(&ihexa,ig_ngll_total);
710 
711  // faces
712  for (i=0; i<NFACE; i++) {
713  int iface = i+1;
714  // 6 faces connectivity
715  int type = curhex->faces[i].type;
716  switch (type) {
717  case HEXA : // hexa
718  ielt_num = elt_tab[curhex->faces[i].num_neigh].localnum + 1; // num commence à 1
719  ielt_face = face2efispec[curhex->faces[i].num_conn];
720  ielt_coty = curhex->faces[i].orientation;
721  ielt_cpu = mesh->part[curhex->faces[i].num_neigh];
722  if (ielt_cpu == iproc) {
723  if (iproc == 0) mod_init_mesh_mp_propagate_gll_nodes_face_(&ihexa, &iface, &ielt_num, &ielt_face, &ielt_coty);
724  else {
725  i_buf[nb_items*7] = FACE*10+HEXA;
726  i_buf[nb_items*7+1] = ihexa;
727  i_buf[nb_items*7+2] = iface;
728  i_buf[nb_items*7+3] = ielt_num;
729  i_buf[nb_items*7+4] = ielt_face;
730  i_buf[nb_items*7+5] = ielt_coty;
731  i_buf[nb_items*7+6] = FACE;
732  nb_items++;
733  }
734  } else {
735  if (iproc == 0) {
736  if (*ig_nneighbor_all_kind >= 26*info->proc[iproc].nb_ext) STOP("ig_nneighbor_all_kind too large\n");
737  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 0] = ihexa;
738  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 1] = iface;
739  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 2] = ielt_cpu;
740  (*ig_nneighbor_all_kind)++;
741  } else {
742  if (nneighbor_all_kind >= 26*info->proc[iproc].nb_ext) STOP("nneighbor_all_kind too large\n");
743  cpu_neighbor_info[nneighbor_all_kind*3 + 0] = ihexa;
744  cpu_neighbor_info[nneighbor_all_kind*3 + 1] = iface;
745  cpu_neighbor_info[nneighbor_all_kind*3 + 2] = ielt_cpu;
746  nneighbor_all_kind++;
747  }
748  }
749  break;
750  case QUAD_P : // quad p
751  if (iproc == 0) {
752  ig_quadp_neighbor_hexa[iparax] = ihexa;
753  ig_quadp_neighbor_hexaface[iparax] = iface;
754  } else {
755  ibuf_quadp_neighbor_hexa[iparax] = ihexa;
756  ibuf_quadp_neighbor_hexaface[iparax] = iface;
757  }
758  if (iproc == 0) {
759  if (mesh->num_node_per_dim*mesh->num_node_per_dim == 4) {
760  switch(iface) {
761  case 1 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
762  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
763  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
764  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
765  break;
766  case 2 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
767  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
768  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
769  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
770  break;
771  case 3 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
772  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
773  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
774  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
775  break;
776  case 4 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
777  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
778  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
779  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
780  break;
781  case 5 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
782  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
783  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
784  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
785  break;
786  case 6 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
787  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
788  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
789  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
790  break;
791  }
792  } else {
793  STOP("not yet ready for 27 nodes elements\n");
794  }
795  }
796  iparax++;
798  if (iproc == 0) mod_init_mesh_mp_propagate_gll_nodes_quad_(&ihexa, &iface, &iparax, ig_quadp_gll_glonum, &info->proc[iproc].nb_quad_p);
799  else {
800  i_buf[nb_items*7] = FACE*10+QUAD_P;
801  i_buf[nb_items*7+1] = ihexa;
802  i_buf[nb_items*7+2] = iface;
803  i_buf[nb_items*7+3] = iparax;
804  i_buf[nb_items*7+4] = 0;
805  i_buf[nb_items*7+5] = info->proc[iproc].nb_quad_p;
806  i_buf[nb_items*7+6] = FACE;
807  nb_items++;
808  // à la construction de ig_quadp_gnode_glonum pour les autres proc, il faudra prendre ihexa-1 et iparax-1 !!!
809  }
810  break;
811  case QUAD_F : // quad f
812  if (iproc == 0) {
813  if (mesh->num_node_per_dim*mesh->num_node_per_dim == 4) {
814  switch(iface) {
815  case 1 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
816  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
817  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
818  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
819  break;
820  case 2 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
821  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
822  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
823  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
824  break;
825  case 3 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+1];
826  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
827  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
828  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
829  break;
830  case 4 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
831  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+2];
832  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
833  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
834  break;
835  case 5 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+0];
836  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+3];
837  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
838  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
839  break;
840  case 6 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+4];
841  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+7];
842  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+6];
843  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*mesh->num_nodes_hexa+5];
844  break;
845  }
846  } else {
847  STOP("not yet ready for 27 nodes elements\n");
848  }
849  }
850  ifsurf++;
852  if (iproc == 0) mod_init_mesh_mp_propagate_gll_nodes_quad_(&ihexa, &iface, &ifsurf, ig_quadf_gll_glonum, &info->proc->nb_quad_f);
853  else {
854  i_buf[nb_items*7] = FACE*10+QUAD_F;
855  i_buf[nb_items*7+1] = ihexa;
856  i_buf[nb_items*7+2] = iface;
857  i_buf[nb_items*7+3] = ifsurf;
858  i_buf[nb_items*7+4] = 0;
859  i_buf[nb_items*7+5] = info->proc[iproc].nb_quad_f;
860  i_buf[nb_items*7+6] = FACE;
861  nb_items++;
862  // à la construction de ig_quadf_gnode_glonum pour les autres proc, il faudra prendre ihexa-1 et ifsurf-1 !!!
863  }
864  break;
865  }
866  }
867  // edges
868  for (i=0; i<NEDGE; i++) { // 12 edges connectivity
869  int iedge = i+1;
870  conn_info* ptr = &curhex->edges[i];
871  int nb_neighbor = getnbneigh(EDGE, curhex, i);
872  for (int ineigh = 0; ineigh < nb_neighbor; ineigh++) {
873  switch(ptr->type) {
874  case HEXA : ielt_number = elt_tab[ptr->num_neigh].localnum + 1; // start from 1
875  ielt_edge = edge2efispec[ptr->num_conn]; // num edge of neighbor
876  ielt_coty = ptr->orientation; // orientation
877  ielt_cpu = mesh->part[ptr->num_neigh]; // proc of neighbour
878  if (ielt_cpu == iproc) {
879  if (ielt_number > ihexa) {
880  if (iproc == 0 ) {
881  mod_init_mesh_mp_propagate_gll_nodes_edge_(&ihexa, &iedge, &ielt_number, &ielt_edge, &ielt_coty);
882  } else {
883  i_buf[nb_items*7] = EDGE*10+HEXA;
884  i_buf[nb_items*7+1] = ihexa;
885  i_buf[nb_items*7+2] = iedge;
886  i_buf[nb_items*7+3] = ielt_number;
887  i_buf[nb_items*7+4] = ielt_edge;
888  i_buf[nb_items*7+5] = ielt_coty;
889  i_buf[nb_items*7+6] = nb_neighbor;
890  nb_items++;
891  }
892  }
893  } else {
894  if (iproc == 0) {
895  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 0] = ihexa;
896  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 1] = NFACE + iedge;
897  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 2] = ielt_cpu;
898  (*ig_nneighbor_all_kind)++;
899  } else {
900  cpu_neighbor_info[nneighbor_all_kind*3 + 0] = ihexa;
901  cpu_neighbor_info[nneighbor_all_kind*3 + 1] = NFACE + iedge;
902  cpu_neighbor_info[nneighbor_all_kind*3 + 2] = ielt_cpu;
903  nneighbor_all_kind++;
904  }
905  }
906  case NONE : break;
907  default : STOP("edge neighbor should be an hexa");
908  }
909  ptr = ptr->next;
910  }
911  }
912  // corner
913  for (j=0; j<NCORNER; j++) { // 8 corners connectivity
914  int inode = j+1;
915  i = cornerefispec2cubit[j]; // written in efispec order
916  conn_info* ptr = &curhex->corners[i];
917  int nb_neighbor = getnbneigh(CORNER, curhex, i); // nb of neighbors connected to corner i
918  for (k=0; k<nb_neighbor; k++) {
919  switch(ptr->type) {
920  case HEXA : ielt_number = elt_tab[ptr->num_neigh].localnum + 1; // start from 1
921  ielt_corner = corner2efispec[ptr->num_conn]; // num edge of neighbor
922  ielt_cpu = mesh->part[ptr->num_neigh]; // proc of neighbour
923  if (ielt_cpu == iproc) {
924  if (ielt_number > ihexa) {
925  if (iproc == 0 ) {
926  mod_init_mesh_mp_propagate_gll_nodes_corner_(&ihexa, &inode, &ielt_number, &ielt_corner);
927  } else {
928  i_buf[nb_items*7] = CORNER*10+HEXA;
929  i_buf[nb_items*7+1] = ihexa;
930  i_buf[nb_items*7+2] = inode;
931  i_buf[nb_items*7+3] = ielt_number;
932  i_buf[nb_items*7+4] = ielt_corner;
933  i_buf[nb_items*7+5] = 0;
934  i_buf[nb_items*7+6] = nb_neighbor;
935  nb_items++;
936  }
937  }
938  } else {
939  if (iproc == 0) {
940  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 0] = ihexa;
941  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 1] = NFACE + NEDGE + inode;
942  ig_cpu_neighbor_info[*ig_nneighbor_all_kind*3 + 2] = ielt_cpu;
943  (*ig_nneighbor_all_kind)++;
944  } else {
945  cpu_neighbor_info[nneighbor_all_kind*3 + 0] = ihexa;
946  cpu_neighbor_info[nneighbor_all_kind*3 + 1] = NFACE + NEDGE + inode;
947  cpu_neighbor_info[nneighbor_all_kind*3 + 2] = ielt_cpu;
948  nneighbor_all_kind++;
949  }
950  }
951  case NONE : break;
952  default : STOP("corner neighbor should be an hexa");
953  }
954  ptr = ptr->next;
955  }
956  }
957 
958  }
959  if (iparax != info->proc[iproc].nb_quad_p || ifsurf != info->proc[iproc].nb_quad_f) STOP("problem with quad number");
960 
961  if (iproc != 0) {
962  MPI_Send(ibuf_hexa_material_number, proc_tab[iproc].nb_elt, MPI_INT, iproc, 21, MPI_COMM_WORLD);
963  MPI_Send(ibuf_hexa_gnode_glonum, proc_tab[iproc].nb_elt*mesh->num_nodes_hexa, MPI_INT, iproc, 22, MPI_COMM_WORLD);
964  MPI_Send(ibuf_quadp_neighbor_hexa, proc_tab[iproc].nb_quad_p, MPI_INT, iproc, 23, MPI_COMM_WORLD);
965  MPI_Send(ibuf_quadp_neighbor_hexaface, proc_tab[iproc].nb_quad_p, MPI_INT, iproc, 24, MPI_COMM_WORLD);
966 
967  ibuf2[0] = nb_items;
968  ibuf2[1] = nneighbor_all_kind;
969  MPI_Send(ibuf2, 2, MPI_INT, iproc, 25, MPI_COMM_WORLD);
970 
971  MPI_Send(i_buf, nb_items*7, MPI_INT, iproc, 26, MPI_COMM_WORLD);
972  MPI_Send(cpu_neighbor_info, 3*nneighbor_all_kind, MPI_INT, iproc, 27, MPI_COMM_WORLD);
973 
974  free(ibuf_hexa_material_number);
975  free(ibuf_hexa_gnode_glonum);
976  free(ibuf_quadp_neighbor_hexa);
977  free(ibuf_quadp_neighbor_hexaface);
978  free(i_buf);
979  free(cpu_neighbor_info);
980  }
981  } else {
982  if (rank == iproc) {
983  i_buf = MALLOC(int, 6*26*nhexa_local);
984  MPI_Recv(ig_hexa_material_number, nhexa_local, MPI_INT, 0, 21, MPI_COMM_WORLD, &status);
985  MPI_Recv(ig_hexa_gnode_glonum, nhexa_local*num_nodes_hexa, MPI_INT, 0, 22, MPI_COMM_WORLD, &status);
986  MPI_Recv(ig_quadp_neighbor_hexa, nquad_parax, MPI_INT, 0, 23, MPI_COMM_WORLD, &status);
987  MPI_Recv(ig_quadp_neighbor_hexaface, nquad_parax, MPI_INT, 0, 24, MPI_COMM_WORLD, &status);
988 
989  MPI_Recv(ibuf2, 2, MPI_INT, 0, 25, MPI_COMM_WORLD, &status);
990  nb_items = ibuf2[0];
991  *ig_nneighbor_all_kind = nneighbor_all_kind = ibuf2[1];
992 
993  MPI_Recv(i_buf, nb_items*7, MPI_INT, 0, 26, MPI_COMM_WORLD, &status);
994  MPI_Recv(ig_cpu_neighbor_info, 3*nneighbor_all_kind, MPI_INT, 0, 27, MPI_COMM_WORLD, &status);
995  int iibuf = 0;
996  for (ielt=0; ielt < nhexa_local; ielt++) {
997 
998  int ihexa = ielt+1;
999  mod_init_mesh_mp_init_gll_number_(&ihexa,ig_ngll_total);
1000  // besoin de savoir combien d'items -> tests avec ihexa & iface
1001  for (i=0; i<NFACE; i++) { // a revérifier ...
1002  int iface = i+1;
1003  if (i_buf[iibuf*7+6] == FACE && i_buf[iibuf*7+1] == ihexa && i_buf[iibuf*7+2] == iface) {
1004  int type = i_buf[iibuf*7];
1005  switch (type) {
1006  case FACE*10+HEXA : mod_init_mesh_mp_propagate_gll_nodes_face_(&i_buf[iibuf*7+1], &i_buf[iibuf*7+2], &i_buf[iibuf*7+3], &i_buf[iibuf*7+4], &i_buf[iibuf*7+5]);
1007  iibuf++;
1008  break;
1009  case FACE*10+QUAD_P : iparax = i_buf[iibuf*7+3]-1;
1010  assert(i_buf[iibuf*7+1]-1 == ielt);
1011  assert(i_buf[iibuf*7+2] == iface);
1012  switch(iface) {
1013  case 1 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1014  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1015  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1016  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1017  break;
1018  case 2 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1019  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1020  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1021  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1022  break;
1023  case 3 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1024  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1025  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1026  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1027  break;
1028  case 4 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1029  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1030  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1031  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1032  break;
1033  case 5 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1034  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1035  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1036  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1037  break;
1038  case 6 : ig_quadp_gnode_glonum[iparax*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1039  ig_quadp_gnode_glonum[iparax*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1040  ig_quadp_gnode_glonum[iparax*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1041  ig_quadp_gnode_glonum[iparax*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1042  break;
1043  }
1044  mod_init_mesh_mp_propagate_gll_nodes_quad_(&i_buf[iibuf*7+1], &i_buf[iibuf*7+2], &i_buf[iibuf*7+3], ig_quadp_gll_glonum, &i_buf[iibuf*7+5]);
1045  iibuf++;
1046  break;
1047  case FACE*10+QUAD_F : ifsurf = i_buf[iibuf*7+3]-1;
1048  assert(i_buf[iibuf*7+1]-1 == ielt);
1049  assert(i_buf[iibuf*7+2] == iface);
1050  switch(iface) {
1051  case 1 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1052  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1053  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1054  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1055  break;
1056  case 2 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1057  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1058  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1059  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1060  break;
1061  case 3 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+1];
1062  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1063  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1064  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1065  break;
1066  case 4 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1067  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+2];
1068  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1069  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1070  break;
1071  case 5 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+0];
1072  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+3];
1073  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1074  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1075  break;
1076  case 6 : ig_quadf_gnode_glonum[ifsurf*4+0] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+4];
1077  ig_quadf_gnode_glonum[ifsurf*4+1] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+7];
1078  ig_quadf_gnode_glonum[ifsurf*4+2] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+6];
1079  ig_quadf_gnode_glonum[ifsurf*4+3] = ig_hexa_gnode_glonum[ielt*num_nodes_hexa+5];
1080  break;
1081  }
1082  mod_init_mesh_mp_propagate_gll_nodes_quad_(&i_buf[iibuf*7+1], &i_buf[iibuf*7+2], &i_buf[iibuf*7+3], ig_quadf_gll_glonum, &i_buf[iibuf*7+5]);
1083  iibuf++;
1084  break;
1085  }
1086  }
1087  }
1088  for (i=0; i<NEDGE; i++) {
1089  int iedge = i+1;
1090  if (i_buf[iibuf*7] == EDGE*10+HEXA && i_buf[iibuf*7+1] == ihexa && i_buf[iibuf*7+2] == iedge) {
1091  int nb_neighbor = i_buf[iibuf*7+6];
1092  for (int ineigh = 0; ineigh < nb_neighbor; ineigh++) {
1093  if (i_buf[iibuf*7] == EDGE*10+HEXA && i_buf[iibuf*7+1] == ihexa && i_buf[iibuf*7+2] == iedge) {
1094  mod_init_mesh_mp_propagate_gll_nodes_edge_(&i_buf[iibuf*7+1], &i_buf[iibuf*7+2], &i_buf[iibuf*7+3], &i_buf[iibuf*7+4], &i_buf[iibuf*7+5]);
1095  iibuf++;
1096  }
1097  }
1098  }
1099  }
1100  for (j=0; j<NCORNER; j++) {
1101  int inode = j+1;
1102  if (i_buf[iibuf*7] == CORNER*10+HEXA && i_buf[iibuf*7+1] == ihexa && i_buf[iibuf*7+2] == inode) {
1103  int nb_neighbor = i_buf[iibuf*7+6];
1104  for (int ineigh = 0; ineigh < nb_neighbor; ineigh++) {
1105  if (i_buf[iibuf*7] == CORNER*10+HEXA && i_buf[iibuf*7+1] == ihexa && i_buf[iibuf*7+2] == inode) {
1106  mod_init_mesh_mp_propagate_gll_nodes_corner_(&i_buf[iibuf*7+1], &i_buf[iibuf*7+2], &i_buf[iibuf*7+3], &i_buf[iibuf*7+4]);
1107  iibuf++;
1108  }
1109  }
1110  }
1111  }
1112  }
1113  free(i_buf);
1114  assert(iibuf == nb_items);
1115  }
1116  }
1117  } // endfor iproc
1118 
1119  return;
1120 }
1121 
1122 
1123 // get all topological informations on all elements. return a filled info_t struct
1124 // write geom nodes in files and overwrite mesh->eind with local node num
1125 info_t* getConnInfo(mesh_t* mesh, int rank, int* ig_mesh_nnode)
1126 {
1127  int inode, ielt, iproc, inodeptr;
1128  idx_t* nodeMask;
1129  MPI_Status status;
1130  int nb_local_nodes;
1131  proc_info_t* proc_tab;
1132  elt_info_t * elt_tab;
1133  info_t * info;
1134 
1135  if (rank == 0) {
1136  // data structure allocation
1137  proc_tab = MALLOC(proc_info_t, mesh->npart); // DAVID2FREE
1138  elt_tab = MALLOC(elt_info_t, mesh->nh); // DAVID2FREE
1139  info = MALLOC(info_t, 1); // DAVID2FREE
1140  info->proc = proc_tab;
1141  info->elt = elt_tab;
1142 
1143  for(iproc = 0; iproc<mesh->npart; iproc++) {
1144  // printf("proc num %d\n", iproc);
1145  proc_tab[iproc].connex = MALLOC(char, mesh->npart); // DAVID2FREE
1146  memset(proc_tab[iproc].connex, 0, mesh->npart*sizeof(char));
1147  proc_tab[iproc].nb_elt = 0;
1148  proc_tab[iproc].nb_conn = 0;
1149  proc_tab[iproc].nb_ext = 0;
1150  proc_tab[iproc].nb_quad_p = 0;
1151  proc_tab[iproc].nb_quad_f = 0;
1152  }
1153 
1154  for(int i=0; i<mesh->nh; i++) {
1155  int j;
1156  for (j=0; j<NFACE; j++) {
1157  elt_tab[i].faces[j].defined = 0;
1158  }
1159  for (j=0; j<NEDGE; j++) {
1160  elt_tab[i].edges[j].defined = 0;
1161  }
1162  for (j=0; j<NCORNER; j++) {
1163  elt_tab[i].corners[j].defined = 0;
1164  }
1165  }
1166 
1167  // get all connectivity info for each element
1168  for(ielt=0; ielt<mesh->nh; ielt++) {
1169  //printf(" \relt : %d / %d", ielt, mesh->nh);
1170  iproc = mesh->part[ielt];
1171  elt_tab[ielt].globalnum = ielt;
1172 
1173  // check if elem is outer & set proc connectivity info
1174  elt_tab[ielt].outer = getProcConn(ielt, iproc, proc_tab, mesh);
1175 
1176  if (elt_tab[ielt].outer) proc_tab[iproc].nb_ext++;
1177 
1178  proc_tab[iproc].nb_elt++;
1179  // element's neighbors connectivity
1180  setEltConn(ielt, iproc, elt_tab, mesh, proc_tab);
1181  }
1182 
1183  // count local geom nodes
1184  nodeMask = MALLOC(idx_t, mesh->nn);
1185  }
1186 
1187  idx_t *p_eptr, *p_eind;
1188  if (mesh->hex27) {
1189  p_eptr = mesh->eptr_27;
1190  p_eind = mesh->eind_27;
1191  } else {
1192  p_eptr = mesh->eptr;
1193  p_eind = mesh->eind;
1194  }
1195  // à l'envers pour terminer avec iproc = 0
1196  float* buff_coord;
1197  for (iproc = mesh->npart-1; iproc >= 0; iproc--) {
1198 
1199  if (rank == 0) {
1200  memset(nodeMask, 0, mesh->nn*sizeof(idx_t));
1201  nb_local_nodes = 0;
1202  for(ielt=0; ielt<mesh->nh; ielt++) {
1203  if (mesh->part[ielt] == iproc) {
1204  for(inodeptr = p_eptr[ielt]; inodeptr<p_eptr[ielt+1]; inodeptr++) {
1205  if (nodeMask[p_eind[inodeptr]] == 0) {
1206  ++nb_local_nodes;
1207  nodeMask[p_eind[inodeptr]]++;
1208  }
1209  }
1210  }
1211  }
1212  if (iproc != 0) {
1213  // send nb_local_nodes to iproc
1214  buff_coord = MALLOC(float, nb_local_nodes*3);
1215  MPI_Send(&nb_local_nodes, 1, MPI_INT, iproc, 0, MPI_COMM_WORLD);
1216  } else {
1217  *ig_mesh_nnode = nb_local_nodes;
1218  // allocate Fortran arrays rg_gnode_x, rg_gnode_y, rg_gnode_z
1219  // the call back function pass_address_back_() is then called to store addresses in xcoord, ycoord & zcoord.
1220  // allocate_node_coord_arrays_(&nb_local_nodes);
1221  allocate_node_coord_arrays_();
1222  }
1223  } else {
1224  if (rank == iproc) {
1225  // recv nb_local_nodes from proc 0
1226  MPI_Recv(ig_mesh_nnode, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
1227  buff_coord = MALLOC(float, *ig_mesh_nnode*3);
1228  allocate_node_coord_arrays_();
1229  }
1230  }
1231 
1232  if (rank == 0) {
1233  // fill arrays
1234  memset(nodeMask, 0, mesh->nn*sizeof(idx_t));
1235  nb_local_nodes = 0;
1236  for(ielt=0; ielt<mesh->nh; ielt++) {
1237  if (mesh->part[ielt] == iproc) {
1238  for(inodeptr = p_eptr[ielt]; inodeptr<p_eptr[ielt+1]; inodeptr++) {
1239  int globalnodenum = p_eind[inodeptr];
1240  if (nodeMask[globalnodenum] == 0) {
1241  nodeMask[globalnodenum] = ++nb_local_nodes;
1242  p_eind[inodeptr] = nb_local_nodes; // overwrite eind since we don't use it anymore
1243  // warning 0 or 1 C/FORTRAN convention
1244  if (iproc == 0) {
1245  xcoord[nb_local_nodes-1] = mesh->ncoords_x[globalnodenum];
1246  ycoord[nb_local_nodes-1] = mesh->ncoords_y[globalnodenum];
1247  zcoord[nb_local_nodes-1] = mesh->ncoords_z[globalnodenum];
1248  } else {
1249  buff_coord[(nb_local_nodes-1)*3] = mesh->ncoords_x[globalnodenum];
1250  buff_coord[(nb_local_nodes-1)*3+1] = mesh->ncoords_y[globalnodenum];
1251  buff_coord[(nb_local_nodes-1)*3+2] = mesh->ncoords_z[globalnodenum];
1252  }
1253  } else {
1254  int localnumber = nodeMask[globalnodenum];
1255  p_eind[inodeptr] = localnumber; // overwrite eind since we don't use it anymore
1256  }
1257  }
1258  }
1259  }
1260  if (iproc != 0) {
1261  MPI_Send(buff_coord, nb_local_nodes*3, MPI_FLOAT, iproc, 1, MPI_COMM_WORLD);
1262  free(buff_coord);
1263  }
1264  } else {
1265  if (rank == iproc) {
1266  MPI_Recv(buff_coord, *ig_mesh_nnode*3, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
1267  for (inode = 0; inode<*ig_mesh_nnode; inode++) {
1268  xcoord[inode] = buff_coord[inode*3];
1269  ycoord[inode] = buff_coord[inode*3+1];
1270  zcoord[inode] = buff_coord[inode*3+2];
1271  }
1272  free(buff_coord);
1273  }
1274  }
1275  }
1276 
1277  if (rank == 0) {
1278  free(nodeMask);
1279  // alloc & settings
1280  int *num_outer;
1281  int *num_inner;
1282  num_outer = MALLOC(int, mesh->npart);
1283  num_inner = MALLOC(int, mesh->npart);
1284  for (iproc=0; iproc<mesh->npart; iproc++) {
1285  num_outer[iproc] = 0;
1286  num_inner[iproc] = proc_tab[iproc].nb_ext;
1287  proc_tab[iproc].local_elts = MALLOC(elt_info_t*, proc_tab[iproc].nb_elt); // DAVID2FREE
1288  }
1289 
1290  // build per proc data struct
1291  for(ielt=0; ielt<mesh->nh; ielt++) {
1292  int elt_proc = mesh->part[ielt];
1293  if (elt_tab[ielt].outer) {
1294  proc_tab[elt_proc].local_elts[num_outer[elt_proc]] = &elt_tab[ielt];
1295  elt_tab[ielt].localnum = num_outer[elt_proc]++; // warning, localnum start from 0
1296  } else {
1297  proc_tab[elt_proc].local_elts[num_inner[elt_proc]] = &elt_tab[ielt];
1298  elt_tab[ielt].localnum = num_inner[elt_proc]++; // warning, localnum start from 0
1299  }
1300  }
1301  free(num_inner);
1302  free(num_outer);
1303  return info;
1304  }
1305  return NULL;
1306 }
1307 
1308 // find topological details on connection between one hexahedral elts and all its neighbors, fill structures
1309 void setEltConn(int ielt, int iproc, elt_info_t* elt_tab, mesh_t* mesh, proc_info_t* proc_tab)
1310 {
1311  int ineig, inode, inode_n, i;
1312  //int num_node_per_face = mesh->num_node_per_dim*mesh->num_node_per_dim;
1313  int num_node_per_face = 4; // DAVID : pour les hexas à 27 noeuds, on ne se sert que des 8 premiers noeuds géom pour etablir la connectivité int* connex_nodes; connex_nodes = MALLOC(int, (2*num_node_per_face)); for (i=0; i<6; i++) { elt_tab[ielt].faces[i].type = NONE; elt_tab[ielt].faces[i].next = NULL; elt_tab[ielt].faces[i].defined = 0; } for (i=0; i<12; i++) { elt_tab[ielt].edges[i].type = NONE; elt_tab[ielt].edges[i].next = NULL; elt_tab[ielt].edges[i].defined = 0; } for (i=0; i<8; i++) { elt_tab[ielt].corners[i].type = NONE; elt_tab[ielt].corners[i].next = NULL; elt_tab[ielt].corners[i].defined = 0; } // loop on hexa neighbors for (ineig = mesh->xadj_hex[ielt]; ineig < mesh->xadj_hex[ielt+1]; ineig++) { int neigh = mesh->adjncy_hex[ineig]; int nb_conn = 0; // get nodes shared by ielt and its neighbor for(inode = mesh->eptr[ielt]; inode < mesh->eptr[ielt+1]; inode++) { for(inode_n = mesh->eptr[neigh]; inode_n < mesh->eptr[neigh+1]; inode_n++) { if (mesh->eind[inode] == mesh->eind[inode_n]) { connex_nodes[nb_conn] = inode - mesh->eptr[ielt]; // contact nodes numbers for current elt connex_nodes[nb_conn+num_node_per_face] = inode_n - mesh->eptr[neigh]; // associated node numbers for its neighbor nb_conn++; } } } // get connection topology if (nb_conn) { setHexaConnDetail(ielt, neigh, nb_conn, connex_nodes, elt_tab, mesh); } else { STOP("neighbor without contact"); } } // loop on quad neighbors for (ineig = mesh->xadj[ielt]; ineig < mesh->xadj[ielt+1]; ineig++) { int neigh = mesh->adjncy[ineig]; if (neigh >= mesh->nh) { // quad if (mesh->types[neigh] == QUAD_P) { proc_tab[iproc].nb_quad_p++; } else { proc_tab[iproc].nb_quad_f++; } int nb_conn = 0; // get nodes shared by ielt and its attached quad for(inode = mesh->eptr[ielt]; inode < mesh->eptr[ielt+1]; inode++) { for(inode_n = mesh->eptr[neigh]; inode_n < mesh->eptr[neigh+1]; inode_n++) { if (mesh->eind[inode] == mesh->eind[inode_n]) { connex_nodes[nb_conn] = inode - mesh->eptr[ielt]; // contact node number for current elt connex_nodes[nb_conn+num_node_per_face] = inode_n - mesh->eptr[neigh]; // associated node number for its neighbor nb_conn++; } } } if (nb_conn) { assert(nb_conn == N_NODES_QUAD); setQuadConnDetail(ielt, neigh, nb_conn, connex_nodes, elt_tab, mesh); } else { STOP("neighbor quad without contact"); } // TODO } } free(connex_nodes); return; } // find topological details on connection between a hexahedral elts and an attached quad, fill structures void setQuadConnDetail(int elt_num, int neigh_num, int nb_conn, int* connex_nodes, elt_info_t* elt_tab, mesh_t* mesh) { int face, inode, inode2; face = -1; for(inode = 0; inode<4; inode++) { if (connex_nodes[inode] == CORNER1) { for(inode2 = 0; inode2<4; inode2++) { switch(connex_nodes[inode2]) { case CORNER3 : face = FACE1; goto face_continue; case CORNER6 : face = FACE2; goto face_continue; case CORNER8 : face = FACE5; goto face_continue; } } } else if(connex_nodes[inode] == CORNER7) { for(inode2 = 0; inode2<4; inode2++) { switch(connex_nodes[inode2]) { case CORNER4 : face = FACE4; goto face_continue; case CORNER5 : face = FACE6; goto face_continue; case CORNER2 : face = FACE3; goto face_continue; } } } } face_continue: if (face == -1) { STOP("error in face search") } add_conn(FACE, elt_tab, elt_num, face, neigh_num, 0, 0, mesh->types[neigh_num]); return; } // find topological details on connection between 2 hexahedral elts, fill structures void setHexaConnDetail(int elt_num, int neigh_num, int nb_conn, int* connex_nodes, elt_info_t* elt_tab, mesh_t* mesh) { int icase, inode, inode2, num_conn_source, num_conn_target; switch(nb_conn) { case 1 : // 1 node connection (corner) num_conn_source = connex_nodes[0]; num_conn_target = connex_nodes[4]; add_conn(CORNER, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, 0, HEXA); break; case 2 : // 2 nodes connection (edge) for(icase = 0; icase<2; icase++) { int edge = -1; int offset = icase*4; for(inode = 0; inode<2; inode++) { switch(connex_nodes[inode+offset]) { case CORNER1 : for(inode2 = 0; inode2<2; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER2 : edge = EDGE1; goto edge_continue; case CORNER4 : edge = EDGE4; goto edge_continue; case CORNER5 : edge = EDGE9; goto edge_continue; } } case CORNER6 : for(inode2 = 0; inode2<2; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER2 : edge = EDGE10; goto edge_continue; case CORNER7 : edge = EDGE6; goto edge_continue; case CORNER5 : edge = EDGE5; goto edge_continue; } } case CORNER3 : for(inode2 = 0; inode2<2; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER2 : edge = EDGE2; goto edge_continue; case CORNER4 : edge = EDGE3; goto edge_continue; case CORNER7 : edge = EDGE11; goto edge_continue; } } case CORNER8 : for(inode2 = 0; inode2<2; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER7 : edge = EDGE7; goto edge_continue; case CORNER4 : edge = EDGE12; goto edge_continue; case CORNER5 : edge = EDGE8; goto edge_continue; } } } } edge_continue: if (edge == -1) { STOP("error in edge search") } if (icase==0) { num_conn_source = edge; } else { num_conn_target = edge; } } add_conn(EDGE, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, findEdgeOrientation(connex_nodes, num_conn_source, num_conn_target), HEXA); break; case 9 : case 4 : // 4 nodes connection (face) for(icase = 0; icase<2; icase++) { int face = -1; int offset = icase*4; for(inode = 0; inode<4; inode++) { if (connex_nodes[inode+offset] == CORNER1) { for(inode2 = 0; inode2<4; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER3 : face = FACE1; goto face_continue; case CORNER6 : face = FACE2; goto face_continue; case CORNER8 : face = FACE5; goto face_continue; } } } else if(connex_nodes[inode+offset] == CORNER7) { for(inode2 = 0; inode2<4; inode2++) { switch(connex_nodes[inode2+offset]) { case CORNER4 : face = FACE4; goto face_continue; case CORNER5 : face = FACE6; goto face_continue; case CORNER2 : face = FACE3; goto face_continue; } } } } face_continue: if (face == -1) { STOP("error in face search") } if (icase == 0) { num_conn_source= face; } else { num_conn_target= face; } } add_conn(FACE, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, findFaceOrientation(connex_nodes, num_conn_source, num_conn_target), HEXA); break; default : // problem ! STOP("Error in nb nodes connected") } return; } // return true if element is on the boundary // as a side effect, build procs adjacency graph int getProcConn(int ielt, int iproc, proc_info_t* proc_tab, mesh_t* mesh) { int ineig; int isOuter = 0; for (ineig = mesh->xadj_hex[ielt]; ineig < mesh->xadj_hex[ielt+1]; ineig++) { int neigh = mesh->adjncy_hex[ineig]; int neig_proc = mesh->part[neigh]; if (neig_proc != iproc) { if (!proc_tab[iproc].connex[neig_proc]) { proc_tab[iproc].connex[neig_proc] = 1; proc_tab[iproc].nb_conn++; } isOuter = 1; } } return isOuter; } // set elts weight for partionning void setHexaWeight(mesh_t* mesh) { idx_t ihex, ineigh; mesh->ncon = 1; mesh->vwgt = MALLOC(idx_t, mesh->nh); // DAVID2FREE for(ihex=0; ihex<mesh->nh; ihex++) { mesh->vwgt[ihex] = HEXA_WEIGHT; for (ineigh = mesh->xadj[ihex]; ineigh<mesh->xadj[ihex+1]; ineigh++) { idx_t neigh = mesh->adjncy[ineigh]; if(mesh->types[neigh] == QUAD_P) mesh->vwgt[ihex] += QUAD_P_WEIGHT; if(mesh->types[neigh] == QUAD_F) mesh->vwgt[ihex] += QUAD_F_WEIGHT; } } return; } void printMemUsage(mesh_t* mesh) { double sum = 0.; int i=1; char* unit; const long un = 1; sum += mesh->npart * sizeof(proc_info_t); sum += mesh->nh * sizeof(elt_info_t); sum += sizeof(info_t); sum += mesh->npart * mesh->npart * sizeof(char); sum += mesh->nh * sizeof(elt_info_t*); sum += (2*mesh->nh + mesh->ne+1 + mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad) * sizeof(idx_t); sum += mesh->nh * sizeof(unsigned char); sum += mesh->nn * 3 * sizeof(real_t); sum += mesh->ne * sizeof(elt_t); sum += mesh->nn * sizeof(idx_t); while(sum > un<<(i*10)) i++; switch(i) { case 1 : unit="Octets"; break; case 2 : unit="Kio"; break; case 3 : unit="Mio"; break; case 4 : unit="Gio"; break; case 5 : unit="Tio"; break; default: unit="Tio"; i=5; } printf("This program will use %d %s of memory\n", (int)(sum/(1<<((i-1)*10))),unit); } void readCubitMeshAbaqus(char* fileinp, mesh_t* mesh, int* number_of_elemnt_block) { char line[LENMAX]; int iblock = 0; FILE *unit_inp; char *hexa8=NULL, *hexa27=NULL; int readquad = 0; int ielt, inode, itmp; int loc_eind = 0; int ielt_glo = 0; char *quadf=NULL, *quadp=NULL; int ihexa=0; unit_inp = fopen(fileinp,"r"); if(unit_inp == NULL) { STOP("can't open input file (readCubitMeshAbaqus (1)") } //First pass on file *.inp: count number of nodes, hexa, quad_parax and quad_fsurf mesh->nh = mesh->nq_parax = mesh->nq_surf = 0; // 4 nodes quads in abaqus mesh->num_nodes_quad = 4; // no 27 nodes hexa in abaqus format ! mesh->hex27 = 0; while (!feof(unit_inp)) { char * cdum = fgets(line, LENMAX, unit_inp); //count number of geometric nodes if (strstr(line,"NSET=ALLNODES")) mesh->nn = count_elt(unit_inp); //count number of hexa if ((hexa8 = strstr(line,HEXA_8_HDR)) || (hexa27 = strstr(line,HEXA_27_HDR))) { // set number of node in hexa elts mesh->num_nodes_hexa = hexa8?8:27; // no 27 nodes hexa in abaqus format ! assert(mesh->num_nodes_hexa == 8); number_of_elemnt_block[iblock] = count_elt(unit_inp); mesh->nh += number_of_elemnt_block[iblock++]; if (readquad) { STOP("unexpected : quads are read before hexa in cubit mesh") } } //count number of quad_parax if (strstr(line,QUAD_P_HDR)) { number_of_elemnt_block[iblock] = count_elt(unit_inp); mesh->nq_parax += number_of_elemnt_block[iblock++]; readquad = 1; } //count number of quad_fsurf if (strstr(line,QUAD_F_HDR)) { number_of_elemnt_block[iblock] = count_elt(unit_inp); mesh->nq_surf += number_of_elemnt_block[iblock++]; readquad = 1; } } fclose(unit_inp); // allocation mesh->ne = mesh->nh + mesh->nq_parax + mesh->nq_surf; mesh->ncon = 0; mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE mesh->eind = MALLOC(idx_t, (mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad)); // DAVID2FREE mesh->layer = MALLOC(int, mesh->nh); // DAVID2FREE mesh->ncoords_x = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->ncoords_y = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->ncoords_z = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->types = MALLOC(elt_t, mesh->ne); // DAVID2FREE mesh->part = MALLOC(idx_t, mesh->nh); // DAVID2FREE mesh->num_node_per_dim = (int) cbrt(mesh->num_nodes_hexa); iblock = 0; unit_inp = fopen(fileinp,"r"); if(unit_inp == NULL) { STOP("can't open input file (readCubitMeshAbaqus (2)") } //Second pass on file *.inp: make eind & eptr iblock = ihexa = 0; while (!feof(unit_inp)) { char *cdum = fgets(line, LENMAX, unit_inp); //fill array gnodes with coordinates of geometric nodes if (strstr(line,"NSET=ALLNODES")) for (inode=0;inode<mesh->nn;inode++) { int idum = fscanf(unit_inp,"%d," "%"SCREAL "," "%"SCREAL "," "%"SCREAL, &itmp, &mesh->ncoords_x[inode], &mesh->ncoords_y[inode], &mesh->ncoords_z[inode]); //SCREAL is defined in metis.h } //fill eptr and eind for hexa if ((hexa8 = strstr(line,HEXA_8_HDR)) || (hexa27 = strstr(line,HEXA_27_HDR))) { // get num of layer char *hdr = hexa8?HEXA_8_HDR:HEXA_27_HDR; int layer = atoi(rtrim(strstr(line,hdr) + strlen(hdr))); for (ielt=0;ielt<number_of_elemnt_block[iblock];ielt++) { int idum = fscanf(unit_inp,"%d",&itmp); mesh->eptr[ielt_glo] = loc_eind; mesh->types[ielt_glo++] = HEXA; mesh->layer[ihexa++] = layer; for (inode=0;inode<mesh->num_nodes_hexa;inode++) { int idum = fscanf(unit_inp,",%d",&mesh->eind[loc_eind]); mesh->eind[loc_eind++]--; // because num starts on 1 in cubit meshes } } iblock++; } //fill eptr and eind for quad_parax and quad_fsurf if ((quadp = strstr(line,QUAD_P_HDR)) || (quadf = strstr(line,QUAD_F_HDR))) { for (ielt=0;ielt<number_of_elemnt_block[iblock];ielt++) { int idum = fscanf(unit_inp,"%d",&itmp); mesh->eptr[ielt_glo] = loc_eind; mesh->types[ielt_glo++] = quadp?QUAD_P:QUAD_F; for (inode=0;inode<mesh->num_nodes_quad;inode++) { int idum = fscanf(unit_inp,",%d",&mesh->eind[loc_eind]); mesh->eind[loc_eind++]--; // because num starts on 1 in cubit meshes } } iblock++; } } mesh->eptr[ielt_glo] = loc_eind; if (mesh->nh == 0) { STOP("not a valid mesh file"); } return; } void readCubitMeshExodus(char *fileinp, mesh_t *mesh, int *number_of_elemnt_block) { int cpu_word_size = sizeof(float); int io_word_size = 0; float version; int exoid; char title[MAX_LINE_LENGTH+1]; int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets; int res; int* idelbs; char elem_type[MAX_STR_LENGTH+1]; char block_name[MAX_STR_LENGTH+1]; int num_elem_this_blk, num_nodes_per_elem, num_attr; int* connect; mesh->nh = 0; mesh->nq_parax = 0; mesh->nq_surf = 0; // open file exoid = ex_open(fileinp, EX_READ, &cpu_word_size, &io_word_size, &version); // get header res = ex_get_init(exoid, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets); // get element blocks' ids idelbs = (int *) calloc(num_elem_blk, sizeof(int)); res = ex_get_elem_blk_ids (exoid, idelbs); // get number of element for each type/block for (int i=0; i<num_elem_blk; i++) { res = ex_get_elem_block (exoid, idelbs[i], elem_type, &num_elem_this_blk, &num_nodes_per_elem, &num_attr); if (!strcmp(elem_type, "HEX8") || !strcmp(elem_type, "HEX") || !strcmp(elem_type, "HEX27")) { mesh->nh += num_elem_this_blk; mesh->num_nodes_hexa = num_nodes_per_elem; mesh->hex27 = (num_nodes_per_elem == 27)?1:0; } else if (!strcmp(elem_type, "SHELL4") || !strcmp(elem_type, "QUAD4") || !strcmp(elem_type, "QUAD") || !strcmp(elem_type, "SHELL9")) { res = ex_get_name(exoid, EX_ELEM_BLOCK, idelbs[i], block_name); mesh->num_nodes_quad = num_nodes_per_elem; if (!strcmp(block_name, "prx")) { // parax mesh->nq_parax += num_elem_this_blk; } else if (!strcmp(block_name, "fsu")){ // surf mesh->nq_surf += num_elem_this_blk; } else { printf("ERROR : unknown quad type : %s\n",block_name); } } else { printf("block %d : ERROR, bad type : %s\n",i,elem_type); } } // allocations mesh->nn = num_nodes; mesh->ne = num_elem; mesh->ncon = 0; mesh->layer = MALLOC(int, mesh->nh); // DAVID2FREE mesh->ncoords_x = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->ncoords_y = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->ncoords_z = MALLOC(real_t, mesh->nn); // DAVID2FREE mesh->types = MALLOC(elt_t, mesh->ne); // DAVID2FREE mesh->part = MALLOC(idx_t, mesh->nh); // DAVID2FREE mesh->num_node_per_dim = (int) cbrt(mesh->num_nodes_hexa); if (mesh->hex27) { mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE mesh->eind = MALLOC(idx_t, (mesh->nh * 8 + (mesh->nq_parax + mesh->nq_surf) * 4)); // DAVID2FREE mesh->eptr_27 = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE mesh->eind_27 = MALLOC(idx_t, (mesh->nh * 27 + (mesh->nq_parax + mesh->nq_surf) * 9)); // DAVID2FREE } else { mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE mesh->eind = MALLOC(idx_t, (mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad)); // DAVID2FREE } // check assert(mesh->ne == mesh->nh + mesh->nq_parax + mesh->nq_surf); // get geometric nodes coordinates ex_get_coord(exoid, mesh->ncoords_x, mesh->ncoords_y, mesh->ncoords_z); // fill eptr and eind int ielem_glob = 0; int ihexa = 0; int offset; int global_offset = 0; int global_offset_27 = 0; for (int iblk=0; iblk<num_elem_blk; iblk++) { res = ex_get_elem_block (exoid, idelbs[iblk], elem_type, &num_elem_this_blk, &num_nodes_per_elem, &num_attr); res = ex_get_name(exoid, EX_ELEM_BLOCK, idelbs[iblk], block_name); // read element connectivity connect = (int *) calloc(num_nodes_per_elem*num_elem_this_blk, sizeof(int)); offset = 0; res = ex_get_elem_conn (exoid, idelbs[iblk], connect); if (!strcmp(elem_type, "HEX8") || !strcmp(elem_type, "HEX") || !strcmp(elem_type, "HEX27")) { int layer = atoi(&block_name[1]); // layers l01 l02 ... for (int ielem=0; ielem<num_elem_this_blk; ielem++) { mesh->eptr[ielem_glob] = global_offset; if (mesh->hex27) { mesh->eptr_27[ielem_glob] = global_offset_27; } mesh->types[ielem_glob] = HEXA; mesh->layer[ihexa++] = layer; for (int inode=0; inode<num_nodes_per_elem; inode++) { // DAVID : il faudra peut-être rajouter l'indirection ex2efi_hexa/quad pour 27/9 noeuds if (mesh->hex27) { mesh->eind_27[global_offset_27 + inode] = connect[offset + inode]-1; } //mesh->eind[global_offset+inode] = connect[offset + ex2efi_hexa[inode]]-1; if (inode<8) mesh->eind[global_offset + inode] = connect[offset + inode]-1; // because num starts on 1 in cubit meshes } offset += num_nodes_per_elem; if (mesh->hex27) { global_offset_27 += 27; global_offset += 8; } else { global_offset += num_nodes_per_elem; } ielem_glob++; } } else if (!strcmp(elem_type, "SHELL4") || !strcmp(elem_type, "QUAD4") || !strcmp(elem_type, "QUAD") || !strcmp(elem_type, "SHELL9")) { for (int ielem=0; ielem<num_elem_this_blk; ielem++) { mesh->eptr[ielem_glob] = global_offset; if (mesh->hex27) { mesh->eptr_27[ielem_glob] = global_offset_27; } if (!strcmp(block_name, "prx")) { mesh->types[ielem_glob] = QUAD_P; } else if (!strcmp(block_name, "fsu")){ mesh->types[ielem_glob] = QUAD_F; } else { printf("ERROR : unknown quad type : %s\n",block_name); } for (int inode=0; inode<num_nodes_per_elem; inode++) { // DAVID : il faudra peut-être rajouter l'indirection ex2efi_hexa/quad pour 27/9 noeuds if (mesh->hex27) { mesh->eind_27[global_offset_27 + inode] = connect[offset + inode]-1; } //mesh->eind[global_offset+inode] = connect[offset + ex2efi_quad[inode]]-1; if (inode<4) mesh->eind[global_offset+inode] = connect[offset + inode]-1; // because num starts on 1 in cubit meshes } offset += num_nodes_per_elem; if (mesh->hex27) { global_offset_27 += 9; global_offset += 4; } else { global_offset += num_nodes_per_elem; } ielem_glob++; } } else { printf("block %d : ERROR, bad type : %s\n",iblk,elem_type); } free(connect); } // dummy ptr on the end of array mesh->eptr[ielem_glob] = global_offset; if (mesh->hex27) mesh->eptr_27[ielem_glob] = global_offset_27; free(idelbs); ex_close(exoid); return; } char* rotateVect(char* ortn, char* ret_or, int cw_quarter_rot) { // ortn[0] -> horizontal vect // ortn[1] -> vertical vect // 1=i, 2=j, 3=k // sign : + : left->right or bottom->up // - : left<-right or bottom<-up // ie ortn[] = {1, -3} -> // +-----> i // | // | k // V char Hsign, Vsign; char H = ortn[0]; char V = ortn[1]; // -H(00) -> +V(11) -> +H(10) -> -V(01) clockwise rotation. H or V : right digit, sign : left digit // 0 3 2 1 char Hrot = H>0?2:0; char Vrot = V>0?3:1; int nbrot = cw_quarter_rot%4; Hrot = (4-nbrot+Hrot)&3; Hsign = (Hrot>>1)?1:-1; Vrot = (4-nbrot+Vrot)&3; Vsign = (Vrot>>1)?1:-1; if (Vrot&1) ret_or[1] = Vsign*abs(V); // impair -> Vert else ret_or[0] = Vsign*abs(V); // pair -> Hor if (Hrot&1) ret_or[1] = Hsign*abs(H); // impair -> Vert else ret_or[0] = Hsign*abs(H); // pair -> Hor return ret_or; } int findFaceOrientation(int *connex_nodes, int num_conn_source, int num_conn_target) { char or[2]; char or_rot[2]; int i, rot; for (i=0; i<4; i++) { if (connex_nodes[4] == face2faceRot[num_conn_source][num_conn_target][i]) { rot = i; break; } } // printf("node s : %d, node t : %d, face2faceRot[%d][%d][%d]=%d\n",connex_nodes[0], connex_nodes[4], num_conn_source+1, num_conn_target+1, i, face2faceRot[num_conn_source][num_conn_target][i]); memcpy(or, face_orientation[num_conn_target], 2*sizeof(char)); // horizontal flip (on voit la face depuis l'arriere) or[0] = -or[0]; // -rot : sens inverse parce que vu a l'envers rotateVect(or, or_rot, -rot); char hsign = or_rot[0]>0?1:0; char vsign = or_rot[1]>0?1:0; char h=abs(or_rot[0]); char v=abs(or_rot[1]); char valret = (v&3)|(vsign<<2)|((h&3)<<3)|(hsign<<5); // printf("rot : %d, face source %d, face target %d : h:%s%s, v:%s%s, " BYTETOBINARYPATTERN , rot, num_conn_source+1, num_conn_target+1, hsign?"+":"-",(abs(or_rot[0])==1)?"i":(abs(or_rot[0])==2)?"j":"k", vsign?"+":"-", (abs(or_rot[1])==1)?"i":(abs(or_rot[1])==2)?"j":"k", BYTETOBINARY(valret)); // printf("\n\n"); return (int) valret; } int findEdgeOrientation(int *connex_nodes, int num_edge_source, int num_edge_target) { int source_order=0; if (edgeOrder[num_edge_source][0] == connex_nodes[0]) source_order = 1; if (edgeOrder[num_edge_target][0] == connex_nodes[4]) { if (source_order) return 1; // same direction } else { if (!source_order) return 1; // same direction } return 0; // opposite directions } void add_conn(topo_t typecon, elt_info_t* elt_tab, int hexa_src, int numcon_src, idx_t num_neigh, int numcon_neigh, int orientation, elt_t type) { conn_info *ptr, *last; switch(typecon) { case FACE : if (elt_tab[hexa_src].faces[numcon_src].defined) { // malloc fprintf(stderr, "ERROR : face %d of elt %d is already connected\n", numcon_src, hexa_src); exit(1); } else { ptr = &elt_tab[hexa_src].faces[numcon_src]; } break; case EDGE : if (elt_tab[hexa_src].edges[numcon_src].defined) { // malloc for (ptr = &elt_tab[hexa_src].edges[numcon_src]; ptr; ptr = ptr->next) last = ptr; ptr = last->next = MALLOC(conn_info, 1); } else { ptr = &elt_tab[hexa_src].edges[numcon_src]; } break; case CORNER : if (elt_tab[hexa_src].corners[numcon_src].defined) { // malloc for (ptr = &elt_tab[hexa_src].corners[numcon_src]; ptr; ptr = ptr->next) last = ptr; ptr = last->next = MALLOC(conn_info, 1); } else { ptr = &elt_tab[hexa_src].corners[numcon_src]; } break; } ptr->num_neigh = num_neigh; ptr->num_conn = numcon_neigh; ptr->orientation = orientation; ptr->type = type; ptr->next = NULL; ptr->defined = 1; return; } int getnbneigh(topo_t typecon, elt_info_t* elt_tab, int numcon_src) { conn_info *ptr; int count=0; switch(typecon) { case FACE : printf("WARNING : funct getnbneigh() should not be called for faces\n"); count++; break; case EDGE : if (elt_tab->edges[numcon_src].defined) { for (ptr = &elt_tab->edges[numcon_src]; ptr; ptr = ptr->next) count++; } break; case CORNER : if (elt_tab->corners[numcon_src].defined) { for (ptr = &elt_tab->corners[numcon_src]; ptr; ptr = ptr->next) count++; } break; } return count; } //This function count the number of lines to compute the number of nodes, etc. int count_elt(FILE *unit) { const char EOL = '\n'; char c; int ne = 0; while(1) { c = getc(unit); if (c == EOL) ne++; else if (c == '*') break; } return ne; } char *ltrim(char *s) { while(isspace(*s)) s++; return s; } char *rtrim(char *s) { char* back = s + strlen(s); while(isspace(*--back)); *(back+1) = '\0'; return s; } char *trim(char *s) { return rtrim(ltrim(s)); }
1314  int* connex_nodes;
1315 
1316  connex_nodes = MALLOC(int, (2*num_node_per_face));
1317 
1318  for (i=0; i<6; i++) {
1319  elt_tab[ielt].faces[i].type = NONE;
1320  elt_tab[ielt].faces[i].next = NULL;
1321  elt_tab[ielt].faces[i].defined = 0;
1322  }
1323  for (i=0; i<12; i++) {
1324  elt_tab[ielt].edges[i].type = NONE;
1325  elt_tab[ielt].edges[i].next = NULL;
1326  elt_tab[ielt].edges[i].defined = 0;
1327  }
1328  for (i=0; i<8; i++) {
1329  elt_tab[ielt].corners[i].type = NONE;
1330  elt_tab[ielt].corners[i].next = NULL;
1331  elt_tab[ielt].corners[i].defined = 0;
1332  }
1333 
1334  // loop on hexa neighbors
1335  for (ineig = mesh->xadj_hex[ielt]; ineig < mesh->xadj_hex[ielt+1]; ineig++) {
1336  int neigh = mesh->adjncy_hex[ineig];
1337  int nb_conn = 0;
1338 
1339  // get nodes shared by ielt and its neighbor
1340  for(inode = mesh->eptr[ielt]; inode < mesh->eptr[ielt+1]; inode++) {
1341  for(inode_n = mesh->eptr[neigh]; inode_n < mesh->eptr[neigh+1]; inode_n++) {
1342  if (mesh->eind[inode] == mesh->eind[inode_n]) {
1343  connex_nodes[nb_conn] = inode - mesh->eptr[ielt]; // contact nodes numbers for current elt
1344  connex_nodes[nb_conn+num_node_per_face] = inode_n - mesh->eptr[neigh]; // associated node numbers for its neighbor
1345  nb_conn++;
1346  }
1347  }
1348  }
1349  // get connection topology
1350  if (nb_conn) {
1351  setHexaConnDetail(ielt, neigh, nb_conn, connex_nodes, elt_tab, mesh);
1352  } else {
1353  STOP("neighbor without contact");
1354  }
1355  }
1356 
1357  // loop on quad neighbors
1358  for (ineig = mesh->xadj[ielt]; ineig < mesh->xadj[ielt+1]; ineig++) {
1359  int neigh = mesh->adjncy[ineig];
1360  if (neigh >= mesh->nh) { // quad
1361  if (mesh->types[neigh] == QUAD_P) {
1362  proc_tab[iproc].nb_quad_p++;
1363  } else {
1364  proc_tab[iproc].nb_quad_f++;
1365  }
1366  int nb_conn = 0;
1367  // get nodes shared by ielt and its attached quad
1368  for(inode = mesh->eptr[ielt]; inode < mesh->eptr[ielt+1]; inode++) {
1369  for(inode_n = mesh->eptr[neigh]; inode_n < mesh->eptr[neigh+1]; inode_n++) {
1370  if (mesh->eind[inode] == mesh->eind[inode_n]) {
1371  connex_nodes[nb_conn] = inode - mesh->eptr[ielt]; // contact node number for current elt
1372  connex_nodes[nb_conn+num_node_per_face] = inode_n - mesh->eptr[neigh]; // associated node number for its neighbor
1373  nb_conn++;
1374  }
1375  }
1376  }
1377  if (nb_conn) {
1378  assert(nb_conn == N_NODES_QUAD);
1379  setQuadConnDetail(ielt, neigh, nb_conn, connex_nodes, elt_tab, mesh);
1380  } else {
1381  STOP("neighbor quad without contact");
1382  }
1383  // TODO
1384  }
1385  }
1386  free(connex_nodes);
1387  return;
1388 }
1389 
1390 // find topological details on connection between a hexahedral elts and an attached quad, fill structures
1391 void setQuadConnDetail(int elt_num, int neigh_num, int nb_conn, int* connex_nodes, elt_info_t* elt_tab, mesh_t* mesh)
1392 {
1393  int face, inode, inode2;
1394 
1395  face = -1;
1396  for(inode = 0; inode<4; inode++) {
1397  if (connex_nodes[inode] == CORNER1) {
1398  for(inode2 = 0; inode2<4; inode2++) {
1399  switch(connex_nodes[inode2]) {
1400  case CORNER3 : face = FACE1; goto face_continue;
1401  case CORNER6 : face = FACE2; goto face_continue;
1402  case CORNER8 : face = FACE5; goto face_continue;
1403  }
1404  }
1405  } else if(connex_nodes[inode] == CORNER7) {
1406  for(inode2 = 0; inode2<4; inode2++) {
1407  switch(connex_nodes[inode2]) {
1408  case CORNER4 : face = FACE4; goto face_continue;
1409  case CORNER5 : face = FACE6; goto face_continue;
1410  case CORNER2 : face = FACE3; goto face_continue;
1411  }
1412  }
1413  }
1414  }
1415 face_continue:
1416  if (face == -1) {
1417  STOP("error in face search")
1418  }
1419  add_conn(FACE, elt_tab, elt_num, face, neigh_num, 0, 0, mesh->types[neigh_num]);
1420  return;
1421 }
1422 
1423 // find topological details on connection between 2 hexahedral elts, fill structures
1424 void setHexaConnDetail(int elt_num, int neigh_num, int nb_conn, int* connex_nodes, elt_info_t* elt_tab, mesh_t* mesh)
1425 {
1426  int icase, inode, inode2, num_conn_source, num_conn_target;
1427 
1428  switch(nb_conn) {
1429  case 1 : // 1 node connection (corner)
1430  num_conn_source = connex_nodes[0];
1431  num_conn_target = connex_nodes[4];
1432  add_conn(CORNER, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, 0, HEXA);
1433  break;
1434  case 2 : // 2 nodes connection (edge)
1435  for(icase = 0; icase<2; icase++) {
1436  int edge = -1;
1437  int offset = icase*4;
1438  for(inode = 0; inode<2; inode++) {
1439  switch(connex_nodes[inode+offset]) {
1440  case CORNER1 :
1441  for(inode2 = 0; inode2<2; inode2++) {
1442  switch(connex_nodes[inode2+offset]) {
1443  case CORNER2 :
1444  edge = EDGE1; goto edge_continue;
1445  case CORNER4 :
1446  edge = EDGE4; goto edge_continue;
1447  case CORNER5 :
1448  edge = EDGE9; goto edge_continue;
1449  }
1450  }
1451  case CORNER6 :
1452  for(inode2 = 0; inode2<2; inode2++) {
1453  switch(connex_nodes[inode2+offset]) {
1454  case CORNER2 :
1455  edge = EDGE10; goto edge_continue;
1456  case CORNER7 :
1457  edge = EDGE6; goto edge_continue;
1458  case CORNER5 :
1459  edge = EDGE5; goto edge_continue;
1460  }
1461  }
1462  case CORNER3 :
1463  for(inode2 = 0; inode2<2; inode2++) {
1464  switch(connex_nodes[inode2+offset]) {
1465  case CORNER2 :
1466  edge = EDGE2; goto edge_continue;
1467  case CORNER4 :
1468  edge = EDGE3; goto edge_continue;
1469  case CORNER7 :
1470  edge = EDGE11; goto edge_continue;
1471  }
1472  }
1473  case CORNER8 :
1474  for(inode2 = 0; inode2<2; inode2++) {
1475  switch(connex_nodes[inode2+offset]) {
1476  case CORNER7 :
1477  edge = EDGE7; goto edge_continue;
1478  case CORNER4 :
1479  edge = EDGE12; goto edge_continue;
1480  case CORNER5 :
1481  edge = EDGE8; goto edge_continue;
1482  }
1483  }
1484  }
1485  }
1486 edge_continue:
1487  if (edge == -1) {
1488  STOP("error in edge search")
1489  }
1490  if (icase==0) {
1491  num_conn_source = edge;
1492  } else {
1493  num_conn_target = edge;
1494  }
1495  }
1496  add_conn(EDGE, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, findEdgeOrientation(connex_nodes, num_conn_source, num_conn_target), HEXA);
1497  break;
1498  case 9 :
1499  case 4 : // 4 nodes connection (face)
1500  for(icase = 0; icase<2; icase++) {
1501  int face = -1;
1502  int offset = icase*4;
1503  for(inode = 0; inode<4; inode++) {
1504  if (connex_nodes[inode+offset] == CORNER1) {
1505  for(inode2 = 0; inode2<4; inode2++) {
1506  switch(connex_nodes[inode2+offset]) {
1507  case CORNER3 : face = FACE1; goto face_continue;
1508  case CORNER6 : face = FACE2; goto face_continue;
1509  case CORNER8 : face = FACE5; goto face_continue;
1510  }
1511  }
1512  } else if(connex_nodes[inode+offset] == CORNER7) {
1513  for(inode2 = 0; inode2<4; inode2++) {
1514  switch(connex_nodes[inode2+offset]) {
1515  case CORNER4 : face = FACE4; goto face_continue;
1516  case CORNER5 : face = FACE6; goto face_continue;
1517  case CORNER2 : face = FACE3; goto face_continue;
1518  }
1519  }
1520  }
1521  }
1522 face_continue:
1523  if (face == -1) {
1524  STOP("error in face search")
1525  }
1526  if (icase == 0) {
1527  num_conn_source= face;
1528  } else {
1529  num_conn_target= face;
1530  }
1531  }
1532  add_conn(FACE, elt_tab, elt_num, num_conn_source, neigh_num, num_conn_target, findFaceOrientation(connex_nodes, num_conn_source, num_conn_target), HEXA);
1533  break;
1534  default : // problem !
1535  STOP("Error in nb nodes connected")
1536  }
1537  return;
1538 }
1539 
1540 // return true if element is on the boundary
1541 // as a side effect, build procs adjacency graph
1542 int getProcConn(int ielt, int iproc, proc_info_t* proc_tab, mesh_t* mesh)
1543 {
1544  int ineig;
1545  int isOuter = 0;
1546 
1547  for (ineig = mesh->xadj_hex[ielt]; ineig < mesh->xadj_hex[ielt+1]; ineig++) {
1548  int neigh = mesh->adjncy_hex[ineig];
1549  int neig_proc = mesh->part[neigh];
1550  if (neig_proc != iproc) {
1551  if (!proc_tab[iproc].connex[neig_proc]) {
1552  proc_tab[iproc].connex[neig_proc] = 1;
1553  proc_tab[iproc].nb_conn++;
1554  }
1555  isOuter = 1;
1556  }
1557  }
1558  return isOuter;
1559 }
1560 
1561 // set elts weight for partionning
1562 void setHexaWeight(mesh_t* mesh)
1563 {
1564  idx_t ihex, ineigh;
1565 
1566  mesh->ncon = 1;
1567  mesh->vwgt = MALLOC(idx_t, mesh->nh); // DAVID2FREE
1568  for(ihex=0; ihex<mesh->nh; ihex++) {
1569  mesh->vwgt[ihex] = HEXA_WEIGHT;
1570  for (ineigh = mesh->xadj[ihex]; ineigh<mesh->xadj[ihex+1]; ineigh++) {
1571  idx_t neigh = mesh->adjncy[ineigh];
1572  if(mesh->types[neigh] == QUAD_P) mesh->vwgt[ihex] += QUAD_P_WEIGHT;
1573  if(mesh->types[neigh] == QUAD_F) mesh->vwgt[ihex] += QUAD_F_WEIGHT;
1574  }
1575  }
1576  return;
1577 }
1578 
1579 void printMemUsage(mesh_t* mesh) {
1580 
1581  double sum = 0.;
1582  int i=1;
1583  char* unit;
1584  const long un = 1;
1585 
1586  sum += mesh->npart * sizeof(proc_info_t);
1587  sum += mesh->nh * sizeof(elt_info_t);
1588  sum += sizeof(info_t);
1589  sum += mesh->npart * mesh->npart * sizeof(char);
1590  sum += mesh->nh * sizeof(elt_info_t*);
1591  sum += (2*mesh->nh + mesh->ne+1 + mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad) * sizeof(idx_t);
1592  sum += mesh->nh * sizeof(unsigned char);
1593  sum += mesh->nn * 3 * sizeof(real_t);
1594  sum += mesh->ne * sizeof(elt_t);
1595  sum += mesh->nn * sizeof(idx_t);
1596 
1597  while(sum > un<<(i*10)) i++;
1598  switch(i) {
1599  case 1 : unit="Octets"; break;
1600  case 2 : unit="Kio"; break;
1601  case 3 : unit="Mio"; break;
1602  case 4 : unit="Gio"; break;
1603  case 5 : unit="Tio"; break;
1604  default: unit="Tio"; i=5;
1605  }
1606  printf("This program will use %d %s of memory\n", (int)(sum/(1<<((i-1)*10))),unit);
1607 }
1608 
1609 
1610 void readCubitMeshAbaqus(char* fileinp, mesh_t* mesh, int* number_of_elemnt_block)
1611 {
1612  char line[LENMAX];
1613  int iblock = 0;
1614  FILE *unit_inp;
1615  char *hexa8=NULL, *hexa27=NULL;
1616  int readquad = 0;
1617  int ielt, inode, itmp;
1618  int loc_eind = 0;
1619  int ielt_glo = 0;
1620  char *quadf=NULL, *quadp=NULL;
1621  int ihexa=0;
1622 
1623  unit_inp = fopen(fileinp,"r");
1624  if(unit_inp == NULL)
1625  {
1626  STOP("can't open input file (readCubitMeshAbaqus (1)")
1627  }
1628 
1629  //First pass on file *.inp: count number of nodes, hexa, quad_parax and quad_fsurf
1630  mesh->nh = mesh->nq_parax = mesh->nq_surf = 0;
1631  // 4 nodes quads in abaqus
1632  mesh->num_nodes_quad = 4;
1633 
1634  // no 27 nodes hexa in abaqus format !
1635  mesh->hex27 = 0;
1636  while (!feof(unit_inp))
1637  {
1638  char * cdum = fgets(line, LENMAX, unit_inp);
1639 
1640  //count number of geometric nodes
1641  if (strstr(line,"NSET=ALLNODES")) mesh->nn = count_elt(unit_inp);
1642 
1643  //count number of hexa
1644  if ((hexa8 = strstr(line,HEXA_8_HDR)) || (hexa27 = strstr(line,HEXA_27_HDR))) {
1645  // set number of node in hexa elts
1646  mesh->num_nodes_hexa = hexa8?8:27;
1647  // no 27 nodes hexa in abaqus format !
1648  assert(mesh->num_nodes_hexa == 8);
1649  number_of_elemnt_block[iblock] = count_elt(unit_inp);
1650  mesh->nh += number_of_elemnt_block[iblock++];
1651  if (readquad) {
1652  STOP("unexpected : quads are read before hexa in cubit mesh")
1653  }
1654  }
1655 
1656 
1657  //count number of quad_parax
1658  if (strstr(line,QUAD_P_HDR)) {
1659  number_of_elemnt_block[iblock] = count_elt(unit_inp);
1660  mesh->nq_parax += number_of_elemnt_block[iblock++];
1661  readquad = 1;
1662  }
1663 
1664  //count number of quad_fsurf
1665  if (strstr(line,QUAD_F_HDR)) {
1666  number_of_elemnt_block[iblock] = count_elt(unit_inp);
1667  mesh->nq_surf += number_of_elemnt_block[iblock++];
1668  readquad = 1;
1669  }
1670 
1671  }
1672  fclose(unit_inp);
1673 
1674  // allocation
1675  mesh->ne = mesh->nh + mesh->nq_parax + mesh->nq_surf;
1676  mesh->ncon = 0;
1677  mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE
1678  mesh->eind = MALLOC(idx_t, (mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad)); // DAVID2FREE
1679  mesh->layer = MALLOC(int, mesh->nh); // DAVID2FREE
1680  mesh->ncoords_x = MALLOC(real_t, mesh->nn); // DAVID2FREE
1681  mesh->ncoords_y = MALLOC(real_t, mesh->nn); // DAVID2FREE
1682  mesh->ncoords_z = MALLOC(real_t, mesh->nn); // DAVID2FREE
1683  mesh->types = MALLOC(elt_t, mesh->ne); // DAVID2FREE
1684  mesh->part = MALLOC(idx_t, mesh->nh); // DAVID2FREE
1685  mesh->num_node_per_dim = (int) cbrt(mesh->num_nodes_hexa);
1686 
1687  iblock = 0;
1688  unit_inp = fopen(fileinp,"r");
1689  if(unit_inp == NULL)
1690  {
1691  STOP("can't open input file (readCubitMeshAbaqus (2)")
1692  }
1693 
1694  //Second pass on file *.inp: make eind & eptr
1695  iblock = ihexa = 0;
1696  while (!feof(unit_inp))
1697  {
1698  char *cdum = fgets(line, LENMAX, unit_inp);
1699 
1700  //fill array gnodes with coordinates of geometric nodes
1701  if (strstr(line,"NSET=ALLNODES"))
1702  for (inode=0;inode<mesh->nn;inode++) {
1703  int idum = fscanf(unit_inp,"%d," "%"SCREAL "," "%"SCREAL "," "%"SCREAL, &itmp, &mesh->ncoords_x[inode], &mesh->ncoords_y[inode], &mesh->ncoords_z[inode]); //SCREAL is defined in metis.h
1704  }
1705 
1706  //fill eptr and eind for hexa
1707  if ((hexa8 = strstr(line,HEXA_8_HDR)) || (hexa27 = strstr(line,HEXA_27_HDR)))
1708  {
1709  // get num of layer
1710  char *hdr = hexa8?HEXA_8_HDR:HEXA_27_HDR;
1711  int layer = atoi(rtrim(strstr(line,hdr) + strlen(hdr)));
1712 
1713  for (ielt=0;ielt<number_of_elemnt_block[iblock];ielt++)
1714  {
1715  int idum = fscanf(unit_inp,"%d",&itmp);
1716  mesh->eptr[ielt_glo] = loc_eind;
1717  mesh->types[ielt_glo++] = HEXA;
1718  mesh->layer[ihexa++] = layer;
1719  for (inode=0;inode<mesh->num_nodes_hexa;inode++)
1720  {
1721  int idum = fscanf(unit_inp,",%d",&mesh->eind[loc_eind]);
1722  mesh->eind[loc_eind++]--; // because num starts on 1 in cubit meshes
1723  }
1724  }
1725 
1726  iblock++;
1727  }
1728 
1729  //fill eptr and eind for quad_parax and quad_fsurf
1730  if ((quadp = strstr(line,QUAD_P_HDR)) || (quadf = strstr(line,QUAD_F_HDR)))
1731  {
1732  for (ielt=0;ielt<number_of_elemnt_block[iblock];ielt++)
1733  {
1734  int idum = fscanf(unit_inp,"%d",&itmp);
1735  mesh->eptr[ielt_glo] = loc_eind;
1736  mesh->types[ielt_glo++] = quadp?QUAD_P:QUAD_F;
1737  for (inode=0;inode<mesh->num_nodes_quad;inode++)
1738  {
1739  int idum = fscanf(unit_inp,",%d",&mesh->eind[loc_eind]);
1740  mesh->eind[loc_eind++]--; // because num starts on 1 in cubit meshes
1741  }
1742  }
1743  iblock++;
1744  }
1745  }
1746  mesh->eptr[ielt_glo] = loc_eind;
1747 
1748  if (mesh->nh == 0) {
1749  STOP("not a valid mesh file");
1750  }
1751  return;
1752 }
1753 
1754 void readCubitMeshExodus(char *fileinp, mesh_t *mesh, int *number_of_elemnt_block)
1755 {
1756  int cpu_word_size = sizeof(float);
1757  int io_word_size = 0;
1758  float version;
1759 
1760  int exoid;
1761  char title[MAX_LINE_LENGTH+1];
1762  int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets;
1763  int res;
1764 
1765  int* idelbs;
1766  char elem_type[MAX_STR_LENGTH+1];
1767  char block_name[MAX_STR_LENGTH+1];
1768  int num_elem_this_blk, num_nodes_per_elem, num_attr;
1769 
1770  int* connect;
1771 
1772  mesh->nh = 0;
1773  mesh->nq_parax = 0;
1774  mesh->nq_surf = 0;
1775 
1776  // open file
1777  exoid = ex_open(fileinp, EX_READ, &cpu_word_size, &io_word_size, &version);
1778 
1779  // get header
1780  res = ex_get_init(exoid, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets);
1781 
1782  // get element blocks' ids
1783  idelbs = (int *) calloc(num_elem_blk, sizeof(int));
1784  res = ex_get_elem_blk_ids (exoid, idelbs);
1785 
1786  // get number of element for each type/block
1787  for (int i=0; i<num_elem_blk; i++) {
1788  res = ex_get_elem_block (exoid, idelbs[i], elem_type, &num_elem_this_blk, &num_nodes_per_elem, &num_attr);
1789  if (!strcmp(elem_type, "HEX8") || !strcmp(elem_type, "HEX") || !strcmp(elem_type, "HEX27")) {
1790  mesh->nh += num_elem_this_blk;
1791  mesh->num_nodes_hexa = num_nodes_per_elem;
1792  mesh->hex27 = (num_nodes_per_elem == 27)?1:0;
1793  } else if (!strcmp(elem_type, "SHELL4") || !strcmp(elem_type, "QUAD4") || !strcmp(elem_type, "QUAD") || !strcmp(elem_type, "SHELL9")) {
1794  res = ex_get_name(exoid, EX_ELEM_BLOCK, idelbs[i], block_name);
1795  mesh->num_nodes_quad = num_nodes_per_elem;
1796  if (!strcmp(block_name, "prx")) {
1797  // parax
1798  mesh->nq_parax += num_elem_this_blk;
1799  } else if (!strcmp(block_name, "fsu")){
1800  // surf
1801  mesh->nq_surf += num_elem_this_blk;
1802  } else {
1803  printf("ERROR : unknown quad type : %s\n",block_name);
1804  }
1805  } else {
1806  printf("block %d : ERROR, bad type : %s\n",i,elem_type);
1807  }
1808 
1809  }
1810 
1811  // allocations
1812  mesh->nn = num_nodes;
1813  mesh->ne = num_elem;
1814  mesh->ncon = 0;
1815  mesh->layer = MALLOC(int, mesh->nh); // DAVID2FREE
1816  mesh->ncoords_x = MALLOC(real_t, mesh->nn); // DAVID2FREE
1817  mesh->ncoords_y = MALLOC(real_t, mesh->nn); // DAVID2FREE
1818  mesh->ncoords_z = MALLOC(real_t, mesh->nn); // DAVID2FREE
1819  mesh->types = MALLOC(elt_t, mesh->ne); // DAVID2FREE
1820  mesh->part = MALLOC(idx_t, mesh->nh); // DAVID2FREE
1821  mesh->num_node_per_dim = (int) cbrt(mesh->num_nodes_hexa);
1822 
1823  if (mesh->hex27) {
1824  mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE
1825  mesh->eind = MALLOC(idx_t, (mesh->nh * 8 + (mesh->nq_parax + mesh->nq_surf) * 4)); // DAVID2FREE
1826  mesh->eptr_27 = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE
1827  mesh->eind_27 = MALLOC(idx_t, (mesh->nh * 27 + (mesh->nq_parax + mesh->nq_surf) * 9)); // DAVID2FREE
1828  } else {
1829  mesh->eptr = MALLOC(idx_t, (mesh->ne+1)); // DAVID2FREE
1830  mesh->eind = MALLOC(idx_t, (mesh->nh * mesh->num_nodes_hexa + (mesh->nq_parax + mesh->nq_surf) * mesh->num_nodes_quad)); // DAVID2FREE
1831  }
1832 
1833  // check
1834  assert(mesh->ne == mesh->nh + mesh->nq_parax + mesh->nq_surf);
1835 
1836  // get geometric nodes coordinates
1837  ex_get_coord(exoid, mesh->ncoords_x, mesh->ncoords_y, mesh->ncoords_z);
1838 
1839  // fill eptr and eind
1840  int ielem_glob = 0;
1841  int ihexa = 0;
1842  int offset;
1843  int global_offset = 0;
1844  int global_offset_27 = 0;
1845  for (int iblk=0; iblk<num_elem_blk; iblk++) {
1846  res = ex_get_elem_block (exoid, idelbs[iblk], elem_type, &num_elem_this_blk, &num_nodes_per_elem, &num_attr);
1847  res = ex_get_name(exoid, EX_ELEM_BLOCK, idelbs[iblk], block_name);
1848 
1849  // read element connectivity
1850  connect = (int *) calloc(num_nodes_per_elem*num_elem_this_blk, sizeof(int));
1851  offset = 0;
1852  res = ex_get_elem_conn (exoid, idelbs[iblk], connect);
1853  if (!strcmp(elem_type, "HEX8") || !strcmp(elem_type, "HEX") || !strcmp(elem_type, "HEX27")) {
1854  int layer = atoi(&block_name[1]); // layers l01 l02 ...
1855  for (int ielem=0; ielem<num_elem_this_blk; ielem++) {
1856  mesh->eptr[ielem_glob] = global_offset;
1857  if (mesh->hex27) {
1858  mesh->eptr_27[ielem_glob] = global_offset_27;
1859  }
1860  mesh->types[ielem_glob] = HEXA;
1861  mesh->layer[ihexa++] = layer;
1862  for (int inode=0; inode<num_nodes_per_elem; inode++) {
1863  // DAVID : il faudra peut-être rajouter l'indirection ex2efi_hexa/quad pour 27/9 noeuds
1864  if (mesh->hex27) {
1865  mesh->eind_27[global_offset_27 + inode] = connect[offset + inode]-1;
1866  }
1867  //mesh->eind[global_offset+inode] = connect[offset + ex2efi_hexa[inode]]-1;
1868  if (inode<8) mesh->eind[global_offset + inode] = connect[offset + inode]-1; // because num starts on 1 in cubit meshes
1869  }
1870  offset += num_nodes_per_elem;
1871  if (mesh->hex27) {
1872  global_offset_27 += 27;
1873  global_offset += 8;
1874  } else {
1875  global_offset += num_nodes_per_elem;
1876  }
1877  ielem_glob++;
1878  }
1879  } else if (!strcmp(elem_type, "SHELL4") || !strcmp(elem_type, "QUAD4") || !strcmp(elem_type, "QUAD") || !strcmp(elem_type, "SHELL9")) {
1880  for (int ielem=0; ielem<num_elem_this_blk; ielem++) {
1881  mesh->eptr[ielem_glob] = global_offset;
1882  if (mesh->hex27) {
1883  mesh->eptr_27[ielem_glob] = global_offset_27;
1884  }
1885  if (!strcmp(block_name, "prx")) {
1886  mesh->types[ielem_glob] = QUAD_P;
1887  } else if (!strcmp(block_name, "fsu")){
1888  mesh->types[ielem_glob] = QUAD_F;
1889  } else {
1890  printf("ERROR : unknown quad type : %s\n",block_name);
1891  }
1892  for (int inode=0; inode<num_nodes_per_elem; inode++) {
1893  // DAVID : il faudra peut-être rajouter l'indirection ex2efi_hexa/quad pour 27/9 noeuds
1894  if (mesh->hex27) {
1895  mesh->eind_27[global_offset_27 + inode] = connect[offset + inode]-1;
1896  }
1897  //mesh->eind[global_offset+inode] = connect[offset + ex2efi_quad[inode]]-1;
1898  if (inode<4) mesh->eind[global_offset+inode] = connect[offset + inode]-1; // because num starts on 1 in cubit meshes
1899  }
1900  offset += num_nodes_per_elem;
1901  if (mesh->hex27) {
1902  global_offset_27 += 9;
1903  global_offset += 4;
1904  } else {
1905  global_offset += num_nodes_per_elem;
1906  }
1907  ielem_glob++;
1908  }
1909  } else {
1910  printf("block %d : ERROR, bad type : %s\n",iblk,elem_type);
1911  }
1912  free(connect);
1913  }
1914  // dummy ptr on the end of array
1915  mesh->eptr[ielem_glob] = global_offset;
1916  if (mesh->hex27) mesh->eptr_27[ielem_glob] = global_offset_27;
1917 
1918  free(idelbs);
1919  ex_close(exoid);
1920  return;
1921 }
1922 
1923 char* rotateVect(char* ortn, char* ret_or, int cw_quarter_rot)
1924 {
1925  // ortn[0] -> horizontal vect
1926  // ortn[1] -> vertical vect
1927  // 1=i, 2=j, 3=k
1928  // sign : + : left->right or bottom->up
1929  // - : left<-right or bottom<-up
1930  // ie ortn[] = {1, -3} ->
1931  // +-----> i
1932  // |
1933  // | k
1934  // V
1935 
1936  char Hsign, Vsign;
1937 
1938  char H = ortn[0];
1939  char V = ortn[1];
1940  // -H(00) -> +V(11) -> +H(10) -> -V(01) clockwise rotation. H or V : right digit, sign : left digit
1941  // 0 3 2 1
1942  char Hrot = H>0?2:0;
1943  char Vrot = V>0?3:1;
1944  int nbrot = cw_quarter_rot%4;
1945 
1946  Hrot = (4-nbrot+Hrot)&3;
1947  Hsign = (Hrot>>1)?1:-1;
1948  Vrot = (4-nbrot+Vrot)&3;
1949  Vsign = (Vrot>>1)?1:-1;
1950 
1951 
1952  if (Vrot&1) ret_or[1] = Vsign*abs(V); // impair -> Vert
1953  else ret_or[0] = Vsign*abs(V); // pair -> Hor
1954 
1955  if (Hrot&1) ret_or[1] = Hsign*abs(H); // impair -> Vert
1956  else ret_or[0] = Hsign*abs(H); // pair -> Hor
1957 
1958  return ret_or;
1959 }
1960 
1961 int findFaceOrientation(int *connex_nodes, int num_conn_source, int num_conn_target)
1962 {
1963  char or[2];
1964  char or_rot[2];
1965  int i, rot;
1966  for (i=0; i<4; i++) {
1967  if (connex_nodes[4] == face2faceRot[num_conn_source][num_conn_target][i]) {
1968  rot = i;
1969  break;
1970  }
1971  }
1972  // printf("node s : %d, node t : %d, face2faceRot[%d][%d][%d]=%d\n",connex_nodes[0], connex_nodes[4], num_conn_source+1, num_conn_target+1, i, face2faceRot[num_conn_source][num_conn_target][i]);
1973  memcpy(or, face_orientation[num_conn_target], 2*sizeof(char));
1974 
1975  // horizontal flip (on voit la face depuis l'arriere)
1976  or[0] = -or[0];
1977  // -rot : sens inverse parce que vu a l'envers
1978  rotateVect(or, or_rot, -rot);
1979  char hsign = or_rot[0]>0?1:0;
1980  char vsign = or_rot[1]>0?1:0;
1981  char h=abs(or_rot[0]);
1982  char v=abs(or_rot[1]);
1983  char valret = (v&3)|(vsign<<2)|((h&3)<<3)|(hsign<<5);
1984  // printf("rot : %d, face source %d, face target %d : h:%s%s, v:%s%s, " BYTETOBINARYPATTERN , rot, num_conn_source+1, num_conn_target+1, hsign?"+":"-",(abs(or_rot[0])==1)?"i":(abs(or_rot[0])==2)?"j":"k", vsign?"+":"-", (abs(or_rot[1])==1)?"i":(abs(or_rot[1])==2)?"j":"k", BYTETOBINARY(valret));
1985  // printf("\n\n");
1986 
1987  return (int) valret;
1988 }
1989 
1990 int findEdgeOrientation(int *connex_nodes, int num_edge_source, int num_edge_target)
1991 {
1992  int source_order=0;
1993  if (edgeOrder[num_edge_source][0] == connex_nodes[0]) source_order = 1;
1994 
1995  if (edgeOrder[num_edge_target][0] == connex_nodes[4]) {
1996  if (source_order) return 1; // same direction
1997  } else {
1998  if (!source_order) return 1; // same direction
1999  }
2000 
2001  return 0; // opposite directions
2002 }
2003 
2004 
2005 void add_conn(topo_t typecon, elt_info_t* elt_tab, int hexa_src, int numcon_src, idx_t num_neigh, int numcon_neigh, int orientation, elt_t type)
2006 {
2007  conn_info *ptr, *last;
2008 
2009  switch(typecon) {
2010  case FACE : if (elt_tab[hexa_src].faces[numcon_src].defined) { // malloc
2011  fprintf(stderr, "ERROR : face %d of elt %d is already connected\n", numcon_src, hexa_src);
2012  exit(1);
2013  } else {
2014  ptr = &elt_tab[hexa_src].faces[numcon_src];
2015  }
2016  break;
2017  case EDGE : if (elt_tab[hexa_src].edges[numcon_src].defined) { // malloc
2018  for (ptr = &elt_tab[hexa_src].edges[numcon_src]; ptr; ptr = ptr->next) last = ptr;
2019  ptr = last->next = MALLOC(conn_info, 1);
2020  } else {
2021  ptr = &elt_tab[hexa_src].edges[numcon_src];
2022  }
2023  break;
2024  case CORNER : if (elt_tab[hexa_src].corners[numcon_src].defined) { // malloc
2025  for (ptr = &elt_tab[hexa_src].corners[numcon_src]; ptr; ptr = ptr->next) last = ptr;
2026  ptr = last->next = MALLOC(conn_info, 1);
2027  } else {
2028  ptr = &elt_tab[hexa_src].corners[numcon_src];
2029  }
2030  break;
2031  }
2032  ptr->num_neigh = num_neigh;
2033  ptr->num_conn = numcon_neigh;
2034  ptr->orientation = orientation;
2035  ptr->type = type;
2036  ptr->next = NULL;
2037  ptr->defined = 1;
2038  return;
2039 }
2040 
2041 int getnbneigh(topo_t typecon, elt_info_t* elt_tab, int numcon_src)
2042 {
2043  conn_info *ptr;
2044  int count=0;
2045  switch(typecon) {
2046  case FACE : printf("WARNING : funct getnbneigh() should not be called for faces\n");
2047  count++;
2048  break;
2049  case EDGE : if (elt_tab->edges[numcon_src].defined) {
2050  for (ptr = &elt_tab->edges[numcon_src]; ptr; ptr = ptr->next) count++;
2051  }
2052  break;
2053  case CORNER : if (elt_tab->corners[numcon_src].defined) {
2054  for (ptr = &elt_tab->corners[numcon_src]; ptr; ptr = ptr->next) count++;
2055  }
2056  break;
2057  }
2058  return count;
2059 }
2060 
2061 //This function count the number of lines to compute the number of nodes, etc.
2062 int count_elt(FILE *unit)
2063 {
2064  const char EOL = '\n';
2065  char c;
2066  int ne = 0;
2067 
2068  while(1)
2069  {
2070  c = getc(unit);
2071  if (c == EOL)
2072  ne++;
2073  else if (c == '*')
2074  break;
2075  }
2076  return ne;
2077 }
2078 
2079 char *ltrim(char *s)
2080 {
2081  while(isspace(*s)) s++;
2082  return s;
2083 }
2084 
2085 char *rtrim(char *s)
2086 {
2087  char* back = s + strlen(s);
2088  while(isspace(*--back));
2089  *(back+1) = '\0';
2090  return s;
2091 }
2092 
2093 char *trim(char *s)
2094 {
2095  return rtrim(ltrim(s));
2096 }
fill_efispec_arrays
void fill_efispec_arrays(info_t *info, mesh_t *mesh, int rank, int *ig_ngll_total, int *ig_nneighbor_all_kind, int *cpu_neighbor, int *ig_cpu_neighbor_info, int *ig_hexa_gnode_glonum, int *ig_quadp_gnode_glonum, int *ig_quadf_gnode_glonum, int *ig_hexa_gll_glonum, int *ig_quadp_gll_glonum, int *ig_quadf_gll_glonum, int *ig_quadp_neighbor_hexa, int *ig_quadp_neighbor_hexaface, int *ig_hexa_material_number)
Definition: partCubitMesh.c:606