root/src/complexfreq.c

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

DEFINITIONS

This source file includes following definitions.
  1. complexfreq

   1 /*     CalculiX - A 3-dimensional finite element program                   */
   2 /*              Copyright (C) 1998-2015 Guido Dhondt                          */
   3 /*     This program is free software; you can redistribute it and/or     */
   4 /*     modify it under the terms of the GNU General Public License as    */
   5 /*     published by the Free Software Foundation(version 2);    */
   6 /*                    */
   7 
   8 /*     This program is distributed in the hope that it will be useful,   */
   9 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  10 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  11 /*     GNU General Public License for more details.                      */
  12 
  13 /*     You should have received a copy of the GNU General Public License */
  14 /*     along with this program; if not, write to the Free Software       */
  15 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  16 
  17 #include <stdio.h>
  18 #include <math.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include "CalculiX.h"
  22 
  23 #ifdef SPOOLES
  24    #include "spooles.h"
  25 #endif
  26 #ifdef SGI
  27    #include "sgi.h"
  28 #endif
  29 #ifdef TAUCS
  30    #include "tau.h"
  31 #endif
  32 #ifdef PARDISO
  33    #include "pardiso.h"
  34 #endif
  35 
  36 void complexfreq(double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, 
  37                ITG **nodebounp, ITG **ndirbounp, double **xbounp, ITG *nboun,
  38                ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp,
  39                ITG *nmpc, ITG *nodeforc,ITG *ndirforc,double *xforc, 
  40                ITG *nforc,ITG *nelemload, char *sideload,double *xload,
  41                ITG *nload, 
  42                ITG **nactdofp,ITG *neq, ITG *nzl,ITG *icol, ITG *irow, 
  43                ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG **ikbounp, 
  44                ITG **ilbounp,double *elcon, ITG *nelcon, double *rhcon, 
  45                ITG *nrhcon,double *cocon, ITG *ncocon,
  46                double *alcon, ITG *nalcon, double *alzero, 
  47                ITG **ielmatp,ITG **ielorienp, ITG *norien, double *orab, 
  48                ITG *ntmat_,double **t0p, 
  49                double **t1p,ITG *ithermal,double *prestr, ITG *iprestr, 
  50                double **voldp,ITG *iperturb, double **stip, ITG *nzs, 
  51                double *timepar, double *xmodal,
  52                double **veoldp, char *amname, double *amta,
  53                ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload,
  54                ITG **iamt1p,ITG *jout,
  55                ITG *kode, char *filab,double **emep, double *xforcold, 
  56                double *xloadold,
  57                double **t1oldp, ITG **iambounp, double **xbounoldp, ITG *iexpl,
  58                double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
  59                double *xstate, ITG *npmat_, char *matname, ITG *mi,
  60                ITG *ncmat_, ITG *nstate_, double **enerp, char *jobnamec,
  61                double *ttime, char *set, ITG *nset, ITG *istartset,
  62                ITG *iendset, ITG **ialsetp, ITG *nprint, char *prlab,
  63                char *prset, ITG *nener, double *trab, 
  64                ITG **inotrp, ITG *ntrans, double **fmpcp, char *cbody, 
  65                ITG *ibody,
  66                double *xbody, ITG *nbody, double *xbodyold, ITG *istep,
  67                ITG *isolver,ITG *jq, char *output, ITG *mcs, ITG *nkon,
  68                ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset,
  69                ITG *idrct, ITG *jmax,
  70                double *ctrl, ITG *itpamp, double *tietol,ITG *nalset,
  71                ITG *ikforc, ITG *ilforc, double *thicke,
  72                char *jobnamef,ITG *mei,ITG *nmat,ITG *ielprop,double *prop){
  73 
  74   char fneig[132]="",description[13]="            ",*lakon=NULL,*labmpc=NULL,
  75     *lakont=NULL;
  76 
  77   ITG nev,i,j,k,idof,*inum=NULL,*ipobody=NULL,inewton=0,id,
  78     iinc=0,l,iout=1,ielas,icmd,ifreebody,mode,m,nherm,
  79     *kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,*islavact=NULL,
  80     *inotr=NULL,*nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,
  81     *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
  82     *ilmpc=NULL,nsectors,nmd,nevd,*nm=NULL,*iamt1=NULL,*islavnode=NULL,
  83     ngraph=1,nkg,neg,ne0,ij,lprev,nope,indexe,ilength,*nslavnode=NULL,
  84     *ipneigh=NULL,*neigh=NULL,index,im,cyclicsymmetry,inode,
  85     *ialset=*ialsetp,mt=mi[1]+1,kmin,kmax,i1,
  86     *iter=NULL,lint,lfin,kk,kkv,kk6,kkx,icomplex,igeneralizedforce,
  87     idir,*inumt=NULL,icntrl,imag,jj,is,l1,*inocs=NULL,ml1,l2,nkt,net,
  88     *ipkont=NULL,*ielmatt=NULL,*inotrt=NULL,*kont=NULL,node,iel,*ielcs=NULL,
  89     ielset,*istartnmd=NULL,*iendnmd=NULL,inmd,neqact,*nshcon=NULL,
  90     *ipev=NULL,icfd=0,*inomat=NULL,mortar=0,*islavsurf=NULL;
  91 
  92   long long i2;
  93 
  94   double *d=NULL, *z=NULL,*stiini=NULL,*cc=NULL,*v=NULL,*zz=NULL,*emn=NULL,
  95     *stn=NULL, *stx=NULL, *een=NULL, *adb=NULL,*xstiff=NULL,*cdn=NULL,
  96     *aub=NULL,*f=NULL, *fn=NULL,*epn=NULL,*xstateini=NULL,
  97     *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
  98     *qfx=NULL, *cgr=NULL, *au=NULL,dtime,reltime,*t0=NULL,*t1=NULL,*t1old=NULL,
  99     sum,qa[3],cam[5],accold[1],bet,gam,*ad=NULL,alpham,betam,
 100     *co=NULL,*xboun=NULL,*xbounold=NULL,*vold=NULL,*emeini=NULL,
 101     *eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,*veold=NULL,
 102     *adc=NULL,*auc=NULL,*zc=NULL,*fnr=NULL,*fni=NULL,setnull,deltmx,dd,
 103     theta,*vini=NULL,*vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,*vmax=NULL,
 104     *stnmax=NULL,*cstr=NULL,*sti=*stip,time0=0.0,time=0.0,zero=0.0,
 105     *springarea=NULL,*eenmax=NULL,*aa=NULL,*bb=NULL,*xx=NULL,
 106     *eiga=NULL,*eigb=NULL,*eigxx=NULL,*temp=NULL,*coefmpcnew=NULL,xreal,
 107     ximag,t[3],*vt=NULL,*t1t=NULL,*stnt=NULL,*eent=NULL,*fnt=NULL,*enernt=NULL,
 108     *stxt=NULL,pi,ctl,stl,*cot=NULL,*qfnt=NULL,vreal,vimag,constant,stnreal,
 109     stnimag,freq,*emnt=NULL,*shcon=NULL,*eig=NULL,*clearini=NULL,
 110     *eigxr=NULL,*eigxi=NULL,*xmac=NULL,*bett=NULL,*betm=NULL,*xmaccpx=NULL,
 111     fmin=0.,fmax=1.e30,*xmr=NULL,*xmi=NULL,*zi=NULL,*eigx=NULL,
 112     *pslavsurf=NULL,*pmastsurf=NULL,*cdnr=NULL,*cdni=NULL,*tinc,*tper,*tmin,*tmax;
 113 
 114   FILE *f1;
 115 
 116 #ifdef SGI
 117   ITG token;
 118 #endif
 119 
 120   pi=4.*atan(1.);
 121   constant=180./pi;
 122 
 123   co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
 124   ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
 125   ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;
 126   xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
 127   vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
 128   coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 129   fmpc=*fmpcp;veold=*veoldp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
 130 
 131   tinc=&timepar[0];
 132   tper=&timepar[1];
 133   tmin=&timepar[2];
 134   tmax=&timepar[3];
 135 
 136   if(ithermal[0]<=1){
 137       kmin=1;kmax=3;
 138   }else if(ithermal[0]==2){
 139       kmin=0;kmax=mi[1];if(kmax>2)kmax=2;
 140   }else{
 141       kmin=0;kmax=3;
 142   }
 143 
 144   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 145 
 146   dtime=*tinc;
 147 
 148   alpham=xmodal[0];
 149   betam=xmodal[1];
 150 
 151   dd=ctrl[16];deltmx=ctrl[26];
 152 
 153   /* determining nzl */
 154 
 155   *nzl=0;
 156   for(i=neq[1];i>0;i--){
 157       if(icol[i-1]>0){
 158           *nzl=i;
 159           break;
 160       }
 161   }
 162 
 163   /* opening the eigenvalue file and checking for cyclic symmetry */
 164 
 165   strcpy(fneig,jobnamec);
 166   strcat(fneig,".eig");
 167 
 168   if((f1=fopen(fneig,"rb"))==NULL){
 169     printf("*ERROR in complexfreq: cannot open eigenvalue file for reading");
 170     exit(0);
 171   }
 172 
 173   printf(" *INFO  in complexfreq: if there are problems reading the .eig file this may be due to:\n");
 174   printf("        1) the nonexistence of the .eig file\n");
 175   printf("        2) other boundary conditions than in the input deck\n");
 176   printf("           which created the .eig file\n\n");
 177 
 178   if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
 179       printf("*ERROR in complexfreq reading the cyclic symmetry flag in the eigenvalue file");
 180       exit(0);
 181   }
 182 
 183   if(fread(&nherm,sizeof(ITG),1,f1)!=1){
 184       printf("*ERROR in complexfreq reading the Hermitian flag in the eigenvalue file");
 185       exit(0);
 186   }
 187 
 188   if(nherm!=1){
 189       printf("*ERROR in complexfreq: the eigenvectors in the .eig-file result\n");
 190       printf("       from a non-Hermitian eigenvalue problem. The complex\n");
 191       printf("       frequency procedure cannot handle that yet\n\n");
 192       FORTRAN(stop,());
 193   }
 194 
 195   nsectors=1;
 196 
 197   if(!cyclicsymmetry){
 198 
 199       nkg=*nk;
 200       neg=*ne;
 201 
 202       if(fread(&nev,sizeof(ITG),1,f1)!=1){
 203           printf("*ERROR in complexfreq reading the number of eigenvalues in the eigenvalue file...");
 204           exit(0);
 205       }
 206       
 207 
 208       if(nherm==1){
 209           NNEW(d,double,nev);
 210           if(fread(d,sizeof(double),nev,f1)!=nev){
 211               printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
 212               exit(0);
 213           }
 214       }else{
 215           NNEW(d,double,2*nev);
 216           if(fread(d,sizeof(double),2*nev,f1)!=2*nev){
 217               printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
 218               exit(0);
 219           }
 220       }
 221       
 222       NNEW(ad,double,neq[1]);
 223       NNEW(adb,double,neq[1]);
 224       NNEW(au,double,nzs[2]);
 225       NNEW(aub,double,nzs[1]);
 226       
 227       /* reading the stiffness matrix */
 228 
 229       if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
 230           printf("*ERROR in complexfreq reading the diagonal of the stiffness matrix in the eigenvalue file...");
 231           exit(0);
 232       }
 233       
 234       if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
 235           printf("*ERROR in complexfreq reading the off-diagonal terms of the stiffness matrix in the eigenvalue file...");
 236           exit(0);
 237       }
 238       
 239       /* reading the mass matrix */
 240 
 241       if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 242           printf("*ERROR in complexfreq reading the diagonal of the mass matrix in eigenvalue file...");
 243           exit(0);
 244       }
 245       
 246       if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 247           printf("*ERROR in complexfreq reading the off-diagonals of the mass matrix in the eigenvalue file...");
 248           exit(0);
 249       }
 250       
 251       /* reading the eigenvectors */
 252 
 253       if(nherm==1){
 254           NNEW(z,double,neq[1]*nev);
 255           if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
 256               printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
 257               exit(0);
 258           }
 259       }else{
 260           NNEW(z,double,2*neq[1]*nev);
 261           if(fread(z,sizeof(double),2*neq[1]*nev,f1)!=2*neq[1]*nev){
 262               printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
 263               exit(0);
 264           }
 265       }
 266 
 267       NNEW(nm,ITG,nev);
 268       for(i=0;i<nev;i++){nm[i]=-1;}
 269   }
 270   else{
 271 
 272       if(*nmethod==6){
 273         printf("*ERROR in complexfreq: Coriolis forces cannot\n");
 274         printf("       be combined with cyclic symmetry\n\n");
 275         FORTRAN(stop,());
 276       }
 277 
 278       nev=0;
 279       do{
 280           if(fread(&nmd,sizeof(ITG),1,f1)!=1){
 281               break;
 282           }
 283           if(fread(&nevd,sizeof(ITG),1,f1)!=1){
 284               printf("*ERROR in complexfreq reading the number of eigenvalues in the eigenvalue file...");
 285               exit(0);
 286               }
 287           if(nev==0){
 288               if(nherm==1){NNEW(d,double,nevd);
 289               }else{NNEW(d,double,2*nevd);}
 290               NNEW(nm,ITG,nevd);
 291           }else{
 292               printf("*ERROR in complexfreq: flutter forces cannot\n");
 293               printf("       be combined with multiple modal diameters\n");
 294               printf("       in cyclic symmetry calculations\n\n");
 295               FORTRAN(stop,());
 296           }
 297           
 298           if(nherm==1){
 299               if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
 300                   printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
 301                   exit(0);
 302               }
 303           }else{
 304               if(fread(&d[nev],sizeof(double),2*nevd,f1)!=2*nevd){
 305                   printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
 306                   exit(0);
 307               }
 308           }
 309           for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
 310           
 311           if(nev==0){
 312               NNEW(adb,double,neq[1]);
 313               NNEW(aub,double,nzs[1]);
 314 
 315               if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 316                   printf("*ERROR in complexfreq reading the diagonal of the mass matrix in the eigenvalue file...");
 317                   exit(0);
 318               }
 319               
 320               if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 321                   printf("*ERROR in complexfreq reading the off-diagonals of the mass matrix in the eigenvalue file...");
 322                   exit(0);
 323               }
 324           }
 325           
 326           if(nev==0){
 327               NNEW(z,double,neq[1]*nevd);
 328           }else{
 329               RENEW(z,double,(long long)neq[1]*(nev+nevd));
 330           }
 331           
 332           if(fread(&z[neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
 333               printf("*ERROR in complexfreq reading eigenvectors in the eigenvalue file...");
 334               exit(0);
 335           }
 336           nev+=nevd;
 337       }while(1);
 338 
 339       /* removing double eigenmodes */
 340 
 341       j=-1;
 342       for(i=0;i<nev;i++){
 343           if((i/2)*2==i){
 344               j++;
 345               if(nherm==1){
 346                   d[j]=d[i];
 347               }else{
 348                   d[2*j]=d[2*i];d[2*j+1]=d[2*i+1];
 349               }
 350               nm[j]=nm[i];
 351               for(k=0;k<neq[1];k++){
 352                   z[j*neq[1]+k]=z[i*neq[1]+k];
 353               }
 354           }
 355       }
 356       nev=j+1;
 357       if(nherm==1){RENEW(d,double,nev);}else{RENEW(d,double,2*nev);}
 358       RENEW(nm,ITG,nev);
 359       RENEW(z,double,neq[1]*nev);
 360 
 361       /* determining the maximum amount of segments */
 362 
 363       for(i=0;i<*mcs;i++){
 364           if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
 365       }
 366 
 367         /* determining the maximum number of sectors to be plotted */
 368 
 369       for(j=0;j<*mcs;j++){
 370           if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
 371       }
 372       nkg=*nk*ngraph;
 373       neg=*ne*ngraph;
 374 
 375   }
 376 
 377   fclose(f1);
 378 
 379   /* assigning nodes and elements to sectors */
 380 
 381   if(cyclicsymmetry){
 382     NNEW(inocs,ITG,*nk);
 383     NNEW(ielcs,ITG,*ne);
 384     ielset=cs[12];
 385     if((*mcs!=1)||(ielset!=0)){
 386       for(i=0;i<*nk;i++) inocs[i]=-1;
 387       for(i=0;i<*ne;i++) ielcs[i]=-1;
 388     }
 389     
 390     for(i=0;i<*mcs;i++){
 391       is=cs[17*i+4];
 392       if((is==1)&&(*mcs==1)) continue;
 393       ielset=cs[17*i+12];
 394       if(ielset==0) continue;
 395       for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
 396         if(ialset[i1]>0){
 397           iel=ialset[i1]-1;
 398           if(ipkon[iel]<0) continue;
 399           ielcs[iel]=i;
 400           indexe=ipkon[iel];
 401           if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 402           else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 403           else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 404           else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 405           else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 406           else {nope=6;}
 407           for(i2=0;i2<nope;++i2){
 408             node=kon[indexe+i2]-1;
 409             inocs[node]=i;
 410           }
 411         }
 412         else{
 413           iel=ialset[i1-2]-1;
 414           do{
 415             iel=iel-ialset[i1];
 416             if(iel>=ialset[i1-1]-1) break;
 417             if(ipkon[iel]<0) continue;
 418             ielcs[iel]=i;
 419             indexe=ipkon[iel];
 420             if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 421             else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 422             else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 423             else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 424             else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 425             else {nope=6;}
 426             for(i2=0;i2<nope;++i2){
 427               node=kon[indexe+i2]-1;
 428               inocs[node]=i;
 429             }
 430           }while(1);
 431         }
 432       } 
 433     }
 434   }
 435   
 436   /* check for rigid body modes 
 437      if there is a jump of 1.e4 in two subsequent eigenvalues
 438      all eigenvalues preceding the jump are considered to
 439      be rigid body modes and their frequency is set to zero */
 440 
 441   if(nherm==1){
 442       setnull=1.;
 443       for(i=nev-2;i>-1;i--){
 444           if(fabs(d[i])<0.0001*fabs(d[i+1])) setnull=0.;
 445           d[i]*=setnull;
 446       }
 447   }else{
 448       setnull=1.;
 449       for(i=nev-2;i>-1;i--){
 450           if(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])<
 451                  0.0001*sqrt(d[2*i+2]*d[2*i+2]+d[2*i+3]*d[2*i+3])) setnull=0.;
 452           d[2*i]*=setnull;
 453           d[2*i+1]*=setnull;
 454       }
 455   }
 456 
 457   /* determining the frequency ranges corresponding to one
 458      and the same nodal diameter */
 459 
 460   if(cyclicsymmetry){
 461       NNEW(istartnmd,ITG,nev);
 462       NNEW(iendnmd,ITG,nev);
 463       nmd=0;
 464       inmd=nm[0];
 465       istartnmd[0]=1;
 466       for(i=1;i<nev;i++){
 467           if(nm[i]==inmd) continue;
 468           iendnmd[nmd]=i;
 469           nmd++;
 470           istartnmd[nmd]=i+1;
 471           inmd=nm[i];
 472       }
 473       iendnmd[nmd]=nev;
 474       nmd++;
 475       RENEW(istartnmd,ITG,nmd);
 476       RENEW(iendnmd,ITG,nmd);
 477   }
 478 
 479   if(*nmethod==6){
 480 
 481       /* Coriolis */
 482 
 483       neqact=neq[1];
 484 
 485   /* assigning the body forces to the elements */ 
 486 
 487       ifreebody=*ne+1;
 488       NNEW(ipobody,ITG,2**ne);
 489       for(k=1;k<=*nbody;k++){
 490           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 491                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 492           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 493       }
 494       RENEW(ipobody,ITG,2*(ifreebody-1));
 495       
 496       if(cyclicsymmetry){
 497           printf("*ERROR in complexfreq: dashpots are not allowed in combination with cyclic symmetry\n");
 498           FORTRAN(stop,());
 499       }
 500 
 501       NNEW(adc,double,neq[1]);
 502       NNEW(auc,double,nzs[1]);
 503       FORTRAN(mafillcorio,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
 504               xboun,nboun,
 505               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 506               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 507               adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
 508               ikmpc,ilmpc,ikboun,ilboun,
 509               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 510               ielorien,norien,orab,ntmat_,
 511               t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 512               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 513               xstiff,npmat_,&dtime,matname,mi,ncmat_,
 514               ttime,&time0,istep,&iinc,ibody));
 515 
 516       /*  zc = damping matrix * eigenmodes */
 517 
 518       NNEW(zc,double,neq[1]*nev);
 519       for(i=0;i<nev;i++){
 520           FORTRAN(op_corio,(&neq[1],&z[i*neq[1]],&zc[i*neq[1]],adc,auc,
 521           jq,irow));
 522       }
 523       SFREE(adc);SFREE(auc);
 524 
 525       /* cc is the reduced damping matrix (damping matrix mapped onto
 526          space spanned by eigenmodes) */
 527 
 528       NNEW(cc,double,nev*nev);
 529       for(i=0;i<nev;i++){
 530           for(j=0;j<=i;j++){
 531               for(k=0;k<neq[1];k++){
 532                   cc[i*nev+j]+=z[j*neq[1]+k]*zc[i*neq[1]+k];
 533               }
 534           }
 535       }
 536 
 537       /* symmetric part of cc matrix */
 538 
 539       for(i=0;i<nev;i++){
 540           for(j=i+1;j<nev;j++){
 541               cc[i*nev+j]=-cc[j*nev+i];
 542           }
 543       }
 544       SFREE(zc);
 545 
 546       /* solving for the complex eigenvalues */
 547 
 548       NNEW(aa,double,4*nev*nev);
 549       NNEW(bb,double,4*nev*nev);
 550       NNEW(xx,double,4*nev*nev);
 551       NNEW(temp,double,4*nev*nev);
 552       NNEW(eiga,double,2*nev);
 553       NNEW(eigb,double,2*nev);
 554       NNEW(eigxx,double,2*nev);
 555       NNEW(iter,ITG,nev);
 556 
 557       FORTRAN(coriolissolve,(cc,&nev,aa,bb,xx,eiga,eigb,eigxx,
 558               iter,d,temp));
 559       
 560       SFREE(aa);SFREE(bb);SFREE(temp);SFREE(eiga);SFREE(eigb);SFREE(iter);SFREE(cc);
 561 
 562   }else{
 563 
 564       /* flutter */
 565 
 566       /* complex force is being read (e.g. due to fluid flow) */
 567 
 568       if(cyclicsymmetry){
 569           neqact=neq[1]/2;
 570       }else{
 571           neqact=neq[1];
 572       }
 573       NNEW(zc,double,2*neqact*nev);
 574       NNEW(aa,double,4*nev*nev);
 575 
 576       FORTRAN(readforce,(zc,&neqact,&nev,nactdof,ikmpc,nmpc,
 577                          ipompc,nodempc,mi,coefmpc,jobnamef,
 578                          aa,&igeneralizedforce));
 579 
 580       NNEW(bb,double,4*nev*nev);
 581       NNEW(xx,double,4*nev*nev);
 582       NNEW(eiga,double,2*nev);
 583       NNEW(eigb,double,2*nev);
 584       NNEW(eigxx,double,2*nev);
 585       NNEW(iter,ITG,nev);
 586       FORTRAN(forcesolve,(zc,&nev,aa,bb,xx,eiga,eigb,eigxx,iter,d,
 587                 &neq[1],z,istartnmd,iendnmd,&nmd,&cyclicsymmetry,
 588                 &neqact,&igeneralizedforce));
 589       SFREE(aa);SFREE(bb);SFREE(eiga);SFREE(eigb);SFREE(iter);SFREE(zc);
 590       
 591  }
 592 
 593 /* sorting the eigenvalues and eigenmodes according to the size of the
 594    eigenvalues */
 595       
 596   NNEW(ipev,ITG,nev);
 597   NNEW(eigxr,double,nev);
 598   NNEW(aa,double,2*nev);
 599   NNEW(bb,double,4*nev*nev);
 600   
 601   FORTRAN(sortev,(&nev,&nmd,eigxx,&cyclicsymmetry,xx,
 602                   eigxr,ipev,istartnmd,iendnmd,aa,bb));
 603   
 604   SFREE(ipev);SFREE(eigxr);SFREE(aa);SFREE(bb);
 605 
 606   /* storing the eigenvalues in the .dat file */
 607 
 608   if(cyclicsymmetry){
 609       FORTRAN(writeevcscomplex,(eigxx,&nev,nm,&fmin,&fmax));
 610   }else{
 611       FORTRAN(writeevcomplex,(eigxx,&nev,&fmin,&fmax));
 612   }
 613 
 614   /* storing the participation factors */
 615 
 616   NNEW(eigxr,double,nev);
 617   NNEW(eigxi,double,nev);
 618   if(nherm==1){
 619       NNEW(eig,double,nev);
 620       for(l=0;l<nev;l++){
 621           if(d[l]<0.){
 622               eig[l]=0.;
 623           }else{
 624               eig[l]=sqrt(d[l]);
 625           }
 626       }
 627   }else{
 628 
 629       NNEW(eig,double,2*nev);
 630       for(l=0;l<nev;l++){
 631           eig[2*l]=sqrt(sqrt(d[2*l]*d[2*l]+d[2*l+1]*d[2*l+1])+d[2*l])/sqrt(2.);
 632           eig[2*l+1]=sqrt(sqrt(d[2*l]*d[2*l]+d[2*l+1]*d[2*l+1])-d[2*l])/sqrt(2.);
 633           if(d[2*l+1]<0.) eig[2*l+1]=-eig[2*l+1];
 634       }
 635   }
 636   for(l=0;l<nev;l++){
 637       mode=l+1;
 638       for(k=0;k<nev;k++){
 639           eigxr[k]=xx[2*l*nev+2*k];
 640           eigxi[k]=xx[2*l*nev+2*k+1];
 641       }
 642       FORTRAN(writepf,(eig,eigxr,eigxi,&zero,&nev,&mode,&nherm));
 643   }
 644   SFREE(eigxr);SFREE(eigxi);SFREE(eig);SFREE(d);
 645 
 646   if(cyclicsymmetry){
 647       
 648        /* assembling the new eigenmodes */
 649 
 650       /* storage in zz: per eigenmode first the complete real part of
 651          the eigenvector, then the complete imaginary part */
 652       
 653       NNEW(zz,double,(long long)2*nev*neqact);
 654       for(l=0;l<nev;l++){
 655           for(i=0;i<neqact;i++){
 656               for(k=0;k<nev;k++){
 657 
 658                   /* real part */
 659 
 660                   zz[(long long)2*l*neqact+i]+=
 661                       (xx[2*l*nev+2*k]*z[(long long)2*k*neqact+i]
 662                        -xx[2*l*nev+2*k+1]*z[(long long)(2*k+1)*neqact+i]);
 663 
 664                   /* imaginary part */
 665 
 666                   zz[(long long)(2*l+1)*neqact+i]+=
 667                       (xx[2*l*nev+2*k]*z[(long long)(2*k+1)*neqact+i]
 668                        +xx[2*l*nev+2*k+1]*z[(long long)2*k*neqact+i]);
 669               }
 670           }
 671       }
 672 
 673       /* calculating the scalar product of all old eigenmodes with
 674          all new eigenmodes => nev x nev matrix */
 675 
 676       NNEW(xmac,double,nev*nev);
 677       NNEW(xmaccpx,double,4*nev*nev);
 678       NNEW(bett,double,nev);
 679       NNEW(betm,double,nev);
 680       FORTRAN(calcmac,(&neq[1],z,zz,&nev,xmac,xmaccpx,istartnmd,
 681                        iendnmd,&nmd,&cyclicsymmetry,&neqact,bett,betm));
 682       FORTRAN(writemaccs,(xmac,&nev,nm));
 683 
 684       SFREE(xmac);SFREE(bett);SFREE(betm);SFREE(xmaccpx);
 685       SFREE(z);
 686       
 687       /* normalizing the eigenmodes */
 688       
 689       NNEW(z,double,neq[1]);
 690       for(l=0;l<nev;l++){
 691           sum=0.;
 692           DMEMSET(z,0,neq[1],0.);
 693           FORTRAN(op,(&neq[1],&zz[l*neq[1]],z,adb,aub,jq,irow));
 694           for(k=0;k<neq[1];k++){
 695               sum+=zz[l*neq[1]+k]*z[k];
 696           }
 697           
 698           sum=sqrt(sum);
 699           for(k=0;k<neq[1];k++){
 700               zz[l*neq[1]+k]/=sum;
 701           }
 702       }
 703       SFREE(z);
 704 
 705       /* calculating the mass-weighted internal products (eigenvectors are not 
 706          necessarily orthogonal, since the matrix of the eigenvalue problem is
 707          not necessarily Hermitian)
 708          = orthogonality matrices */
 709 
 710       if(mei[3]==1){
 711           
 712           NNEW(xmr,double,nev*nev);
 713           NNEW(xmi,double,nev*nev);
 714           NNEW(z,double,neq[1]);
 715           NNEW(zi,double,neq[1]);
 716           
 717           for(l=0;l<nev;l++){
 718               DMEMSET(z,0,neq[1],0.);
 719               FORTRAN(op,(&neq[1],&zz[l*neq[1]],z,adb,aub,jq,irow));
 720               for(m=l;m<nev;m++){
 721                   for(k=0;k<neq[1];k++){
 722                       xmr[l*nev+m]+=zz[m*neq[1]+k]*z[k];
 723                   }
 724               }
 725               
 726               memcpy(&zi[0],&zz[(2*l+1)*neqact],sizeof(double)*neqact);
 727               for(k=0;k<neqact;k++){zi[neqact+k]=-zz[2*l*neqact+k];}
 728               DMEMSET(z,0,neq[1],0.);
 729               FORTRAN(op,(&neq[1],zi,z,adb,aub,jq,irow));
 730               for(m=l;m<nev;m++){
 731                   for(k=0;k<neq[1];k++){
 732                       xmi[l*nev+m]+=zz[m*neq[1]+k]*z[k];
 733               }
 734               }
 735           }
 736           
 737           /* Hermitian part of the matrix */
 738           
 739           for(l=0;l<nev;l++){
 740               for(m=0;m<l;m++){
 741                   xmr[l*nev+m]=xmr[m*nev+l];
 742                   xmi[l*nev+m]=-xmi[m*nev+l];
 743               }
 744           }
 745           
 746           for(l=0;l<nev;l++){
 747               for(m=0;m<nev;m++){
 748                   printf(" %f",xmr[m*nev+l]);
 749           }
 750               printf("\n");
 751           }
 752           printf("\n");
 753           for(l=0;l<nev;l++){
 754               for(m=0;m<nev;m++){
 755                   printf(" %f",xmi[m*nev+l]);
 756               }
 757               printf("\n");
 758           }
 759           SFREE(z);SFREE(zi);
 760       }
 761 
 762   }else{
 763       
 764       /* no cyclic symmmetry */
 765       
 766       /* assembling the new eigenmodes */
 767       
 768       NNEW(zz,double,2*nev*neq[1]);
 769       for(l=0;l<nev;l++){
 770           for(j=0;j<2;j++){
 771               for(i=0;i<neq[1];i++){
 772                   for(k=0;k<nev;k++){
 773                       zz[(2*l+j)*neq[1]+i]+=xx[2*l*nev+2*k+j]*
 774                             z[(long long)k*neq[1]+i];
 775                   }
 776               }
 777           }
 778       }
 779 
 780       /* calculating the scalar product of all old eigenmodes with
 781          all new eigenmodes => nev x nev matrix */
 782 
 783       NNEW(xmac,double,nev*nev);
 784       NNEW(xmaccpx,double,4*nev*nev);
 785       NNEW(bett,double,nev);
 786       NNEW(betm,double,nev);
 787       FORTRAN(calcmac,(&neq[1],z,zz,&nev,xmac,xmaccpx,istartnmd,
 788                      iendnmd,&nmd,&cyclicsymmetry,&neqact,bett,betm));
 789       FORTRAN(writemac,(xmac,&nev));
 790       SFREE(xmac);SFREE(bett);SFREE(betm);SFREE(xmaccpx);
 791 
 792       SFREE(z);
 793       
 794       /* normalizing the eigenmodes */
 795       
 796       NNEW(z,double,neq[1]);
 797       for(l=0;l<nev;l++){
 798           sum=0.;
 799           
 800           /* Ureal^T*M*Ureal */
 801           
 802           DMEMSET(z,0,neq[1],0.);
 803           FORTRAN(op,(&neq[1],&zz[2*l*neq[1]],z,adb,aub,jq,irow));
 804           for(k=0;k<neq[1];k++){
 805               sum+=zz[2*l*neq[1]+k]*z[k];
 806           }
 807           
 808           /* Uimag^T*M*Uimag */
 809           
 810           DMEMSET(z,0,neq[1],0.);
 811           FORTRAN(op,(&neq[1],&zz[(2*l+1)*neq[1]],z,adb,aub,jq,irow));
 812           for(k=0;k<neq[1];k++){
 813               sum+=zz[(2*l+1)*neq[1]+k]*z[k];
 814           }
 815           
 816           sum=sqrt(sum);
 817           for(k=0;k<2*neq[1];k++){
 818               zz[2*l*neq[1]+k]/=sum;
 819           }
 820       }
 821       SFREE(z);
 822 
 823       /* calculating the mass-weighted internal products (eigenvectors are not 
 824          necessarily orthogonal, since the matrix of the eigenvalue problem is
 825          not necessarily symmetric)
 826          = orthogonality matrices */
 827 
 828       if(mei[3]==1){
 829           
 830           NNEW(xmr,double,nev*nev);
 831           NNEW(xmi,double,nev*nev);
 832           NNEW(z,double,neq[1]);
 833           
 834           for(l=0;l<nev;l++){
 835               sum=0.;
 836               
 837               /* M*Ureal */
 838               
 839               DMEMSET(z,0,neq[1],0.);
 840               FORTRAN(op,(&neq[1],&zz[2*l*neq[1]],z,adb,aub,jq,irow));
 841               
 842               /* Ureal^T*M*Ureal and Uimag^T*M*Ureal */
 843               
 844               for(m=l;m<nev;m++){
 845                   for(k=0;k<neq[1];k++){
 846                       xmr[l*nev+m]+=zz[2*m*neq[1]+k]*z[k];
 847                   }
 848                   for(k=0;k<neq[1];k++){
 849                       xmi[l*nev+m]-=zz[(2*m+1)*neq[1]+k]*z[k];
 850                   }
 851               }
 852               
 853               /* M*Uimag */
 854               
 855               DMEMSET(z,0,neq[1],0.);
 856               FORTRAN(op,(&neq[1],&zz[(2*l+1)*neq[1]],z,adb,aub,jq,irow));
 857               
 858               /* Ureal^T*M*Uimag and Uimag^T*M*Uimag */
 859               
 860               for(m=l;m<nev;m++){
 861                   for(k=0;k<neq[1];k++){
 862                       xmr[l*nev+m]+=zz[(2*m+1)*neq[1]+k]*z[k];
 863                   }
 864                   for(k=0;k<neq[1];k++){
 865                       xmi[l*nev+m]+=zz[2*m*neq[1]+k]*z[k];
 866                   }
 867               }
 868           }
 869           
 870           /* Hermitian part of the matrix */
 871           
 872           for(l=0;l<nev;l++){
 873               for(m=0;m<l;m++){
 874                   xmr[l*nev+m]=xmr[m*nev+l];
 875                   xmi[l*nev+m]=-xmi[m*nev+l];
 876               }
 877           }
 878           
 879           for(l=0;l<nev;l++){
 880               for(m=0;m<nev;m++){
 881                   printf(" %f",xmr[m*nev+l]);
 882               }
 883               printf("\n");
 884           }
 885           printf("\n");
 886           for(l=0;l<nev;l++){
 887               for(m=0;m<nev;m++){
 888                   printf(" %f",xmi[m*nev+l]);
 889               }
 890               printf("\n");
 891           }
 892           SFREE(z);
 893       }
 894 
 895   }
 896 
 897  /*storing new eigenmodes and eigenvalues to *.eig-file for later use in 
 898  steady states dynamic analysis*/
 899 
 900   if(mei[3]==1){
 901 
 902       nherm=0;
 903       
 904       if(!cyclicsymmetry){
 905           if((f1=fopen(fneig,"wb"))==NULL){
 906               printf("*ERROR in complexfreq: cannot open eigenvalue file for writing...");
 907               exit(0);
 908           }
 909           
 910           /* storing a zero as indication that this was not a
 911              cyclic symmetry calculation */
 912           
 913           if(fwrite(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
 914               printf("*ERROR in complexfreq saving the cyclic symmetry flag to the eigenvalue file...");
 915               exit(0);
 916           }
 917           
 918           /* not Hermitian */
 919           
 920           if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
 921               printf("*ERROR in complexfreq saving the Hermitian flag to the eigenvalue file...");
 922               exit(0);
 923           }
 924           
 925           /* storing the number of eigenvalues */
 926           
 927           if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
 928               printf("*ERROR in complexfreq saving the number of eigenfrequencies to the eigenvalue file...");
 929               exit(0);
 930           }
 931           
 932           /* the eigenfrequencies are stored as (radians/time)**2
 933              squaring the complexe eigenvalues first */
 934           
 935           NNEW(eigx,double,2*nev);
 936           for(i=0;i<nev;i++){
 937               eigx[2*i]=eigxx[2*i]*eigxx[2*i]-eigxx[2*i+1]*eigxx[2*i+1];
 938               eigx[2*i+1]=2.*eigxx[2*i]*eigxx[2*i+1];
 939           }
 940           
 941           if(fwrite(eigx,sizeof(double),2*nev,f1)!=2*nev){
 942               printf("*ERROR in complexfreq saving the eigenfrequencies to the eigenvalue file...");
 943               exit(0);
 944           }
 945           
 946           SFREE(eigx);
 947           
 948           /* storing the stiffness matrix */
 949           
 950           if(fwrite(ad,sizeof(double),neq[1],f1)!=neq[1]){
 951               printf("*ERROR in complexfreq saving the diagonal of the stiffness matrix to the eigenvalue file...");
 952               exit(0);
 953           }
 954           if(fwrite(au,sizeof(double),nzs[2],f1)!=nzs[2]){
 955               printf("*ERROR in complexfreq saving the off-diagonal entries of the stiffness matrix to the eigenvalue file...");
 956               exit(0);
 957           }
 958           
 959           /* storing the mass matrix */
 960           
 961           if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
 962               printf("*ERROR in complexfreq saving the diagonal of the mass matrix to the eigenvalue file...");
 963               exit(0);
 964           }
 965           if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 966               printf("*ERROR in complexfreq saving the off-diagonal entries of the mass matrix to the eigenvalue file...");
 967               exit(0);
 968           }
 969           
 970           /* storing the eigenvectors */
 971           
 972           lfin=0;
 973           lint=0;
 974           for(j=0;j<nev;++j){
 975               lint=lfin;
 976               lfin=lfin+2*neq[1];
 977               if(fwrite(&zz[lint],sizeof(double),2*neq[1],f1)!=2*neq[1]){
 978                   printf("*ERROR in complexfreq saving the eigenvectors to the eigenvalue file...");
 979                   exit(0);
 980               }
 981           }
 982           
 983           /* storing the orthogonality matrices */
 984           
 985           if(fwrite(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
 986               printf("*ERROR in complexfreq saving the real orthogonality matrix to the eigenvalue file...");
 987               exit(0);
 988           }
 989           
 990           if(fwrite(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
 991               printf("*ERROR in complexfreq saving the imaginary orthogonality matrix to the eigenvalue file...");
 992               exit(0);
 993           }
 994           
 995       }else{
 996           
 997           /* opening .eig file for writing */
 998           
 999           if((f1=fopen(fneig,"wb"))==NULL){
1000               printf("*ERROR in complexfreq: cannot open eigenvalue file for writing...");
1001               exit(0);
1002           }
1003           /* storing a one as indication that this was a
1004              cyclic symmetry calculation */
1005           
1006           if(fwrite(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
1007               printf("*ERROR in complexfreq saving the cyclic symmetry flag to the eigenvalue file...");
1008               exit(0);
1009           }
1010           
1011           /* not Hermitian */
1012           
1013           if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
1014               printf("*ERROR in complexfreq saving the Hermitian flag to the eigenvalue file...");
1015               exit(0);
1016           }
1017           
1018           if(fwrite(&nm[0],sizeof(ITG),1,f1)!=1){
1019               printf("*ERROR in complexfreq saving the nodal diameter to the eigenvalue file...");
1020               exit(0);
1021           }
1022           if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
1023               printf("*ERROR in complexfreq saving the number of eigenvalues to the eigenvalue file...");
1024               exit(0);
1025           }
1026           
1027           /* the eigenfrequencies are stored as (radians/time)**2
1028              squaring the complexe eigenvalues first */
1029           
1030           NNEW(eigx,double,2*nev);
1031           for(i=0;i<nev;i++){
1032               eigx[2*i]=eigxx[2*i]*eigxx[2*i]-eigxx[2*i+1]*eigxx[2*i+1];
1033               eigx[2*i+1]=2.*eigxx[2*i]*eigxx[2*i+1];
1034           }
1035           
1036           if(fwrite(eigx,sizeof(double),2*nev,f1)!=2*nev){
1037               printf("*ERROR in complexfreq saving the eigenfrequencies to the eigenvalue file...");
1038               exit(0);
1039           }
1040           
1041           SFREE(eigx);
1042           
1043           /* storing the mass matrix */
1044           
1045           if(fwrite(adb,sizeof(double),*neq,f1)!=*neq){
1046               printf("*ERROR in complexfreq saving the diagonal of the mass matrix to the eigenvalue file...");
1047               exit(0);
1048           }
1049           if(fwrite(aub,sizeof(double),*nzs,f1)!=*nzs){
1050               printf("*ERROR in complexfreq saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
1051               exit(0);
1052           }
1053           
1054           /* storing the eigenvectors */
1055           
1056           lfin=0;
1057           for(j=0;j<nev;++j){
1058               lint=lfin;
1059               lfin=lfin+neq[1];
1060               if(fwrite(&zz[lint],sizeof(double),neq[1],f1)!=neq[1]){
1061                   printf("*ERROR in complexfreq saving the eigenvectors to the eigenvalue file...");
1062                   exit(0);
1063               }
1064           }
1065           
1066           /* storing the orthogonality matrices */
1067           
1068           if(fwrite(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
1069               printf("*ERROR in complexfreq saving the real orthogonality matrix to the eigenvalue file...");
1070               exit(0);
1071           }
1072           
1073           if(fwrite(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
1074               printf("*ERROR in complexfreq saving the imaginary orthogonality matrix to the eigenvalue file...");
1075               exit(0);
1076           }
1077       }
1078       SFREE(xmr);SFREE(xmi);
1079       
1080       fclose(f1);
1081   }
1082 
1083   SFREE(adb);SFREE(aub);
1084   if(!cyclicsymmetry){SFREE(ad);SFREE(au);}
1085 
1086   /* calculating the displacements and the stresses and storing */
1087   /* the results in frd format for each valid eigenmode */
1088 
1089   NNEW(v,double,2*mt**nk);
1090   NNEW(fn,double,2*mt**nk);
1091   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1653],"MAXS")==0)|| 
1092      (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1093      (strcmp1(&filab[1044],"ERR ")==0)) 
1094       NNEW(stn,double,12**nk);
1095 
1096   if((strcmp1(&filab[261],"E   ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) 
1097       NNEW(een,double,12**nk);
1098   if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,2**nk);
1099   if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,12**nk);
1100 
1101   NNEW(inum,ITG,*nk);
1102   NNEW(stx,double,2*6*mi[0]**ne);
1103   
1104   NNEW(coefmpcnew,double,*mpcend);
1105 
1106   NNEW(cot,double,3**nk*ngraph);
1107   if(*ntrans>0){NNEW(inotrt,ITG,2**nk*ngraph);}
1108   if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0))
1109 
1110 // real and imaginary part of the displacements
1111 
1112     NNEW(vt,double,2*mt**nk*ngraph);
1113   if(strcmp1(&filab[87],"NT  ")==0)
1114     NNEW(t1t,double,*nk*ngraph);
1115   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1116      (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1117 
1118 // real and imaginary part of the stresses
1119 
1120     NNEW(stnt,double,2*6**nk*ngraph);
1121   if(strcmp1(&filab[261],"E   ")==0) 
1122       NNEW(eent,double,2*6**nk*ngraph);
1123   if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0))
1124 
1125 // real and imaginary part of the forces
1126 
1127     NNEW(fnt,double,2*mt**nk*ngraph);
1128   if(strcmp1(&filab[522],"ENER")==0)
1129     NNEW(enernt,double,*nk*ngraph);
1130   if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1131     NNEW(stxt,double,2*6*mi[0]**ne*ngraph);
1132   if(strcmp1(&filab[2697],"ME  ")==0) 
1133       NNEW(emnt,double,2*6**nk*ngraph);
1134 
1135   NNEW(kont,ITG,*nkon*ngraph);
1136   NNEW(ipkont,ITG,*ne*ngraph);
1137   for(l=0;l<*ne*ngraph;l++)ipkont[l]=-1;
1138   NNEW(lakont,char,8**ne*ngraph);
1139   NNEW(inumt,ITG,*nk*ngraph);
1140   NNEW(ielmatt,ITG,mi[2]**ne*ngraph);
1141 
1142   nkt=ngraph**nk;
1143   net=ngraph**ne;
1144 
1145   /* copying the coordinates of the first sector */
1146 
1147   for(l=0;l<3**nk;l++){cot[l]=co[l];}
1148   if(*ntrans>0){for(l=0;l<*nk;l++){inotrt[2*l]=inotr[2*l];}}
1149   for(l=0;l<*nkon;l++){kont[l]=kon[l];}
1150   for(l=0;l<*ne;l++){ipkont[l]=ipkon[l];}
1151   for(l=0;l<8**ne;l++){lakont[l]=lakon[l];}
1152   for(l=0;l<*ne;l++){ielmatt[mi[2]*l]=ielmat[mi[2]*l];}
1153 
1154   /* generating the coordinates for the other sectors */
1155 
1156   if(cyclicsymmetry){
1157 
1158     icntrl=1;
1159     
1160     FORTRAN(rectcyl,(cot,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1161     
1162     for(jj=0;jj<*mcs;jj++){
1163       is=cs[17*jj+4];
1164       for(i=1;i<is;i++){
1165         
1166         theta=i*2.*pi/cs[17*jj];
1167         
1168         for(l=0;l<*nk;l++){
1169           if(inocs[l]==jj){
1170             cot[3*l+i*3**nk]=cot[3*l];
1171             cot[1+3*l+i*3**nk]=cot[1+3*l]+theta;
1172             cot[2+3*l+i*3**nk]=cot[2+3*l];
1173             if(*ntrans>0){inotrt[2*l+i*2**nk]=inotrt[2*l];}
1174           }
1175         }
1176         for(l=0;l<*nkon;l++){kont[l+i**nkon]=kon[l]+i**nk;}
1177         for(l=0;l<*ne;l++){
1178           if(ielcs[l]==jj){
1179             if(ipkon[l]>=0){
1180               ipkont[l+i**ne]=ipkon[l]+i**nkon;
1181               ielmatt[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
1182               for(l1=0;l1<8;l1++){
1183                 l2=8*l+l1;
1184                 lakont[l2+i*8**ne]=lakon[l2];
1185               }
1186             }
1187           }
1188         }
1189       }
1190     }
1191     
1192     icntrl=-1;
1193     
1194     FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1195                      &imag,mi,emnt));
1196   }
1197 
1198   /* check that the tensor fields which are extrapolated from the
1199      integration points are requested in global coordinates */
1200 
1201   if(strcmp1(&filab[174],"S   ")==0){
1202       if((strcmp1(&filab[179],"L")==0)&&(*norien>0)){
1203       printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculations\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1204       strcpy1(&filab[179],"G",1);
1205     }
1206   }
1207 
1208   if(strcmp1(&filab[261],"E   ")==0){
1209       if((strcmp1(&filab[266],"L")==0)&&(*norien>0)){
1210       printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1211       strcpy1(&filab[266],"G",1);
1212     }
1213   }
1214 
1215   if(strcmp1(&filab[1479],"PHS ")==0){
1216       if((strcmp1(&filab[1484],"L")==0)&&(*norien>0)){
1217       printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1218       strcpy1(&filab[1484],"G",1);
1219     }
1220   }
1221 
1222   if(strcmp1(&filab[1653],"MAXS")==0){
1223       if((strcmp1(&filab[1658],"L")==0)&&(*norien>0)){
1224       printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1225       strcpy1(&filab[1658],"G",1);
1226     }
1227   }   
1228 
1229   if(strcmp1(&filab[2523],"MAXE")==0){
1230       if((strcmp1(&filab[2528],"L")==0)&&(*norien>0)){
1231       printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1232       strcpy1(&filab[1658],"G",1);
1233     }
1234   }   
1235 
1236   /* allocating fields for magnitude and phase information of
1237      displacements and stresses */
1238 
1239   if(strcmp1(&filab[870],"PU")==0){
1240     NNEW(vr,double,mt*nkt);
1241     NNEW(vi,double,mt*nkt);
1242   }
1243 
1244   if(strcmp1(&filab[1479],"PHS")==0){
1245     NNEW(stnr,double,6*nkt);
1246     NNEW(stni,double,6*nkt);
1247   }
1248 
1249   if(strcmp1(&filab[1566],"MAXU")==0){
1250     NNEW(vmax,double,4*nkt);
1251   }
1252 
1253   if(strcmp1(&filab[1653],"MAXS")==0){
1254     NNEW(stnmax,double,7*nkt);
1255   }
1256 
1257   if(strcmp1(&filab[2523],"MAXE")==0){
1258     NNEW(eenmax,double,7*nkt);
1259   }
1260 
1261   /* storing the results */
1262 
1263   if(!cyclicsymmetry) (neq[1])*=2;
1264 
1265   lfin=0;
1266   for(j=0;j<nev;++j){
1267     lint=lfin;
1268     lfin=lfin+neq[1];
1269 
1270     /* calculating the cosine and sine */
1271 
1272     for(k=0;k<*mcs;k++){
1273         theta=nm[j]*2.*pi/cs[17*k];
1274         cs[17*k+14]=cos(theta);
1275         cs[17*k+15]=sin(theta);
1276     }
1277 
1278     if(*nprint>0)FORTRAN(writehe,(&j));
1279 
1280     NNEW(eei,double,6*mi[0]**ne);
1281     if(*nener==1){
1282         NNEW(stiini,double,6*mi[0]**ne);
1283         NNEW(enerini,double,mi[0]**ne);}
1284 
1285     DMEMSET(v,0,2*mt**nk,0.);
1286 
1287     for(k=0;k<neq[1];k+=neq[1]/2){
1288 
1289       for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
1290 
1291       if(k==0) {kk=0;kkv=0;kk6=0;kkx=0;if(*nprint>0)FORTRAN(writere,());}
1292       else {kk=*nk;kkv=mt**nk;kk6=6**nk;kkx=6*mi[0]**ne;
1293             if(*nprint>0)FORTRAN(writeim,());}
1294 
1295     /* generating the cyclic MPC's (needed for nodal diameters
1296        different from 0 */
1297 
1298       for(i=0;i<*nmpc;i++){
1299         index=ipompc[i]-1;
1300         /* check whether thermal mpc */
1301         if(nodempc[3*index+1]==0) continue;
1302         coefmpcnew[index]=coefmpc[index];
1303         while(1){
1304           index=nodempc[3*index+2];
1305           if(index==0) break;
1306           index--;
1307 
1308           icomplex=0;
1309           inode=nodempc[3*index];
1310           if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
1311             icomplex=atoi(&labmpc[20*i+6]);}
1312           else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
1313             for(ij=0;ij<*mcs;ij++){
1314               lprev=cs[ij*17+13];
1315               ilength=cs[ij*17+3];
1316               FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
1317               if(id!=0){
1318                 if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
1319               }
1320             }
1321           }
1322 
1323           if(icomplex!=0){
1324             idir=nodempc[3*index+1];
1325             idof=nactdof[mt*(inode-1)+idir]-1;
1326             if(idof==-1){xreal=1.;ximag=1.;}
1327             else{xreal=zz[lint+idof];ximag=zz[lint+idof+neq[1]/2];}
1328             if(k==0) {
1329               if(fabs(xreal)<1.e-30)xreal=1.e-30;
1330               coefmpcnew[index]=coefmpc[index]*
1331                 (cs[17*(icomplex-1)+14]+ximag/xreal*cs[17*(icomplex-1)+15]);}
1332             else {
1333               if(fabs(ximag)<1.e-30)ximag=1.e-30;
1334               coefmpcnew[index]=coefmpc[index]*
1335                 (cs[17*(icomplex-1)+14]-xreal/ximag*cs[17*(icomplex-1)+15]);}
1336           }
1337           else{coefmpcnew[index]=coefmpc[index];}
1338         }
1339       }
1340 
1341       if(*iperturb==0){
1342         results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1343             &stx[kkx],elcon,
1344             nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1345             norien,orab,ntmat_,t0,t0,ithermal,
1346             prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1347             f,&fn[kkv],nactdof,&iout,qa,vold,&zz[lint+k],
1348             nodeboun,ndirboun,xboun,nboun,ipompc,
1349             nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1350             &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1351             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1352             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1353             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1354             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1355             nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1356             &ne0,xforc,nforc,thicke,shcon,nshcon,
1357             sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1358             &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1359             islavsurf,ielprop,prop);}
1360       else{
1361         results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1362             &stx[kkx],elcon,
1363             nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1364             norien,orab,ntmat_,t0,t1old,ithermal,
1365             prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1366             f,&fn[kkv],nactdof,&iout,qa,vold,&zz[lint+k],
1367             nodeboun,ndirboun,xboun,nboun,ipompc,
1368             nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1369             &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1370             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1371             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1372             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1373             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1374             nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1375             &ne0,xforc,nforc,thicke,shcon,nshcon,
1376             sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1377             &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1378             islavsurf,ielprop,prop);
1379       }
1380 
1381     }
1382     SFREE(eei);
1383     if(*nener==1){SFREE(stiini);SFREE(enerini);}
1384    
1385     /* changing the basic results into cylindrical coordinates
1386        (only for cyclic symmetry structures */
1387 
1388     for(l=0;l<*nk;l++){inumt[l]=inum[l];}
1389 
1390     if(cyclicsymmetry){
1391       icntrl=2;imag=1;      
1392       FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1393     }
1394 
1395     /* copying the basis results (real and imaginary) into
1396        a new field */
1397 
1398     if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1399       for(l=0;l<mt**nk;l++){vt[l]=v[l];}
1400       for(l=0;l<mt**nk;l++){vt[l+mt**nk*ngraph]=v[l+mt**nk];}}
1401     if(strcmp1(&filab[87],"NT  ")==0)
1402       for(l=0;l<*nk;l++){t1t[l]=t1[l];};
1403     if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1404         for(l=0;l<6**nk;l++){stnt[l]=stn[l];}
1405         for(l=0;l<6**nk;l++){stnt[l+6**nk*ngraph]=stn[l+6**nk];}}
1406     if(strcmp1(&filab[261],"E   ")==0){
1407         for(l=0;l<6**nk;l++){eent[l]=een[l];};
1408         for(l=0;l<6**nk;l++){eent[l+6**nk*ngraph]=een[l+6**nk];}}
1409     if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1410       for(l=0;l<mt**nk;l++){fnt[l]=fn[l];}
1411       for(l=0;l<mt**nk;l++){fnt[l+mt**nk*ngraph]=fn[l+mt**nk];}}
1412     if(strcmp1(&filab[522],"ENER")==0)
1413       for(l=0;l<*nk;l++){enernt[l]=enern[l];};
1414     if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)){
1415       for(l=0;l<6*mi[0]**ne;l++){stxt[l]=stx[l];}
1416       for(l=0;l<6*mi[0]**ne;l++){stxt[l+6*mi[0]**ne*ngraph]=stx[l+6*mi[0]**ne];}}
1417     if(strcmp1(&filab[2697],"ME  ")==0){
1418         for(l=0;l<6**nk;l++){emnt[l]=emn[l];};
1419         for(l=0;l<6**nk;l++){emnt[l+6**nk*ngraph]=emn[l+6**nk];}}
1420     
1421     /* mapping the results to the other sectors
1422        (only for cyclic symmetric structures */
1423 
1424     if(cyclicsymmetry){
1425 
1426       for(jj=0;jj<*mcs;jj++){
1427         ilength=cs[17*jj+3];
1428         is=cs[17*jj+4];
1429         lprev=cs[17*jj+13];
1430         for(i=1;i<is;i++){
1431           
1432           for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
1433           
1434           theta=i*nm[j]*2.*pi/cs[17*jj];
1435           ctl=cos(theta);
1436           stl=sin(theta);
1437           
1438           if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1439             for(l1=0;l1<*nk;l1++){
1440               if(inocs[l1]==jj){
1441                 
1442                 /* check whether node lies on axis */
1443                 
1444                 ml1=-l1-1;
1445                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1446                 if(id!=0){
1447                   if(ics[lprev+id-1]==ml1){
1448                     for(l2=0;l2<4;l2++){
1449                       l=mt*l1+l2;
1450                       vt[l+mt**nk*i]=v[l];
1451                     }
1452                     continue;
1453                   }
1454                 }
1455                 for(l2=0;l2<4;l2++){
1456                   l=mt*l1+l2;
1457                   vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
1458                 }
1459                 
1460               }
1461             }
1462           }
1463           
1464           /* imaginary part of the displacements in cylindrical
1465              coordinates */
1466           
1467           if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1468             for(l1=0;l1<*nk;l1++){
1469               if(inocs[l1]==jj){
1470                 
1471                 /* check whether node lies on axis */
1472                 
1473                 ml1=-l1-1;
1474                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1475                 if(id!=0){
1476                   if(ics[lprev+id-1]==ml1){
1477                     for(l2=0;l2<4;l2++){
1478                       l=mt*l1+l2;
1479                       vt[l+mt**nk*(i+ngraph)]=v[l+mt**nk];
1480                     }
1481                     continue;
1482                   }
1483                 }
1484                 for(l2=0;l2<4;l2++){
1485                   l=mt*l1+l2;
1486                   vt[l+mt**nk*(i+ngraph)]=stl*v[l]+ctl*v[l+mt**nk];
1487                 }
1488               }
1489             }
1490           }
1491           
1492           if(strcmp1(&filab[87],"NT  ")==0){
1493             for(l=0;l<*nk;l++){
1494               if(inocs[l]==jj) t1t[l+*nk*i]=t1[l];
1495             }
1496           }
1497           
1498           if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1499             for(l1=0;l1<*nk;l1++){
1500               if(inocs[l1]==jj){
1501                 
1502                 /* check whether node lies on axis */
1503                 
1504                 ml1=-l1-1;
1505                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1506                 if(id!=0){
1507                   if(ics[lprev+id-1]==ml1){
1508                     for(l2=0;l2<6;l2++){
1509                       l=6*l1+l2;
1510                       stnt[l+6**nk*i]=stn[l];
1511                     }
1512                     continue;
1513                   }
1514                 }
1515                 for(l2=0;l2<6;l2++){
1516                   l=6*l1+l2;
1517                   stnt[l+6**nk*i]=ctl*stn[l]-stl*stn[l+6**nk];
1518                 }
1519               }
1520             }
1521           }
1522           
1523           /* imaginary part of the stresses in cylindrical
1524              coordinates */
1525           
1526           if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1527             for(l1=0;l1<*nk;l1++){
1528               if(inocs[l1]==jj){
1529                 
1530                 /* check whether node lies on axis */
1531                 
1532                 ml1=-l1-1;
1533                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1534                 if(id!=0){
1535                   if(ics[lprev+id-1]==ml1){
1536                     for(l2=0;l2<6;l2++){
1537                       l=6*l1+l2;
1538                       stnt[l+6**nk*(i+ngraph)]=stn[l+6**nk];
1539                     }
1540                     continue;
1541                   }
1542                 }
1543                 for(l2=0;l2<6;l2++){
1544                   l=6*l1+l2;
1545                   stnt[l+6**nk*(i+ngraph)]=stl*stn[l]+ctl*stn[l+6**nk];
1546                 }
1547               }
1548             }
1549           }
1550           
1551           if(strcmp1(&filab[261],"E   ")==0){
1552             for(l1=0;l1<*nk;l1++){
1553               if(inocs[l1]==jj){
1554                 
1555                 /* check whether node lies on axis */
1556                 
1557                 ml1=-l1-1;
1558                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1559                 if(id!=0){
1560                   if(ics[lprev+id-1]==ml1){
1561                     for(l2=0;l2<6;l2++){
1562                       l=6*l1+l2;
1563                       eent[l+6**nk*i]=een[l];
1564                     }
1565                     continue;
1566                   }
1567                 }
1568                 for(l2=0;l2<6;l2++){
1569                   l=6*l1+l2;
1570                   eent[l+6**nk*i]=ctl*een[l]-stl*een[l+6**nk];
1571                 }
1572               }
1573             }
1574           }
1575         
1576         /* imaginary part of the strains in cylindrical
1577            coordinates */
1578         
1579           if(strcmp1(&filab[261],"E   ")==0){
1580               for(l1=0;l1<*nk;l1++){
1581                   if(inocs[l1]==jj){
1582                       
1583                       /* check whether node lies on axis */
1584                       
1585                       ml1=-l1-1;
1586                       FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1587                       if(id!=0){
1588                           if(ics[lprev+id-1]==ml1){
1589                               for(l2=0;l2<6;l2++){
1590                                   l=6*l1+l2;
1591                                   eent[l+6**nk*(i+ngraph)]=een[l+6**nk];
1592                               }
1593                               continue;
1594                           }
1595                       }
1596                       for(l2=0;l2<6;l2++){
1597                           l=6*l1+l2;
1598                           eent[l+6**nk*(i+ngraph)]=stl*een[l]+ctl*een[l+6**nk];
1599                       }
1600                   }
1601               }
1602           }
1603           
1604           if(strcmp1(&filab[2697],"ME  ")==0){
1605             for(l1=0;l1<*nk;l1++){
1606               if(inocs[l1]==jj){
1607                 
1608                 /* check whether node lies on axis */
1609                 
1610                 ml1=-l1-1;
1611                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1612                 if(id!=0){
1613                   if(ics[lprev+id-1]==ml1){
1614                     for(l2=0;l2<6;l2++){
1615                       l=6*l1+l2;
1616                       emnt[l+6**nk*i]=emn[l];
1617                     }
1618                     continue;
1619                   }
1620                 }
1621                 for(l2=0;l2<6;l2++){
1622                   l=6*l1+l2;
1623                   emnt[l+6**nk*i]=ctl*emn[l]-stl*emn[l+6**nk];
1624                 }
1625               }
1626             }
1627           }
1628         
1629         /* imaginary part of the mechanical strains in cylindrical
1630            coordinates */
1631         
1632           if(strcmp1(&filab[2697],"ME  ")==0){
1633               for(l1=0;l1<*nk;l1++){
1634                   if(inocs[l1]==jj){
1635                       
1636                       /* check whether node lies on axis */
1637                       
1638                       ml1=-l1-1;
1639                       FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1640                       if(id!=0){
1641                           if(ics[lprev+id-1]==ml1){
1642                               for(l2=0;l2<6;l2++){
1643                                   l=6*l1+l2;
1644                                   emnt[l+6**nk*(i+ngraph)]=emn[l+6**nk];
1645                               }
1646                               continue;
1647                           }
1648                       }
1649                       for(l2=0;l2<6;l2++){
1650                           l=6*l1+l2;
1651                           emnt[l+6**nk*(i+ngraph)]=stl*emn[l]+ctl*emn[l+6**nk];
1652                       }
1653                   }
1654               }
1655           }
1656           
1657           if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1658             for(l1=0;l1<*nk;l1++){
1659               if(inocs[l1]==jj){
1660                 
1661                 /* check whether node lies on axis */
1662                 
1663                 ml1=-l1-1;
1664                 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1665                 if(id!=0){
1666                   if(ics[lprev+id-1]==ml1){
1667                     for(l2=0;l2<4;l2++){
1668                       l=mt*l1+l2;
1669                       fnt[l+mt**nk*i]=fn[l];
1670                     }
1671                     continue;
1672                   }
1673                 }
1674                 for(l2=0;l2<4;l2++){
1675                   l=mt*l1+l2;
1676                   fnt[l+mt**nk*i]=ctl*fn[l]-stl*fn[l+mt**nk];
1677                 }
1678               }
1679             }
1680           }
1681         
1682         /* imaginary part of the forces in cylindrical
1683            coordinates */
1684 
1685           if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1686               for(l1=0;l1<*nk;l1++){
1687                   if(inocs[l1]==jj){
1688                       
1689                       /* check whether node lies on axis */
1690                       
1691                       ml1=-l1-1;
1692                       FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1693                       if(id!=0){
1694                           if(ics[lprev+id-1]==ml1){
1695                               for(l2=0;l2<4;l2++){
1696                                   l=mt*l1+l2;
1697                                   fnt[l+mt**nk*(i+ngraph)]=fn[l+mt**nk];
1698                               }
1699                               continue;
1700                           }
1701                       }
1702                       for(l2=0;l2<4;l2++){
1703                           l=mt*l1+l2;
1704                           fnt[l+mt**nk*(i+ngraph)]=stl*fn[l]+ctl*fn[l+mt**nk];
1705                       }
1706                   }
1707               }
1708           }
1709           
1710           if(strcmp1(&filab[522],"ENER")==0){
1711             for(l=0;l<*nk;l++){
1712               if(inocs[l]==jj) enernt[l+*nk*i]=0.;
1713             }
1714           }
1715         }
1716       }
1717       
1718       icntrl=-2;imag=0;
1719       
1720       FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1721                        &imag,mi,emnt));
1722       
1723       FORTRAN(rectcylvi,(cot,&vt[mt**nk*ngraph],&fnt[mt**nk*ngraph],
1724                          &stnt[6**nk*ngraph],qfnt,&eent[6**nk*ngraph],
1725                          cs,&nkt,&icntrl,t,filab,&imag,mi,&emnt[6**nk*ngraph]));
1726       
1727     }
1728 
1729     /* determining magnitude and phase angle for the displacements */
1730 
1731     if(strcmp1(&filab[870],"PU")==0){
1732       for(l1=0;l1<nkt;l1++){
1733         for(l2=0;l2<4;l2++){
1734           l=mt*l1+l2;
1735           vreal=vt[l];
1736           vimag=vt[l+mt**nk*ngraph];
1737           vr[l]=sqrt(vreal*vreal+vimag*vimag);
1738           if(fabs(vreal)<1.e-10){
1739             if(vimag>0){vi[l]=90.;}
1740             else{vi[l]=-90.;}
1741           }
1742           else{
1743             vi[l]=atan(vimag/vreal)*constant;
1744             if(vreal<0) vi[l]+=180.;
1745           }
1746         }
1747       }
1748     }
1749 
1750     /* determining magnitude and phase for the stress */
1751 
1752     if(strcmp1(&filab[1479],"PHS")==0){
1753       for(l1=0;l1<nkt;l1++){
1754         for(l2=0;l2<6;l2++){
1755           l=6*l1+l2;
1756           stnreal=stnt[l];
1757           stnimag=stnt[l+6**nk*ngraph];
1758           stnr[l]=sqrt(stnreal*stnreal+stnimag*stnimag);
1759           if(fabs(stnreal)<1.e-10){
1760             if(stnimag>0){stni[l]=90.;}
1761             else{stni[l]=-90.;}
1762           }
1763           else{
1764             stni[l]=atan(stnimag/stnreal)*constant;
1765             if(stnreal<0) stni[l]+=180.;
1766           }
1767         }
1768       }
1769     }
1770 
1771     ++*kode;
1772 
1773     /* storing the real part of the eigenfrequencies in freq */
1774 
1775     freq=eigxx[2*j]/6.283185308;
1776     if(strcmp1(&filab[1044],"ZZS")==0){
1777         NNEW(neigh,ITG,40*net);
1778         NNEW(ipneigh,ITG,nkt);
1779     }
1780     frd(cot,&nkt,kont,ipkont,lakont,&net,vt,stnt,inumt,nmethod,
1781             kode,filab,eent,t1t,fnt,&freq,epn,ielmatt,matname,enernt,xstaten,
1782             nstate_,istep,&iinc,ithermal,qfn,&j,&nm[j],trab,inotrt,
1783             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1784             mi,stxt,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&net,
1785             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emnt,
1786             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1787     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1788   }
1789 
1790   SFREE(xstiff);if(*nbody>0) SFREE(ipobody);
1791   SFREE(cstr);SFREE(zz);SFREE(eigxx);SFREE(xx);
1792 
1793   if(cyclicsymmetry){
1794       SFREE(istartnmd);SFREE(iendnmd);
1795   }else{
1796       (neq[1])/=2;
1797   }
1798 
1799   SFREE(nm);SFREE(coefmpcnew);
1800 
1801   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1653],"MAXS")==0)|| 
1802      (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1803      (strcmp1(&filab[1044],"ERR ")==0)) 
1804      SFREE(stn);
1805 
1806   SFREE(v);SFREE(fn);SFREE(inum);SFREE(stx);
1807 
1808   if((strcmp1(&filab[261],"E   ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) SFREE(een);
1809   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1810   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
1811 
1812   if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)) SFREE(vt);
1813   if(strcmp1(&filab[87],"NT  ")==0) SFREE(t1t);
1814   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1815      (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stnt);
1816   if(strcmp1(&filab[261],"E   ")==0) SFREE(eent);
1817   if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)) SFREE(fnt);
1818   if(strcmp1(&filab[522],"ENER")==0) SFREE(enernt);
1819   if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stxt);
1820   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emnt);
1821 
1822   SFREE(cot);SFREE(kont);SFREE(ipkont);SFREE(lakont);SFREE(inumt);SFREE(ielmatt);
1823   if(*ntrans>0){SFREE(inotrt);}
1824 
1825   *ialsetp=ialset;
1826   *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
1827   *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
1828   *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;
1829   *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
1830   *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
1831   *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
1832   *fmpcp=fmpc;*veoldp=veold;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
1833   *stip=sti;
1834 
1835   return;
1836 }
1837 

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