root/src/exo.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. exo

   1 /*     Calculix - A 3-dimensional finite element program                 */
   2 /*              Copyright (C) 1998-2014 Guido Dhondt                     */
   3 /*     This subroutine                                                   */
   4 /*              Copyright (C) 2013-2014 Peter A. Gustafson               */
   5 /*                                                                       */
   6 /*     This program is free software; you can redistribute it and/or     */
   7 /*     modify it under the terms of the GNU General Public License as    */
   8 /*     published by the Free Software Foundation(version 2);    */
   9 /*                                                                       */
  10 
  11 /*     This program is distributed in the hope that it will be useful,   */
  12 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  13 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  14 /*     GNU General Public License for more details.                      */
  15 
  16 /*     You should have received a copy of the GNU General Public License */
  17 /*     along with this program; if not, write to the Free Software       */
  18 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  19 
  20 #include <stdlib.h>
  21 #include <math.h>
  22 #include <stdio.h>
  23 #include <string.h>
  24 #include "CalculiX.h"
  25 #include "exodusII.h"
  26 #include "exo.h"
  27 
  28 #define min(a,b) ((a) <= (b) ? (a) : (b))
  29 #define max(a,b) ((a) >= (b) ? (a) : (b))
  30 
  31 void exo(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne0,
  32          double *v,double *stn,ITG *inum,ITG *nmethod,ITG *kode,
  33          char *filab,double *een,double *t1,double *fn,double *time,
  34          double *epn,ITG *ielmat,char *matname,double *enern,
  35          double *xstaten,ITG *nstate_,ITG *istep,ITG *iinc,
  36          ITG *ithermal,double *qfn,ITG *mode,ITG *noddiam,
  37          double *trab,ITG *inotr,ITG *ntrans,double *orab,
  38          ITG *ielorien,ITG *norien,char *description,ITG *ipneigh,
  39          ITG *neigh,ITG *mi,double *stx,double *vr,double *vi,
  40          double *stnr,double *stni,double *vmax,double *stnmax,
  41          ITG *ngraph,double *veold,double *ener,ITG *ne,double *cs,
  42          char *set,ITG *nset,ITG *istartset,ITG *iendset,ITG *ialset,
  43          double *eenmax,double *fnr,double *fni,double *emn,
  44          double *thicke,char *jobnamec,char *output,double *qfx,
  45          double *cdn,ITG *mortar,double *cdnr,double *cdni,ITG *nmat){
  46 
  47   /* stores the results in exo format
  48 
  49      iselect selects which nodes are to be stored:
  50      iselect=-1 means only those nodes for which inum negative
  51      ist, i.e. network nodes
  52      iselect=+1 means only those nodes for which inum positive
  53      ist, i.e. structural nodes
  54      iselect=0  means both of the above */
  55   
  56   /* Note for frd strcmp1(output,"asc")==0 defines ascii file, otherwise binary */
  57   
  58   char fneig[132]="", material[6]="     ",text[2]=" ";
  59   
  60   static ITG nkcoords,nout,noutmin,noutplus;
  61   ITG nterms;
  62 
  63   ITG i,j,k,l,m,n,o,indexe,nemax,nlayer,noutloc,iset,iselect,ncomp,nope,
  64     nodes,ifield[7],nfield[2],icomp[7],ifieldstate[*nstate_],icompstate[*nstate_],nelout;
  65   
  66   ITG ncompscalar=1,ifieldscalar[1]={1},icompscalar[1]={0},nfieldscalar[2]={1,0};
  67   ITG ncompvector=3,ifieldvector[3]={1,1,1},icompvector[3]={0,1,2},nfieldvector1[2]={3,0},nfieldvector0[2]={mi[1]+1,0};
  68   ITG ncomptensor=6,ifieldtensor[6]={1,1,1,1,1,1},icomptensor[6]={0,1,2,3,5,4},nfieldtensor[2]={6,0};
  69   ITG ncompscalph=2,ifieldscalph[2]={1,2},icompscalph[2]={0,0},nfieldscalph[2]={0,0};
  70   ITG ncompvectph=6,ifieldvectph[6]={1,1,1,2,2,2},icompvectph[6]={1,2,3,1,2,3},nfieldvectph[2]={mi[1]+1,mi[1]+1};
  71   ITG ncomptensph=12,ifieldtensph[12]={1,1,1,1,1,1,2,2,2,2,2,2},icomptensph[12]={0,1,2,3,5,4,0,1,2,3,5,4},nfieldtensph[2]={6,6};
  72 
  73   int errr=0, exoid=0;
  74   ITG num_dim, num_elem;
  75   ITG num_elem_blk; /* Node element blocks.  One eltype per block*/
  76   ITG num_ns, num_ss, num_es, num_fs; /* Node sets, side sets, element sets, face sets */
  77   int CPU_word_size = sizeof(float);
  78   int IO_word_size = sizeof(float);
  79 
  80   /* Filename */
  81   strcpy (fneig, jobnamec);
  82   strcat (fneig, ".exo");
  83 
  84   /* nkcoords is the number of nodes at the time when 
  85      the nodal coordinates are stored in the exo file */
  86   nkcoords = *nk;
  87   ITG num_nodes = nkcoords;
  88 
  89   /* determining nout, noutplus and noutmin 
  90      nout: number of structural and network nodes
  91      noutplus: number of structural nodes
  92      noutmin: number of network nodes */
  93   if(*nmethod!=0){
  94     nout=0;
  95     noutplus=0;
  96     noutmin=0;
  97     for(i=0;i<*nk;i++){
  98       if(inum[i]==0) continue;
  99       if(inum[i]>0) noutplus++;
 100       if(inum[i]<0) noutmin++;
 101       nout++;
 102     }
 103   }else{
 104     nout=*nk;
 105   }
 106 
 107   // Allocate memory for node positions
 108   float *x, *y, *z;
 109   x = (float *) calloc(nout, sizeof(float));
 110   y = (float *) calloc(nout, sizeof(float));
 111   z = (float *) calloc(nout, sizeof(float));
 112 
 113   // Write optional node map
 114   j = 0; // Counter for the exo order of the nodes
 115   ITG *node_map,*node_map_inv;
 116   node_map = (ITG *) calloc(nout, sizeof(ITG));
 117   node_map_inv = (ITG *) calloc(nkcoords, sizeof(ITG));
 118   /* storing the coordinates of the nodes */
 119   if(*nmethod!=0){
 120     for(i=0;i<*nk;i++){
 121       if(inum[i]==0){continue;}
 122       // The difference between i and j is that not all values of i
 123       // increment j.
 124       node_map[j] = i+1;
 125       node_map_inv[i] = j+1;
 126       x[j]   = co[3*i];
 127       y[j]   = co[3*i+1];
 128       z[j++] = co[3*i+2];
 129     }
 130   }else{
 131     for(i=0;i<*nk;i++){
 132       node_map[j] = i+1;
 133       node_map_inv[i] = j+1;
 134       x[j]   = co[3*i];
 135       y[j]   = co[3*i+1];
 136       z[j++] = co[3*i+2];
 137     }    
 138   }    
 139 
 140   /* first time something is written in the exo-file: store
 141      computational metadata, the nodal coordinates and the
 142      topology */
 143   if(*kode==1){
 144     num_dim = 3;  
 145     num_elem_blk = 19;
 146     
 147     // Find the number of sets
 148     exosetfind(set, nset, ialset, istartset, iendset, 
 149                &num_ns, &num_ss, &num_es, &num_fs, NULL, exoid, (int) 0, nk);
 150     printf ("Side sets to exo file not implemented.\n");
 151     num_ss=0;
 152 
 153 #ifdef LONGLONG
 154     // This handles LONGLONG transparently
 155     remove(fneig);
 156     exoid = ex_create (fneig, /*Filename*/
 157                        EX_ALL_INT64_API,        /* create mode */
 158                        &CPU_word_size,  /* CPU float word size in bytes */
 159                        &IO_word_size);  /* I/O float word size in bytes */
 160 #else
 161     exoid = ex_create (fneig, /*Filename*/
 162                        EX_CLOBBER,      /* create mode */
 163                        &CPU_word_size,  /* CPU float word size in bytes */
 164                        &IO_word_size);  /* I/O float word size in bytes */
 165 #endif
 166     
 167     /* determining the number of elements */
 168     if(*nmethod!=0){
 169       nelout=0;
 170       for(i=0;i<*ne0;i++){
 171         if(ipkon[i]<-1){
 172           continue;
 173         }else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
 174           continue;
 175         }else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
 176           continue;
 177         }else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
 178           continue;
 179         }else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
 180           // Count the number of layers
 181           nlayer=0;
 182           for(k=0;k<mi[2];k++){
 183             if(ielmat[i*mi[2]+k]==0) break;
 184             nlayer++;
 185           }
 186           // Allow an element for each layer
 187           for(k=0;k<nlayer;k++){
 188             nelout++;
 189           }
 190         }else{
 191           nelout++;
 192         }
 193       }
 194     }else{
 195       nelout=*ne;
 196     }
 197     num_elem = nelout;
 198     
 199     /* initialize file with parameters */
 200     printf("\nData writen to the .exo file\n");
 201     num_nodes=nout;
 202     printf("Number of nodes: %" ITGFORMAT "\n", num_nodes);
 203     printf("Number of elements %" ITGFORMAT "\n", num_elem);
 204     printf("Number of element blocks %" ITGFORMAT "\n", num_elem_blk);
 205     printf("Number of node sets %" ITGFORMAT "\n", num_ns);
 206     printf("Number of side sets %" ITGFORMAT "\n", num_ss);
 207 
 208     errr = ex_put_init (exoid, "CalculiX EXO File", 
 209                         num_dim, num_nodes, 
 210                         num_elem, num_elem_blk, 
 211                         num_ns, num_ss);
 212     
 213     // Write values to database
 214     errr = ex_put_coord (exoid, x, y, z);
 215     if(errr)printf("*ERROR in exo: failed node positions");
 216     errr = ex_put_node_num_map (exoid, node_map);
 217     if(errr)printf("*ERROR in exo: failed node map");
 218     
 219     // Deallocate
 220     free (node_map);
 221     free (x);
 222     free (y);
 223     free (z);
 224 
 225     // Write coordinate names
 226     char *coord_names[3];
 227     coord_names[0] = "x";
 228     coord_names[1] = "y";
 229     coord_names[2] = "z";
 230     errr = ex_put_coord_names (exoid, coord_names);
 231     if(errr){printf("*ERROR in exo: failed coordinate names");}
 232     
 233     // Initialize enough memory to store the element numbers
 234     ITG *elem_map;
 235     elem_map = (ITG *) calloc(num_elem, sizeof(ITG));
 236     char curblk[6]="     ";
 237     
 238     ITG *blkassign;
 239     blkassign = (ITG *) calloc(num_elem, sizeof(ITG));
 240     
 241     l=0;
 242     for(i=0;i<*ne0;i++){ // For each element.  Composite elements are
 243                          // one increment in this loop and all layers
 244                          // have the same element number.
 245       if(ipkon[i]<=-1){
 246         continue;
 247       }else if(strcmp1(&lakon[8*i],"F")==0){
 248         continue;
 249       }else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
 250         continue;
 251       }else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
 252         continue;
 253       }else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
 254         continue;
 255       }else{
 256         indexe=ipkon[i];
 257       }
 258       
 259       elem_map[l] = i+1;
 260       
 261       strcpy1(curblk,&lakon[8*i],5);
 262       strcpy1(material,&matname[80*(ielmat[i*mi[2]]-1)],5);
 263       if(strcmp1(&lakon[8*i+3],"2")==0){
 264         /* 20-node brick element */
 265         if(((strcmp1(&lakon[8*i+6]," ")==0)||
 266             (strcmp1(&filab[4],"E")==0)||
 267             (strcmp1(&lakon[8*i+6],"I")==0))&&
 268            (strcmp2(&lakon[8*i+6],"LC",2)!=0)){
 269           blkassign[l++]=1;
 270         }else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
 271           /* composite material */
 272           nlayer=0;
 273           for(k=0;k<mi[2];k++){
 274             if(ielmat[i*mi[2]+k]==0) break;
 275             nlayer++;
 276           }
 277           for(k=0;k<nlayer;k++){
 278             nemax++;
 279             elem_map[l] = i+1;
 280             blkassign[l++]=2;
 281           }
 282         }else if(strcmp1(&lakon[8*i+6],"B")==0){
 283           /* 3-node beam element */
 284           blkassign[l++]=3;
 285         }else{
 286           /* 8-node 2d element */
 287           blkassign[l++]=4;
 288         }
 289       }else if(strcmp1(&lakon[8*i+3],"8")==0){
 290         if((strcmp1(&lakon[8*i+6]," ")==0)||
 291            (strcmp1(&filab[4],"E")==0)){
 292           /* 8-node brick element */
 293           blkassign[l++]=5;
 294         }else if(strcmp1(&lakon[8*i+6],"B")==0){
 295           /* 2-node 1d element */
 296           if(strcmp1(&lakon[8*i+4],"R")==0){
 297             blkassign[l++]=6;
 298           }else if(strcmp1(&lakon[8*i+4],"I")==0){
 299             blkassign[l++]=7;
 300           }
 301         }else{
 302           /* 4-node 2d element */
 303           if(strcmp1(&lakon[8*i+4],"R")==0){
 304             blkassign[l++]=8;
 305           }else if(strcmp1(&lakon[8*i+4],"I")==0){
 306             blkassign[l++]=9;
 307           }
 308         }
 309       }else if((strcmp1(&lakon[8*i+3],"10")==0)||
 310                (strcmp1(&lakon[8*i+3],"14")==0)){
 311         /* 10-node tetrahedral element */
 312         blkassign[l++]=10;
 313       }else if(strcmp1(&lakon[8*i+3],"4")==0){
 314         /* 4-node tetrahedral element */
 315         blkassign[l++]=11;
 316       }else if(strcmp1(&lakon[8*i+3],"15")==0){
 317         if((strcmp1(&lakon[8*i+6]," ")==0)||
 318            (strcmp1(&filab[4],"E")==0)){
 319           /* 15-node wedge element */
 320           blkassign[l++]=12;
 321         }else{
 322           /* 6-node 2d element */
 323           blkassign[l++]=13;
 324         }
 325       }else if(strcmp1(&lakon[8*i+3],"6")==0){
 326         if((strcmp1(&lakon[8*i+6]," ")==0)||
 327            (strcmp1(&filab[4],"E")==0)){
 328           /* 6-node wedge element */
 329           blkassign[l++]=14;
 330         }else{
 331           /* 3-node 2d element */
 332           blkassign[l++]=15;
 333         }
 334         //      }else if((strcmp1(&lakon[8*i],"D")==0)&&
 335         //             (strcmp1(&lakon[8*i],"DCOUP3D")!=0)){
 336       }else if(strcmp1(&lakon[8*i],"D")==0){
 337         if(kon[indexe]==0){
 338           /* 2-node 1d element (network entry element) */
 339           blkassign[l++]=16;
 340         }else if(kon[indexe+2]==0){
 341           /* 2-node 1d element (network exit element) */
 342           blkassign[l++]=17;
 343         }else{
 344           /* 3-node 1d element (genuine network element) */
 345           blkassign[l++]=18;
 346         }
 347       }else if((strcmp1(&lakon[8*i],"E")==0)&&
 348                (strcmp1(&lakon[8*i+6],"A")==0)){
 349         /* 2-node 1d element (spring element) */
 350         blkassign[l++]=19;
 351       }
 352     }
 353 
 354 
 355     
 356     int num_nodes_per_elem[num_elem_blk];
 357     char *blknames[num_elem_blk];
 358     j=0;
 359     num_nodes_per_elem[j]=1;   blknames[j++]="PNT";
 360     num_nodes_per_elem[j]=20;  blknames[j++]="C3D20 or C3D20R or S8R";
 361     num_nodes_per_elem[j]=20;  blknames[j++]="COMPOSITE LAYER C3D20";
 362     num_nodes_per_elem[j]=3;   blknames[j++]="Beam B32 or B32R";
 363     num_nodes_per_elem[j]=8;   blknames[j++]="CPS8 or CPE8 or CAX8";
 364     num_nodes_per_elem[j]=8;   blknames[j++]="C3D8 or C3D8R";
 365     num_nodes_per_elem[j]=2;   blknames[j++]="TRUSS2";
 366     num_nodes_per_elem[j]=2;   blknames[j++]="TRUSS2";
 367     num_nodes_per_elem[j]=4;   blknames[j++]="CPS4R or CPE4R";
 368     num_nodes_per_elem[j]=4;   blknames[j++]="CPS4I or CPE4I";
 369     num_nodes_per_elem[j]=10;  blknames[j++]="C3D10";
 370     num_nodes_per_elem[j]=4;   blknames[j++]="C3D4";
 371     num_nodes_per_elem[j]=15;  blknames[j++]="C3D15";
 372     num_nodes_per_elem[j]=6;   blknames[j++]="CPS6 or CPE6";
 373     num_nodes_per_elem[j]=6;   blknames[j++]="C3D6";
 374     num_nodes_per_elem[j]=3;   blknames[j++]="CPS3 or CPE3";
 375     num_nodes_per_elem[j]=2;   blknames[j++]="2-node 1d network entry elem";
 376     num_nodes_per_elem[j]=2;   blknames[j++]="2-node 1d network exit elem";
 377     num_nodes_per_elem[j]=3;   blknames[j++]="2-node 1d genuine network elem";
 378     num_nodes_per_elem[j]=2;   blknames[j++]="2-node 1d spring elem";
 379 
 380     errr = ex_put_names (exoid, EX_ELEM_BLOCK, blknames);
 381     if(errr){printf("*ERROR in exo: cannot write block names");}
 382     
 383 
 384     /* write element connectivity */
 385     ITG *connect;
 386     ITG num_elem_in_blk;
 387     ITG blksize[num_elem_blk];
 388 
 389     for(l=0;l<num_elem_blk;l++){
 390       // First determine the size of the block
 391       j=0;
 392       m=0;
 393       for(i=0;i<*ne0;i++){
 394         if(ipkon[i]<0) continue;
 395         if(blkassign[m]==l){
 396           if(blkassign[m]==2){//composite
 397             for(k=0;k<mi[2];k++){
 398               if(ielmat[i*mi[2]+k]==0) break;
 399               j++;
 400             }
 401           }else{
 402             j++;
 403           }
 404         }
 405         m++;
 406       }
 407 
 408       blksize[l]=j;
 409       num_elem_in_blk=blksize[l];
 410       
 411       connect = (ITG *) calloc (num_elem_in_blk*num_nodes_per_elem[l], sizeof(ITG));
 412       // printf ("Size of connect %" ITGFORMAT "\n", num_elem_in_blk*num_nodes_per_elem[l]*sizeof(ITG));
 413       k=0; o=0;
 414       // Now connectivity
 415       for(i=0;i<*ne0;i++){
 416         if(ipkon[i]<0) continue;
 417         indexe=ipkon[i];
 418         if (blkassign[o]==l){
 419           // printf ("block assignment %" ITGFORMAT "\n", blkassign[o]);
 420           if(blkassign[o]==1){ // S8R and all? C3D20 promotions?
 421             for(m=0;m<12;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 422             for(m=16;m<20;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 423             for(m=12;m<16;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 424           }else if(blkassign[o]==12){ // C3D15
 425             for(m=0;m<9;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 426             for(m=12;m<15;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 427             for(m=9;m<12;m++){connect[k++] = node_map_inv[kon[indexe+m]-1];}
 428           }else if (blkassign[o]==2){ // Composite
 429             nlayer=0;
 430             for(l=0;l<mi[2];l++){
 431               if(ielmat[i*mi[2]+l]==0) break;
 432               nlayer++;
 433             }
 434             for(n=0;n<nlayer;n++){
 435               for(m=0;m<12;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
 436               for(m=16;m<20;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
 437               for(m=12;m<16;m++){connect[k++] = node_map_inv[kon[indexe+28+20*n+m]-1];}
 438             }
 439           }else if(blkassign[o]==4){ // 8 Node 2D elements CAX8
 440             for (j = 0; j <num_nodes_per_elem[l]; j++){
 441               connect[k++] = node_map_inv[kon[indexe+20+j]-1];
 442             }
 443           }else {
 444             for (j = 0; j <num_nodes_per_elem[l]; j++){
 445               connect[k++] = node_map_inv[kon[indexe+j]-1];
 446             }
 447           }
 448         }
 449         o++;
 450       }
 451       
 452       int num_attr=0;
 453       switch (l)
 454         {
 455         case 0:
 456           errr = ex_put_elem_block (exoid, l, "SPHERE", num_elem_in_blk, num_nodes_per_elem[l], num_attr);
 457           break;
 458         case 1:
 459           errr = ex_put_elem_block (exoid, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_attr);   
 460           break;
 461         case 2:
 462           errr = ex_put_elem_block (exoid, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_attr);   
 463           break;
 464         case 3:
 465           errr = ex_put_elem_block (exoid, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 466           break;
 467         case 4:
 468           errr = ex_put_elem_block (exoid, l, "QUAD", num_elem_in_blk, num_nodes_per_elem[l], num_attr);          
 469           break;
 470         case 5:
 471           errr = ex_put_elem_block (exoid, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_attr);   
 472           break;
 473         case 6:
 474           errr = ex_put_elem_block (exoid, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 475           break;
 476         case 7:
 477           errr = ex_put_elem_block (exoid, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 478           break;
 479         case 8:
 480           errr = ex_put_elem_block (exoid, l, "SHELL", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 481           break;
 482         case 9:
 483           errr = ex_put_elem_block (exoid, l, "SHELL", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 484           break;
 485         case 10:
 486           errr = ex_put_elem_block (exoid, l, "TETRA", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 487           break;
 488         case 11:
 489           errr = ex_put_elem_block (exoid, l, "TETRA", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 490           break;
 491         case 12:
 492           errr = ex_put_elem_block (exoid, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 493           break;
 494         case 13:
 495           errr = ex_put_elem_block (exoid, l, "HEX", num_elem_in_blk, num_nodes_per_elem[l], num_attr);   
 496           break;
 497         case 14:
 498           errr = ex_put_elem_block (exoid, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 499           break;
 500         case 15:
 501           errr = ex_put_elem_block (exoid, l, "WEDGE", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 502           break;
 503         default:
 504           // case 16:
 505           // case 17:
 506           // case 18:
 507           // case 19:
 508           errr = ex_put_elem_block (exoid, l, "TRUSS", num_elem_in_blk, num_nodes_per_elem[l], num_attr);         
 509           break;
 510         };
 511           
 512 
 513 
 514       if (num_elem_in_blk>0){
 515         errr = ex_put_elem_conn (exoid, l, connect);
 516         if (errr)
 517           printf ("ERROR in ex_put_elem_conn %i\n", errr);
 518       }
 519       free (connect);
 520     }
 521     
 522     // Write the element map into the file
 523     errr = ex_put_elem_num_map (exoid, elem_map); 
 524     if (errr)
 525       printf ("ERROR in ex_put_elem_num_map %i\n", errr);
 526     
 527     // Write the node sets into the file
 528     exosetfind(set, nset, ialset, istartset, iendset, 
 529                &num_ns, &num_ss, &num_es, &num_fs, node_map_inv, exoid, (int) 1, nk);
 530     
 531     // Free up memory which is gathering dust
 532     free (elem_map); 
 533     free (blkassign);
 534 
 535     // Close files
 536     ex_update (exoid);  
 537     ex_close (exoid);
 538 
 539     if(*nmethod==0){return;}
 540     
 541     /* End of if(*kode==1) */
 542   }
 543 
 544 
 545 
 546 
 547   /* Start the data storage section */
 548   float version;
 549   exoid = ex_open (fneig, /*Filename*/
 550                    EX_WRITE,    /* create mode */
 551                    &CPU_word_size,  /* CPU float word size in bytes */
 552                    &IO_word_size, /* I/O float word size in bytes */
 553                    &version);
 554     
 555   // Define an output variable for general use  
 556   float *nodal_var_vals;
 557 
 558   int num_time_steps;
 559   float fdum;
 560   char cdum = 0;
 561   
 562   errr = ex_inquire (exoid, EX_INQ_TIME, &num_time_steps, &fdum, &cdum);
 563   printf ("Time periods existing in exo file: %i\n", num_time_steps);
 564   ++num_time_steps;
 565   float timet = (float) *time;
 566   // printf ("Writing time period %" ITGFORMAT " at time %f\n", num_time_steps, *time);
 567   errr = ex_put_time (exoid, num_time_steps, &timet);
 568   if (errr) printf ("Error storing time into exo file.\n");
 569 
 570   // Statically allocate an array of pointers.  Can't figure out how to do this dynamically.
 571   // 100 should be enough to store the variable names.
 572   char *var_names[100];
 573 
 574   int countvars=0;
 575   int countbool=3;
 576 
 577   while(countbool>0){
 578     // First time count
 579     // Second time assign names
 580     // Last time save data
 581     
 582     /*  for cyclic symmetry frequency calculations only results for
 583         even numbers (= odd modes, numbering starts at 0) are stored */
 584     if((*nmethod==2)&&(((*mode/2)*2!=*mode)&&(*noddiam>=0))){ex_close(exoid);return;}
 585     
 586     /* storing the displacements in the nodes */
 587     if(strcmp1(filab,"U ")==0){
 588       if (countbool==3){
 589         countvars+=3;
 590       }else if(countbool==2){
 591         var_names[countvars++]="Ux";
 592         var_names[countvars++]="Uy";
 593         var_names[countvars++]="Uz";
 594       }else{
 595         iselect=1;
 596         
 597         frdset(filab,set,&iset,istartset,iendset,ialset,
 598                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 599                ngraph);
 600         
 601         exovector(v,&iset,ntrans,filab,&nkcoords,inum,inotr,
 602                   trab,co,istartset,iendset,ialset,mi,ngraph,exoid,
 603                   num_time_steps,countvars,nout);
 604         countvars+=3;
 605       }
 606     }
 607 
 608     /*     storing the imaginary part of displacements in the nodes
 609            for the odd modes of cyclic symmetry calculations */
 610     if(*noddiam>=0){
 611       if(strcmp1(filab,"U ")==0){
 612         if (countbool==3){
 613           countvars+=3;
 614         }else if(countbool==2){
 615           var_names[countvars++]="U-imag-x";
 616           var_names[countvars++]="U-imag-y";
 617           var_names[countvars++]="U-imag-z";
 618         }else{
 619           exovector(&v[*nk*(mi[1]+1)],&iset,ntrans,filab,&nkcoords,inum,inotr,
 620                     trab,co,istartset,iendset,ialset,mi,ngraph,exoid,
 621                     num_time_steps,countvars,nout);
 622 
 623           countvars+=3;
 624         }
 625       }
 626     }
 627     
 628     /* storing the velocities in the nodes */
 629     if(strcmp1(&filab[1740],"V   ")==0){
 630       if (countbool==3){
 631         countvars+=3;
 632       }else if(countbool==2){
 633         var_names[countvars++]="Vx";
 634         var_names[countvars++]="Vy";
 635         var_names[countvars++]="Vz";
 636       }else{
 637         iselect=1;
 638 
 639         frdset(&filab[1740],set,&iset,istartset,iendset,ialset,
 640                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 641                ngraph);
 642     
 643         exovector(veold,&iset,ntrans,&filab[1740],&nkcoords,inum,inotr,
 644                   trab,co,istartset,iendset,ialset,mi,ngraph,exoid,
 645                   num_time_steps,countvars,nout);
 646 
 647         countvars+=3;
 648       }
 649     }
 650 
 651     /* storing the temperatures in the nodes */
 652     if(strcmp1(&filab[87],"NT  ")==0){
 653       if (countbool==3){
 654         countvars+=1;
 655       }else if(countbool==2){
 656         var_names[countvars++]="NT";
 657       }else{
 658         iselect=0;
 659         
 660         frdset(&filab[87],set,&iset,istartset,iendset,ialset,
 661                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 662                ngraph);
 663         
 664         if(*ithermal<=1){
 665           exoselect(t1,t1,&iset,&nkcoords,inum,istartset,iendset,
 666                     ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
 667                     nfieldscalar,&iselect,exoid,num_time_steps,countvars,nout);
 668         }else{
 669           exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
 670                     ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
 671                     nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
 672         }
 673         countvars+=1;   
 674       }
 675     }
 676     
 677     /* storing the stresses in the nodes */
 678     if(strcmp1(&filab[174],"S   ")==0){
 679       if (countbool==3){
 680         countvars+=6;
 681       }else if(countbool==2){
 682         // Note reordered (done in exoselect) relative to frd file
 683         // Order must be xx  yy  zz  xy  xz  yz
 684         var_names[countvars++]="Sxx";
 685         var_names[countvars++]="Syy";
 686         var_names[countvars++]="Szz";
 687         var_names[countvars++]="Sxy";
 688         var_names[countvars++]="Sxz"; 
 689         var_names[countvars++]="Syz";
 690       }else{
 691         iselect=1;
 692     
 693         frdset(&filab[174],set,&iset,istartset,iendset,ialset,
 694                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 695                ngraph);
 696     
 697         exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
 698                   ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 699                   nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 700         countvars+=6;
 701       }
 702     }
 703 
 704     /* storing the imaginary part of the stresses in the nodes
 705        for the odd modes of cyclic symmetry calculations */
 706     if(*noddiam>=0){
 707       if (countbool==3){
 708         countvars+=6;
 709       }else if(countbool==2){
 710         // Note reordered (done in exoselect) relative to frd file
 711         // Order must be xx  yy  zz  xy  xz  yz
 712         var_names[countvars++]="S-imagxx";
 713         var_names[countvars++]="S-imagyy";
 714         var_names[countvars++]="S-imagzz";
 715         var_names[countvars++]="S-imagxy";
 716         var_names[countvars++]="S-imagzx";
 717         var_names[countvars++]="S-imagyz";
 718       }else{
 719         if(strcmp1(&filab[174],"S   ")==0){      
 720           exoselect(&stn[6**nk],stn,&iset,&nkcoords,inum,istartset,iendset,
 721                     ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 722                     nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 723           countvars+=6;
 724         }
 725       }
 726     }
 727 
 728     /* storing the total strains in the nodes */
 729     if(strcmp1(&filab[261],"E   ")==0){
 730       if (countbool==3){
 731         countvars+=6;
 732       }else if(countbool==2){
 733         // Note reordered (done in exoselect) relative to frd file
 734         // Order must be xx  yy  zz  xy  xz  yz
 735         var_names[countvars++]="Exx";
 736         var_names[countvars++]="Eyy";
 737         var_names[countvars++]="Ezz";
 738         var_names[countvars++]="Exy";
 739         var_names[countvars++]="Exz";
 740         var_names[countvars++]="Eyz";
 741       }else{
 742         iselect=1;
 743         
 744     
 745         frdset(&filab[261],set,&iset,istartset,iendset,ialset,
 746                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 747                ngraph);
 748         
 749         exoselect(een,een,&iset,&nkcoords,inum,istartset,iendset,
 750                   ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 751                   nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 752 
 753         countvars+=6;
 754       }
 755     }
 756     
 757     /* storing the imaginary part of the total strains in the nodes
 758        for the odd modes of cyclic symmetry calculations */
 759     if(*noddiam>=0){
 760       if(strcmp1(&filab[261],"E   ")==0){
 761         if (countbool==3){
 762           countvars+=6;
 763         }else if(countbool==2){
 764           // Note reordered (done in exoselect) relative to frd file
 765           // Order must be xx  yy  zz  xy  xz  yz
 766           var_names[countvars++]="E-imagxx";
 767           var_names[countvars++]="E-imagyy";
 768           var_names[countvars++]="E-imagzz";
 769           var_names[countvars++]="E-imagxy";
 770           var_names[countvars++]="E-imagxz";
 771           var_names[countvars++]="E-imagyz";
 772         }else{
 773           exoselect(&een[6**nk],een,&iset,&nkcoords,inum,istartset,iendset,
 774                     ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 775                     nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 776 
 777           countvars+=6;
 778         }
 779       }
 780     }
 781     
 782     /* storing the mechanical strains in the nodes */
 783     if(strcmp1(&filab[2697],"ME  ")==0){
 784       if (countbool==3){
 785         countvars+=6;
 786       }else if(countbool==2){
 787         // Note reordered (done in exoselect) relative to frd file
 788         // Order must be xx  yy  zz  xy  xz  yz
 789         var_names[countvars++]="MExx";
 790         var_names[countvars++]="MEyy";
 791         var_names[countvars++]="MEzz";
 792         var_names[countvars++]="MExy";
 793         var_names[countvars++]="MExz";
 794         var_names[countvars++]="MEyz";
 795       }else{
 796         iselect=1;
 797         frdset(&filab[2697],set,&iset,istartset,iendset,ialset,
 798                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 799                ngraph);
 800         exoselect(emn,emn,&iset,&nkcoords,inum,istartset,iendset,
 801                   ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 802                   nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 803         countvars+=6;
 804       }    
 805     }
 806 
 807     /* storing the imaginary part of the mechanical strains in the nodes
 808        for the odd modes of cyclic symmetry calculations */
 809     if(*noddiam>=0){
 810       if(strcmp1(&filab[2697],"ME  ")==0){
 811         if (countbool==3){
 812           countvars+=6;
 813         }else if(countbool==2){
 814           // Note reordered (done in exoselect) relative to frd file
 815           // Order must be xx  yy  zz  xy  xz  yz
 816           var_names[countvars++]="ME-imagxx";
 817           var_names[countvars++]="ME-imagyy";
 818           var_names[countvars++]="ME-imagzz";
 819           var_names[countvars++]="ME-imagxy";
 820           var_names[countvars++]="ME-imagxz";
 821           var_names[countvars++]="ME-imagyz";
 822         }else{
 823           exoselect(&emn[6**nk],een,&iset,&nkcoords,inum,istartset,iendset,
 824                     ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
 825                     nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
 826           countvars+=6;      
 827         }
 828       }
 829     }
 830 
 831     /* storing the forces in the nodes */
 832     
 833     if(strcmp1(&filab[348],"RF  ")==0){
 834       if (countbool==3){
 835         countvars+=3;
 836       }else if(countbool==2){
 837         var_names[countvars++]="RFx";
 838         var_names[countvars++]="RFy";
 839         var_names[countvars++]="RFz";
 840       }else{
 841         iselect=1;
 842         frdset(&filab[348],set,&iset,istartset,iendset,ialset,
 843                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 844                ngraph);
 845         exovector(fn,&iset,ntrans,&filab[348],&nkcoords,inum,inotr,
 846                   trab,co,istartset,iendset,ialset,mi,ngraph,exoid,
 847                   num_time_steps,countvars,nout);
 848         countvars+=3;
 849       }
 850     }
 851     
 852     /*     storing the imaginary part of the forces in the nodes
 853            for the odd modes of cyclic symmetry calculations */
 854     if(*noddiam>=0){
 855       if(strcmp1(&filab[348],"RF  ")==0){
 856         if (countbool==3){
 857           countvars+=3;
 858         }else if(countbool==2){
 859           var_names[countvars++]="RF-imagx";
 860           var_names[countvars++]="RF-imagy";
 861           var_names[countvars++]="RF-imagz";
 862         }else{
 863           exovector(&fn[*nk*(mi[1]+1)],&iset,ntrans,filab,&nkcoords,inum,inotr,
 864                     trab,co,istartset,iendset,ialset,mi,ngraph,exoid,
 865                     num_time_steps,countvars,nout);
 866           countvars+=3;
 867         }
 868       }
 869     }
 870 
 871     /* storing the equivalent plastic strains in the nodes */
 872     if(strcmp1(&filab[435],"PEEQ")==0){
 873       if (countbool==3){
 874         countvars+=1;
 875       }else if(countbool==2){
 876         var_names[countvars++]="PEEQ";
 877       }else{
 878         iselect=1;
 879         frdset(&filab[435],set,&iset,istartset,iendset,ialset,
 880                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 881                ngraph);
 882         exoselect(epn,epn,&iset,&nkcoords,inum,istartset,iendset,
 883                   ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
 884                   nfieldscalar,&iselect,exoid,num_time_steps,countvars,nout);
 885         printf ("Warning: export PEEQ to exo not tested.\n");
 886         countvars+=1;
 887       }
 888     }
 889       
 890 
 891     /* storing the energy in the nodes */
 892     if(strcmp1(&filab[522],"ENER")==0){
 893       if (countbool==3){
 894         countvars+=1;
 895       }else if(countbool==2){
 896         var_names[countvars++]="ENER";
 897       }else{
 898         iselect=1;
 899         frdset(&filab[522],set,&iset,istartset,iendset,ialset,
 900                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 901                ngraph);
 902         exoselect(enern,enern,&iset,&nkcoords,inum,istartset,iendset,
 903                   ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
 904                   nfieldscalar,&iselect,exoid,num_time_steps,countvars,nout);
 905         printf ("Warning: export ENER to exo not tested.\n");
 906         countvars+=1;
 907       }
 908     }
 909   
 910     /* storing the contact displacements and stresses at the slave nodes */
 911     if(strcmp1(&filab[2175],"CONT")==0){
 912       if (countbool==3){
 913         countvars+=6;
 914       }else if(countbool==2){
 915         var_names[countvars++]="COPEN";
 916         var_names[countvars++]="CSLIP1";
 917         var_names[countvars++]="CSLIP2";
 918         var_names[countvars++]="CPRESS";
 919         var_names[countvars++]="CSHEAR1";
 920         var_names[countvars++]="CSHEAR2";
 921       }else{
 922         
 923         nodal_var_vals = (float *) calloc (nout, sizeof(float));
 924         
 925         for(i=*ne-1;i>=0;i--){
 926           if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
 927             break;
 928         }
 929         noutloc=*ne-i-1;
 930         
 931         for(j=0;j<6;j++){
 932           for(i=*ne-1;i>=0;i--){
 933             if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
 934               break;
 935             strcpy1(text,&lakon[8*i+7],1);
 936             nope=atoi(text)+1;
 937             nodes=node_map_inv[kon[ipkon[i]+nope-1]-1]-1;
 938             nodal_var_vals[nodes]=stx[6*mi[0]*i+j]; 
 939           }
 940           
 941           errr = ex_put_nodal_var (exoid, num_time_steps, 1+countvars++, nout, nodal_var_vals);
 942           if (errr) printf ("ERROR storing contact data into exo file.\n");
 943         }
 944         
 945         free(nodal_var_vals);
 946       }
 947     }
 948     
 949     /* storing the contact energy at the slave nodes */
 950     if(strcmp1(&filab[2262],"CELS")==0){
 951       if (countbool==3){
 952         countvars+=1;
 953       }else if(countbool==2){
 954         var_names[countvars++]="CELS";
 955       }else{
 956         for(i=*ne-1;i>=0;i--){
 957           if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
 958             break;
 959         }
 960         noutloc=*ne-i-1;
 961             
 962         nodal_var_vals = (float *) calloc (nkcoords, sizeof(float));
 963         for(i=*ne-1;i>=0;i--){
 964           if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
 965             break;
 966           nope=atoi(&lakon[8*i+7])+1;
 967           nodes=node_map_inv[kon[ipkon[i]+nope-1]-1]-1;
 968           nodal_var_vals[nodes]=ener[i*mi[0]];
 969         }
 970         
 971         countvars+=1;
 972         int errr = ex_put_nodal_var (exoid, num_time_steps, countvars, nout, nodal_var_vals);
 973         if (errr) printf ("ERROR storing CELS data into exo file.\n");
 974 
 975         free(nodal_var_vals);
 976       }
 977     }
 978   
 979     /* storing the internal state variables in the nodes */
 980     if(strcmp1(&filab[609],"SDV ")==0){
 981       if (countbool==3){
 982         countvars+=*nstate_;
 983       }else if(countbool==2){
 984         for(j=1;j<=*nstate_;j++){
 985           sprintf(var_names[countvars++],"SDV%" ITGFORMAT,j);
 986         }
 987       }else{
 988         iselect=1;
 989         
 990         frdset(&filab[609],set,&iset,istartset,iendset,ialset,
 991                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
 992                ngraph);
 993         
 994         for(i=0;i<*nstate_;i++){
 995           ifieldstate[i]=1;icompstate[i]=i;
 996         }
 997         nfield[0]=*nstate_;
 998         
 999         exoselect(xstaten,xstaten,&iset,&nkcoords,inum,istartset,iendset,
1000                   ialset,ngraph,nstate_,ifieldstate,icompstate,
1001                   nfield,&iselect,exoid,num_time_steps,countvars,nout);
1002         countvars+=*nstate_;
1003       }
1004     }
1005   
1006     /* storing the heat flux in the nodes
1007        the heat flux has been extrapolated from the integration points
1008        in subroutine extropolate.f, taking into account whether the 
1009        results are requested in the global system or in a local system.
1010        Therefore, subroutine frdvector cannot be used, since it assumes
1011        the values are stored in the global system */
1012     if((strcmp1(&filab[696],"HFL ")==0)&&(*ithermal>1)){
1013       if (countbool==3){
1014         countvars+=3;
1015       }else if(countbool==2){
1016         var_names[countvars++]="HFLx";
1017         var_names[countvars++]="HFLy";
1018         var_names[countvars++]="HFLz";
1019       }else{
1020         iselect=1;
1021       
1022         frdset(&filab[696],set,&iset,istartset,iendset,ialset,
1023                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1024                ngraph);
1025     
1026         exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
1027                   ialset,ngraph,&ncompvector,ifieldvector,icompvector,
1028                   nfieldvector1,&iselect,exoid,num_time_steps,countvars,nout);
1029         countvars+=3;
1030       }
1031     }
1032 
1033     /* storing the heat generation in the nodes */
1034     if((strcmp1(&filab[783],"RFL ")==0)&&(*ithermal>1)){
1035       if (countbool==3){
1036         countvars+=1;
1037       }else if(countbool==2){
1038         var_names[countvars++]="RFL";
1039       }else{
1040         iselect=1;
1041       
1042         frdset(&filab[783],set,&iset,istartset,iendset,ialset,
1043                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1044                ngraph);
1045         exoselect(fn,fn,&iset,&nkcoords,inum,istartset,iendset,
1046                   ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1047                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1048         countvars+=1;
1049       }
1050     }
1051   
1052     /* storing the Zienkiewicz-Zhu improved stresses in the nodes */
1053     if(strcmp1(&filab[1044],"ZZS")==0){
1054       if (countbool==3){
1055         countvars+=6;
1056       }else if(countbool==2){
1057         // Note reordered (done in exoselect) relative to frd file
1058         // Order must be xx  yy  zz  xy  xz  yz
1059         var_names[countvars++]="ZZSxx";
1060         var_names[countvars++]="ZZSyy";
1061         var_names[countvars++]="ZZSzz";
1062         var_names[countvars++]="ZZSxy";
1063         var_names[countvars++]="ZZSxz";
1064         var_names[countvars++]="ZZSyz";
1065       }else{
1066         FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
1067                          stx,&mi[0]));
1068 
1069         iselect=1;
1070     
1071         frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
1072                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1073                ngraph);
1074     
1075         exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
1076                   ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1077                   nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
1078         printf ("Warning: export of ZZSTR to exo not tested.\n");
1079         countvars+=6;
1080       }
1081     }
1082 
1083     /* storing the imaginary part of the Zienkiewicz-Zhu 
1084        improved stresses in the nodes
1085        for the odd modes of cyclic symmetry calculations */
1086     
1087     if(*noddiam>=0){
1088       if(strcmp1(&filab[1044],"ZZS")==0){
1089         if (countbool==3){
1090           countvars+=6;
1091         }else if(countbool==2){
1092           // Note reordered (done in exoselect) relative to frd file
1093           // Order must be xx  yy  zz  xy  xz  yz
1094           var_names[countvars++]="ZZS-imagxx";
1095           var_names[countvars++]="ZZS-imagyy";
1096           var_names[countvars++]="ZZS-imagzz";
1097           var_names[countvars++]="ZZS-imagxy";
1098           var_names[countvars++]="ZZS-imagxz";
1099           var_names[countvars++]="ZZS-imagyz";
1100         }else{  
1101 
1102           FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
1103                            &stx[6*mi[0]**ne],&mi[0]));
1104           
1105           exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
1106                     ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1107                     nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
1108           printf ("Warning: export of ZZSTR imaginary to exo not tested.\n");
1109           countvars+=6;
1110         }
1111       }
1112     }
1113 
1114     /* storing the error estimator in the nodes */
1115   
1116     if(strcmp1(&filab[1044],"ERR")==0){
1117       if (countbool==3){
1118         countvars+=2;
1119       }else if(countbool==2){
1120         var_names[countvars++]="PSTDERROR";
1121         var_names[countvars++]="VMSTDERROR";
1122       }else{  
1123         
1124         nterms=6;
1125     
1126         FORTRAN(errorestimator,(stx,stn,ipkon,inum,kon,lakon,nk,ne,
1127                                 mi,ielmat,thicke,&nterms));
1128         
1129         iselect=1;
1130         
1131         frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
1132                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1133                ngraph);
1134         
1135         ncomp=2;
1136         ifield[0]=1;ifield[1]=1;
1137         icomp[0]=2;icomp[1]=4;
1138         
1139         exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
1140                   ialset,ngraph,&ncomp,ifield,icomp,
1141                   nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
1142 
1143         countvars+=2;
1144       }
1145     }
1146 
1147 
1148     /* storing the imaginary part of the error estimator in the nodes
1149        for the odd modes of cyclic symmetry calculations */
1150     if(*noddiam>=0){
1151       if(strcmp1(&filab[1044],"ERR")==0){
1152         if (countbool==3){
1153           countvars+=2;
1154         }else if(countbool==2){
1155           var_names[countvars++]="PSTD-imag-ERROR";
1156           var_names[countvars++]="VMSTD-imag-ERROR";
1157         }else{
1158           nterms=6;
1159           FORTRAN(errorestimator,(&stx[6*mi[0]**ne],stn,ipkon,inum,kon,lakon,nk,ne,
1160                                   mi,ielmat,thicke,&nterms));
1161           
1162           ncomp=2;
1163           ifield[0]=1;ifield[1]=1;
1164           icomp[0]=2;icomp[1]=4;
1165         
1166           exoselect(stn,stn,&iset,&nkcoords,inum,istartset,iendset,
1167                     ialset,ngraph,&ncomp,ifield,icomp,
1168                     nfieldtensor,&iselect,exoid,num_time_steps,countvars,nout);
1169 
1170           countvars+=2;
1171         }
1172       }
1173     }
1174 
1175     /* storing the thermal error estimator in the nodes */
1176     if(strcmp1(&filab[2784],"HER")==0){
1177       if (countbool==3){
1178         countvars+=1;
1179       }else if(countbool==2){
1180         var_names[countvars++]="HFLSTD";
1181       }else{
1182         nterms=3;
1183         FORTRAN(errorestimator,(qfx,qfn,ipkon,inum,kon,lakon,nk,ne,
1184                                 mi,ielmat,thicke,&nterms));
1185         
1186         iselect=1;
1187         frdset(&filab[2784],set,&iset,istartset,iendset,ialset,
1188                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1189                ngraph);
1190         
1191         ncomp=1;
1192         ifield[0]=1;
1193         icomp[0]=2;
1194         
1195         exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
1196                   ialset,ngraph,&ncomp,ifield,icomp,
1197                   nfieldvector1,&iselect,exoid,num_time_steps,countvars,nout);
1198         countvars+=1;
1199       }
1200     }
1201 
1202     /* storing the imaginary part of the thermal error estimator in the nodes
1203        for the odd modes of cyclic symmetry calculations */
1204   
1205     if(*noddiam>=0){
1206       if(strcmp1(&filab[2784],"HER")==0){
1207         if (countbool==3){
1208           countvars+=1;
1209         }else if(countbool==2){
1210           var_names[countvars++]="HFLSTD-IMG";
1211         }else{
1212           nterms=3;
1213           FORTRAN(errorestimator,(&qfx[3*mi[0]**ne],qfn,ipkon,inum,kon,lakon,nk,ne,
1214                                   mi,ielmat,thicke,&nterms));
1215           
1216           ncomp=1;
1217           ifield[0]=1;
1218           icomp[0]=2;
1219           
1220           exoselect(qfn,qfn,&iset,&nkcoords,inum,istartset,iendset,
1221                     ialset,ngraph,&ncomp,ifield,icomp,
1222                     nfieldvector1,&iselect,exoid,num_time_steps,countvars,nout);
1223           
1224           countvars+=1;
1225         }
1226       }
1227     }
1228 
1229     /* storing the total temperatures in the network nodes */
1230     if(strcmp1(&filab[1131],"TT  ")==0){
1231       if (countbool==3){
1232         countvars+=1;
1233       }else if(countbool==2){
1234         var_names[countvars++]="TT";
1235       }else{
1236         
1237         iselect=-1;
1238         frdset(&filab[1131],set,&iset,istartset,iendset,ialset,
1239                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1240                ngraph);
1241         
1242         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1243                   ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1244                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1245         printf ("Warning: export of TT to exo not tested.\n");
1246         countvars+=1;
1247       }
1248     }
1249 
1250     /* storing the total temperatures in the network nodes */
1251     if(strcmp1(&filab[1218],"MF  ")==0){
1252       if (countbool==3){
1253         countvars+=1;
1254       }else if(countbool==2){
1255         var_names[countvars++]="MF";
1256       }else{
1257         
1258         iselect=-1;
1259         frdset(&filab[1218],set,&iset,istartset,iendset,ialset,
1260                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1261                ngraph);
1262         
1263         icomp[0]=1;
1264         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1265                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1266                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1267         countvars+=1;
1268       }
1269     }
1270 
1271     /* storing the total pressure in the network nodes */
1272     if(strcmp1(&filab[1305],"PT  ")==0){
1273       if (countbool==3){
1274         countvars+=1;
1275       }else if(countbool==2){
1276         var_names[countvars++]="PT";
1277       }else{
1278         
1279         iselect=-1;
1280         frdset(&filab[1305],set,&iset,istartset,iendset,ialset,
1281                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1282                ngraph);
1283         
1284         icomp[0]=2;
1285         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1286                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1287                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1288         countvars+=1;
1289       }
1290     }
1291     
1292     /* storing the static pressure in the liquid network nodes */
1293     if(strcmp1(&filab[1827],"PS  ")==0){
1294       if (countbool==3){
1295         countvars+=1;
1296       }else if(countbool==2){
1297         var_names[countvars++]="STPRESS PS";
1298       }else{
1299         
1300         iselect=-1;
1301         frdset(&filab[1827],set,&iset,istartset,iendset,ialset,
1302                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1303                ngraph);
1304         
1305         icomp[0]=2;
1306         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1307                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1308                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1309         printf ("Warning: export of PS to exo not tested.\n");
1310         countvars+=1;
1311       }    
1312     }
1313 
1314     /* storing the liquid depth in the channel nodes */
1315     if(strcmp1(&filab[2349],"PS  ")==0){
1316       if (countbool==3){
1317         countvars+=1;
1318       }else if(countbool==2){
1319         var_names[countvars++]="DEPTH";
1320       }else{
1321         
1322         iselect=-1;
1323         frdset(&filab[2349],set,&iset,istartset,iendset,ialset,
1324                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1325                ngraph);
1326         
1327         icomp[0]=2;
1328         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1329                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1330                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1331         printf ("Warning: export of DEPTH to exo not tested.\n");
1332         countvars+=1;
1333       }
1334     }
1335 
1336     /* storing the critical depth in the channel nodes */
1337     if(strcmp1(&filab[2436],"HCRI")==0){
1338       if (countbool==3){
1339         countvars+=1;
1340       }else if(countbool==2){
1341         var_names[countvars++]="HCRIT";
1342       }else{
1343 
1344         iselect=-1;
1345         frdset(&filab[2436],set,&iset,istartset,iendset,ialset,
1346                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1347                ngraph);
1348     
1349         icomp[0]=3;
1350         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1351                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1352                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1353         printf ("Warning: export of HCRIT to exo not tested.\n");
1354         countvars+=1;
1355       }
1356     }
1357 
1358     /* storing the static temperature in the network nodes */
1359     if(strcmp1(&filab[1392],"TS  ")==0){
1360       if (countbool==3){
1361         countvars+=1;
1362       }else if(countbool==2){
1363         var_names[countvars++]="STTEMP";
1364       }else{
1365 
1366         iselect=-1;
1367         frdset(&filab[1392],set,&iset,istartset,iendset,ialset,
1368                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1369                ngraph);
1370     
1371         icomp[0]=3;
1372         exoselect(v,v,&iset,&nkcoords,inum,istartset,iendset,
1373                   ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
1374                   nfieldvector0,&iselect,exoid,num_time_steps,countvars,nout);
1375         printf ("Warning: export of STTEMP to exo not tested.\n");
1376         countvars+=1;
1377       }
1378     }
1379 
1380     /*  the remaining lines only apply to frequency calculations
1381         with cyclic symmetry, complex frequency and steady state calculations */
1382   
1383     if((*nmethod!=2)&&(*nmethod!=5)&&(*nmethod!=6)&&(*nmethod!=7)){goto WRITENAMES;}
1384     if((*nmethod==5)&&(*mode==-1)){goto WRITENAMES;}
1385   
1386     /* storing the displacements in the nodes (magnitude, phase) */
1387     if(strcmp1(&filab[870],"PU  ")==0){
1388       if (countbool==3){
1389         countvars+=6;
1390       }else if(countbool==2){
1391         // Note reordered (done in exoselect) relative to frd file
1392         // Order must be xx  yy  zz  xy  xz  yz
1393         var_names[countvars++]="PDISP MAG1";
1394         var_names[countvars++]="PDISP MAG2";
1395         var_names[countvars++]="PDISP MAG3";
1396         var_names[countvars++]="PDISP PHA1";
1397         var_names[countvars++]="PDISP PHA3";
1398         var_names[countvars++]="PDISP PHA2";
1399       }else{
1400         iselect=1;
1401         
1402         frdset(&filab[870],set,&iset,istartset,iendset,ialset,
1403                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1404                ngraph);
1405         
1406         exoselect(vr,vi,&iset,&nkcoords,inum,istartset,iendset,
1407                   ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
1408                   nfieldvectph,&iselect,exoid,num_time_steps,countvars,nout);
1409         countvars+=6;
1410       }
1411     }
1412 
1413     /* storing the temperatures in the nodes (magnitude, phase) */
1414     if(strcmp1(&filab[957],"PNT ")==0){
1415       if (countbool==3){
1416         countvars+=2;
1417       }else if(countbool==2){
1418         var_names[countvars++]="PNDTEMP MAG1";
1419         var_names[countvars++]="PNDTEMP PHA2";
1420       }else{
1421         iselect=1;
1422         
1423         frdset(&filab[957],set,&iset,istartset,iendset,ialset,
1424                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1425                ngraph);
1426     
1427         exoselect(vr,vi,&iset,&nkcoords,inum,istartset,iendset,
1428                   ialset,ngraph,&ncompscalph,ifieldscalph,icompscalph,
1429                   nfieldscalph,&iselect,exoid,num_time_steps,countvars,nout);
1430         printf ("Warning: export of PNDTEMP to exo not tested.\n");
1431         countvars+=6;
1432       }
1433     }
1434 
1435     /* storing the stresses in the nodes (magnitude, phase) */
1436           
1437     if(strcmp1(&filab[1479],"PHS ")==0){
1438       if (countbool==3){
1439         countvars+=12;
1440       }else if(countbool==2){
1441         // Note reordered (done in exoselect) relative to frd file
1442         // Order must be xx  yy  zz  xy  xz  yz
1443         var_names[countvars++]="PSTRESS MAGXX";
1444         var_names[countvars++]="PSTRESS MAGYY";
1445         var_names[countvars++]="PSTRESS MAGZZ";
1446         var_names[countvars++]="PSTRESS MAGXY";
1447         var_names[countvars++]="PSTRESS MAGXZ";
1448         var_names[countvars++]="PSTRESS MAGYZ";
1449         var_names[countvars++]="PSTRESS PHAXX";
1450         var_names[countvars++]="PSTRESS PHAYY";
1451         var_names[countvars++]="PSTRESS PHAZZ";
1452         var_names[countvars++]="PSTRESS PHAXY";
1453         var_names[countvars++]="PSTRESS PHAXZ";
1454         var_names[countvars++]="PSTRESS PHAYZ";
1455       }else{
1456         iselect=1;
1457         
1458         frdset(&filab[1479],set,&iset,istartset,iendset,ialset,
1459                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1460                ngraph);
1461         
1462         exoselect(stnr,stni,&iset,&nkcoords,inum,istartset,iendset,
1463                   ialset,ngraph,&ncomptensph,ifieldtensph,icomptensph,
1464                   nfieldtensph,&iselect,exoid,num_time_steps,countvars,nout);
1465         countvars+=12;
1466       }
1467     }
1468 
1469     /* storing the displacements in the nodes (magnitude, phase) */
1470     if(strcmp1(&filab[2610],"PRF ")==0){
1471       if (countbool==3){
1472         countvars+=6;
1473       }else if(countbool==2){
1474         // Note reordered (done in exoselect) relative to frd file
1475         // Order must be xx  yy  zz  xy  xz  yz
1476         var_names[countvars++]="PFORC MAG1";
1477         var_names[countvars++]="PFORC MAG2";
1478         var_names[countvars++]="PFORC MAG3";
1479         var_names[countvars++]="PFORC PHA1";
1480         var_names[countvars++]="PFORC PHA3";
1481         var_names[countvars++]="PFORC PHA2";
1482       }else{
1483         iselect=1;
1484     
1485         frdset(&filab[2610],set,&iset,istartset,iendset,ialset,
1486                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1487                ngraph);
1488     
1489         exoselect(fnr,fni,&iset,&nkcoords,inum,istartset,iendset,
1490                   ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
1491                   nfieldvectph,&iselect,exoid,num_time_steps,countvars,nout);
1492         printf ("Warning: export of PFORC to exo not tested.\n");
1493         countvars+=6;
1494       }
1495     }
1496 
1497     /* the remaining parts are for frequency calculations with cyclic symmetry only */
1498     if(*nmethod!=2){goto WRITENAMES;}
1499     
1500     /* storing the maximum displacements of the nodes in the base sector
1501        (components, magnitude) */
1502           
1503     if(strcmp1(&filab[1566],"MAXU")==0){
1504       if (countbool==3){
1505         countvars+=4;
1506       }else if(countbool==2){
1507         var_names[countvars++]="MDISP DX";
1508         var_names[countvars++]="MDISP DY";
1509         var_names[countvars++]="MDISP DZ";
1510         var_names[countvars++]="MDISP ANG";
1511       }else{
1512         iselect=1;
1513       
1514         frdset(&filab[1566],set,&iset,istartset,iendset,ialset,
1515                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1516                ngraph);
1517     
1518         ncomp=4;
1519         ifield[0]=1;icomp[0]=1;
1520         ifield[1]=1;icomp[1]=2;
1521         ifield[2]=1;icomp[2]=3;
1522         ifield[3]=1;icomp[3]=0;
1523         nfield[0]=4;nfield[1]=4;
1524 
1525         exoselect(vmax,vmax,&iset,&nkcoords,inum,istartset,iendset,
1526                   ialset,ngraph,&ncomp,ifield,icomp,
1527                   nfield,&iselect,exoid,num_time_steps,countvars,nout);
1528         printf ("Warning: export of MDISP to exo not tested.\n");
1529         countvars+=4;
1530       }
1531     }
1532 
1533     //  /* storing the worst principal stress at the nodes
1534     //     in the basis sector (components, magnitude)
1535     //
1536     //     the worst principal stress is the maximum of the
1537     //     absolute value of all principal stresses, times
1538     //     its original sign */
1539     //    
1540     if(strcmp1(&filab[1653],"MAXS")==0){
1541       if (countbool==3){
1542         countvars+=7;
1543       }else if(countbool==2){
1544         // Note reordered (done in exoselect) relative to frd file
1545         // Order must be xx  yy  zz  xy  xz  yz
1546         var_names[countvars++]="MSTRESS SXX";
1547         var_names[countvars++]="MSTRESS SYY";
1548         var_names[countvars++]="MSTRESS SZZ";
1549         var_names[countvars++]="MSTRESS SXY";
1550         var_names[countvars++]="MSTRESS SXZ";
1551         var_names[countvars++]="MSTRESS SYZ";
1552         var_names[countvars++]="MSTRESS MAG";
1553       }else{
1554         iselect=1;
1555       
1556         frdset(&filab[1653],set,&iset,istartset,iendset,ialset,
1557                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1558                ngraph);
1559     
1560         ncomp=7;
1561         ifield[0]=1;icomp[0]=1;
1562         ifield[1]=1;icomp[1]=2;
1563         ifield[2]=1;icomp[2]=3;
1564         ifield[3]=1;icomp[3]=4;
1565         ifield[4]=1;icomp[4]=6;
1566         ifield[5]=1;icomp[5]=5;
1567         ifield[6]=1;icomp[6]=0;
1568         nfield[0]=7;nfield[1]=7;
1569 
1570         exoselect(stnmax,stnmax,&iset,&nkcoords,inum,istartset,iendset,
1571                   ialset,ngraph,&ncomp,ifield,icomp,
1572                   nfield,&iselect,exoid,num_time_steps,countvars,nout);
1573         printf ("Warning: export of MSTRESS to exo not tested.\n");
1574         countvars+=7;
1575       }
1576     }
1577 
1578     /* storing the worst principal strain at the nodes
1579        in the basis sector (components, magnitude)
1580 
1581        the worst principal strain is the maximum of the
1582        absolute value of all principal strains, times
1583        its original sign */
1584     if(strcmp1(&filab[2523],"MAXE")==0){
1585       if (countbool==3){
1586         countvars+=7;
1587       }else if(countbool==2){
1588         // Note reordered (done in exoselect) relative to frd file
1589         // Order must be xx  yy  zz  xy  xz  yz
1590         var_names[countvars++]="MAXE EXX";
1591         var_names[countvars++]="MAXE EYY";
1592         var_names[countvars++]="MAXE EZZ";
1593         var_names[countvars++]="MAXE EXY";
1594         var_names[countvars++]="MAXE EXZ";
1595         var_names[countvars++]="MAXE EYZ";
1596         var_names[countvars++]="MAXE MAG";
1597       }else{
1598         iselect=1;
1599     
1600         frdset(&filab[2523],set,&iset,istartset,iendset,ialset,
1601                inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1602                ngraph);
1603     
1604         ncomp=7;
1605         ifield[0]=1;icomp[0]=1;
1606         ifield[1]=1;icomp[1]=2;
1607         ifield[2]=1;icomp[2]=3;
1608         ifield[3]=1;icomp[3]=4;
1609         ifield[4]=1;icomp[4]=6;
1610         ifield[5]=1;icomp[5]=5;
1611         ifield[6]=1;icomp[6]=0;
1612         nfield[0]=7;nfield[1]=7;
1613 
1614         exoselect(eenmax,eenmax,&iset,&nkcoords,inum,istartset,iendset,
1615                   ialset,ngraph,&ncomp,ifield,icomp,
1616                   nfield,&iselect,exoid,num_time_steps,countvars,nout);
1617         printf ("Warning: export of MAXE to exo not tested.\n");
1618         countvars+=7;
1619       }
1620     }
1621     
1622   WRITENAMES:
1623     if (countbool==3){
1624       errr = ex_put_var_param (exoid, "n", countvars);
1625       ex_update (exoid);  
1626       // var_names = (char *) calloc (countvars, sizeof (char));
1627       // var_names = (char *) calloc (countvars, MAX_STR_LENGTH);
1628       int ii=0; // 4 bit integer
1629       for (ii=0; ii<countvars; ii++)
1630         {
1631           var_names[ii] = (char *) calloc ((MAX_STR_LENGTH+1), sizeof(char));
1632         }
1633     }else if(countbool==2){
1634       errr = ex_put_var_names (exoid, "n", countvars, var_names);
1635       if (errr) printf ("Unable to update variable names.  Was output data requested?\n");
1636       //  for (i=0; i<countvars; i++)
1637       //        {
1638       //          free(var_names[i]);
1639       //        }
1640       ex_update (exoid);  
1641     }
1642   
1643     countvars=0;
1644     // printf ("Decrement countbool to %i\n",--countbool);
1645     --countbool;
1646   };
1647 
1648   // free (var_names);
1649   free (node_map_inv);
1650   ex_update (exoid);  
1651   ex_close(exoid);
1652   return;
1653 }

/* [<][>][^][v][top][bottom][index][help] */