root/src/steadystate.c

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

DEFINITIONS

This source file includes following definitions.
  1. steadystate

   1 /*     CalculiX - A 3-dimensional finite element program                   */
   2 /*              Copyright (C) 1998-2015 Guido Dhondt                          */
   3 
   4 /*     This program is free software; you can redistribute it and/or     */
   5 /*     modify it under the terms of the GNU General Public License as    */
   6 /*     published by the Free Software Foundation(version 2);    */
   7 /*                    */
   8 
   9 /*     This program is distributed in the hope that it will be useful,   */
  10 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  11 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  12 /*     GNU General Public License for more details.                      */
  13 
  14 /*     You should have received a copy of the GNU General Public License */
  15 /*     along with this program; if not, write to the Free Software       */
  16 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  17 
  18 #include <stdio.h>
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include "CalculiX.h"
  23 
  24 #ifdef SPOOLES
  25    #include "spooles.h"
  26 #endif
  27 #ifdef SGI
  28    #include "sgi.h"
  29 #endif
  30 #ifdef TAUCS
  31    #include "tau.h"
  32 #endif
  33 #ifdef PARDISO
  34    #include "pardiso.h"
  35 #endif
  36 
  37 void steadystate(double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, 
  38                ITG **nodebounp, ITG **ndirbounp, double **xbounp, ITG *nboun,
  39                ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp,
  40                ITG *nmpc, ITG *nodeforc,ITG *ndirforc,double *xforc, 
  41                ITG *nforc,ITG *nelemload, char *sideload,double *xload,
  42                ITG *nload, 
  43                ITG **nactdofp,ITG *neq, ITG *nzl,ITG *icol, ITG *irow, 
  44                ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG **ikbounp, 
  45                ITG **ilbounp,double *elcon, ITG *nelcon, double *rhcon, 
  46                ITG *nrhcon,double *cocon, ITG *ncocon,
  47                double *alcon, ITG *nalcon, double *alzero, 
  48                ITG **ielmatp,ITG **ielorienp, ITG *norien, double *orab, 
  49                ITG *ntmat_,double **t0p, 
  50                double **t1p,ITG *ithermal,double *prestr, ITG *iprestr, 
  51                double **voldp,ITG *iperturb, double *sti, ITG *nzs, 
  52                double *timepar, double *xmodal,
  53                double **veoldp, char *amname, double *amta,
  54                ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload,
  55                ITG **iamt1p,ITG *jout,
  56                ITG *kode, char *filab,double **emep, double *xforcold, 
  57                double *xloadold,
  58                double **t1oldp, ITG **iambounp, double **xbounoldp, ITG *iexpl,
  59                double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
  60                double *xstate, ITG *npmat_, char *matname, ITG *mi,
  61                ITG *ncmat_, ITG *nstate_, double **enerp, char *jobnamec,
  62                double *ttime, char *set, ITG *nset, ITG *istartset,
  63                ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  64                char *prset, ITG *nener, double *trab, 
  65                ITG **inotrp, ITG *ntrans, double **fmpcp, char *cbody, ITG *ibody,
  66                double *xbody, ITG *nbody, double *xbodyold, ITG *istep,
  67                ITG *isolver, ITG *jq, char *output, ITG *mcs,ITG *nkon, 
  68                ITG *ics, double *cs, ITG *mpcend,double *ctrl,
  69                ITG *ikforc, ITG *ilforc, double *thicke,ITG *nmat,
  70                char *typeboun,ITG *ielprop,double *prop){
  71 
  72   char fneig[132]="",description[13]="            ",*lakon=NULL,*labmpc=NULL,
  73     *labmpcold=NULL,cflag[1]=" ";
  74 
  75   ITG nev,i,j,k, *inum=NULL,*ipobody=NULL,inewton=0,nsectors,im,
  76     iinc=0,l,iout,ielas,icmd,iprescribedboundary,ndata,nmd,nevd,
  77     ndatatot,*iphaseforc=NULL,*iphaseload=NULL,*iphaseboun=NULL,
  78     *isave=NULL,nfour,ii,ir,ic,mode,noddiam=-1,*nm=NULL,*islavact=NULL,
  79     *kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,*inotr=NULL,
  80     *nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,jj,
  81     *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
  82     *ilmpc=NULL,*ipompcold=NULL,*nodempcold=NULL,*ikmpcold=NULL,
  83     *ilmpcold=NULL,nmpcold,mpcendold,kflag=2,*iamt1=NULL,ifreebody,
  84     *itg=NULL,ntg=0,symmetryflag=0,inputformat=0,dashpot,nrhs=1,
  85     *ipiv=NULL,info,nev2,ngraph=1,nkg,neg,iflag=1,idummy=1,imax,
  86     nzse[3],mt=mi[1]+1,*ikactmech=NULL,nactmech,id,nasym=0,
  87     *imddof=NULL,nmddof,*imdnode=NULL,nmdnode,*imdboun=NULL,nmdboun,
  88     *imdmpc=NULL,nmdmpc,*izdof=NULL,nzdof,cyclicsymmetry,mortar=0,
  89     *ikactmechr=NULL,*ikactmechi=NULL,nactmechr,nactmechi,intpointvar,
  90     iforc,iload,ne0,*iponoel=NULL,*inoel=NULL,*imdelem=NULL,
  91     nmdelem,*integerglob=NULL,*nshcon=NULL,nherm,icfd=0,*inomat=NULL,
  92     *islavnode=NULL,*nslavnode=NULL,*islavsurf=NULL;
  93 
  94   long long i2;
  95 
  96   double *d=NULL, *z=NULL,*stiini=NULL,*vini=NULL,*freqnh=NULL,
  97     *xforcact=NULL, *xloadact=NULL,y,*fr=NULL,*fi=NULL,*cc=NULL,
  98     *t1act=NULL,*ampli=NULL, *aa=NULL,*bb=NULL,*vr=NULL,*vi=NULL,
  99     *stn=NULL,*stx=NULL,*een=NULL,*adb=NULL,*xstiff=NULL,ptime,
 100     *aub=NULL,*bjr=NULL, *bji=NULL,*xbodyr=NULL,*cdn=NULL,
 101     *f=NULL, *fn=NULL, *xbounact=NULL,*epn=NULL,*xstateini=NULL,
 102     *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
 103     *qfx=NULL, *xbodyact=NULL, *cgr=NULL, *au=NULL,*xbodyi=NULL,
 104     time,dtime,reltime,*co=NULL,*xboun=NULL,*xbounold=NULL,
 105     physcon[1],qa[3],cam[5],accold[1],bet,gam,*emn=NULL,timem,
 106     *ad=NULL,sigma=0.,alpham,betam,*fnr=NULL,*fni=NULL,*emeini=NULL,
 107     fmin,fmax,bias,*freq=NULL,*xforcr=NULL,dd,pi,vreal,constant,
 108     *xforci=NULL,*xloadr=NULL,*xloadi=NULL,*xbounr=NULL,*xbouni=NULL,
 109     *br=NULL,*bi=NULL,*ubr=NULL,*ubi=NULL,*mubr=NULL,*mubi=NULL,
 110     *wsave=NULL,*r=NULL,*xbounacttime=NULL,*btot=NULL,breal,tmin,tmax,
 111     *vold=NULL,*eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,
 112     *coefmpcold=NULL,*t0=NULL,*t1=NULL,*t1old=NULL,*adc=NULL,*auc=NULL,
 113     *am=NULL,*bm=NULL,*zc=NULL,*e=NULL,*stnr=NULL,*stni=NULL,
 114     *vmax=NULL,*stnmax=NULL,*va=NULL,*vp=NULL,*fric=NULL,*springarea=NULL,
 115     *stna=NULL,*stnp=NULL,*bp=NULL,*eenmax=NULL,*clearini=NULL,
 116     *doubleglob=NULL,*shcon=NULL,*veold=NULL,*xmr=NULL,*xmi=NULL,*eig=NULL,
 117     *ax=NULL,*bx=NULL,*pslavsurf=NULL,*pmastsurf=NULL,xnull=0.,
 118     *cdnr=NULL,*cdni=NULL,*tinc,*tper;
 119 
 120   /* dummy arguments for the call of expand*/
 121 
 122   char* tieset=NULL;
 123   ITG *jqe=NULL,*icole=NULL,*irowe=NULL,ntie=0;
 124   double *adbe=NULL,*aube=NULL;
 125 
 126   FILE *f1;
 127 
 128   ITG *ipneigh=NULL,*neigh=NULL;
 129 
 130 #ifdef SGI
 131   ITG token;
 132 #endif
 133 
 134   co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
 135   ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
 136   ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;veold=*veoldp;
 137   xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
 138   vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
 139   coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 140   fmpc=*fmpcp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
 141 
 142   tinc=&timepar[0];
 143   tper=&timepar[1];
 144 
 145   pi=4.*atan(1.);
 146   iout=2;
 147 
 148   alpham=xmodal[0];
 149   betam=xmodal[1];
 150 
 151   fmin=2.*pi*xmodal[2];
 152   fmax=2.*pi*xmodal[3];
 153   ndata=floor(xmodal[4]);
 154   bias=xmodal[5];
 155   nfour=floor(xmodal[6]);
 156   if(nfour>0){
 157       tmin=xmodal[7];
 158       tmax=xmodal[8];
 159   }
 160 
 161   /* determining nzl */
 162 
 163   *nzl=0;
 164   for(i=neq[1];i>0;i--){
 165       if(icol[i-1]>0){
 166           *nzl=i;
 167           break;
 168       }
 169   }
 170 
 171   /* opening the eigenvalue file and checking for cyclic symmetry */
 172 
 173   strcpy(fneig,jobnamec);
 174   strcat(fneig,".eig");
 175 
 176   if((f1=fopen(fneig,"rb"))==NULL){
 177     printf("*ERROR in steadystate: cannot open eigenvalue file for reading");
 178     exit(0);
 179   }
 180 
 181   printf(" *INFO  in steadystate: if there are problems reading the .eig file this may be due to:\n");
 182   printf("        1) the nonexistence of the .eig file\n");
 183   printf("        2) other boundary conditions than in the input deck\n");
 184   printf("           which created the .eig file\n\n");
 185 
 186   if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
 187       printf("*ERROR in steadystate reading the cyclic symmetry flag in the eigenvalue file");
 188       exit(0);
 189   }
 190 
 191   if(fread(&nherm,sizeof(ITG),1,f1)!=1){
 192       printf("*ERROR in steadystate reading the Hermitian flag in the eigenvalue file");
 193       exit(0);
 194   }
 195 
 196   /* creating imddof containing the degrees of freedom
 197      retained by the user and imdnode containing the nodes */
 198 
 199   nmddof=0;nmdnode=0;nmdboun=0;nmdmpc=0;nmdelem=0;
 200 
 201   NNEW(imddof,ITG,*nk*3);
 202   NNEW(imdnode,ITG,*nk);
 203   NNEW(imdboun,ITG,*nboun);
 204   NNEW(imdmpc,ITG,*nmpc);
 205   FORTRAN(createmddof,(imddof,&nmddof,istartset,iendset,
 206                        ialset,nactdof,ithermal,mi,imdnode,&nmdnode,
 207                        ikmpc,ilmpc,ipompc,nodempc,nmpc,
 208                        imdmpc,&nmdmpc,imdboun,&nmdboun,ikboun,
 209                        nboun,nset,&ntie,tieset,set,lakon,kon,ipkon,labmpc,
 210                        ilboun,filab,prlab,prset,nprint,ne,&cyclicsymmetry));
 211 
 212   /* if results are requested in too many nodes, it is faster to 
 213      calculate the results in all nodes */
 214 
 215   if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
 216       nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 217   }
 218   
 219   if(nmdnode!=0){
 220       if(!cyclicsymmetry){
 221           for(i=0;i<*nload;i++){
 222               iload=i+1;
 223               FORTRAN(addimdnodedload,(nelemload,sideload,ipkon,kon,lakon,
 224                       &iload,imdnode,&nmdnode,ikmpc,ilmpc,ipompc,nodempc,nmpc,
 225                       imddof,&nmddof,nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 226                       ikboun,nboun,ilboun,ithermal));
 227           }
 228           
 229           for(i=0;i<*nforc;i++){
 230               iforc=i+1;
 231               FORTRAN(addimdnodecload,(nodeforc,&iforc,imdnode,&nmdnode,xforc,
 232                    ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
 233                    nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 234                    ikboun,nboun,ilboun,ithermal));
 235           }
 236       }
 237           
 238       /* determining the elements belonging to a given node */
 239       
 240       NNEW(iponoel,ITG,*nk);
 241       NNEW(inoel,ITG,2**nkon);
 242       FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
 243       NNEW(imdelem,ITG,*ne);
 244 
 245       /* storing the elements in which integration point results
 246          are needed; storing the nodes which are needed to
 247          calculate these results */
 248 
 249       FORTRAN(createmdelem,(imdnode,&nmdnode,xforc,
 250                    ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
 251                    nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 252                    ikboun,nboun,ilboun,ithermal,imdelem,&nmdelem,
 253                    iponoel,inoel,prlab,prset,nprint,lakon,set,nset,
 254                    ialset,ipkon,kon,istartset,iendset,nforc,
 255                    ikforc,ilforc));
 256 
 257       RENEW(imdelem,ITG,nmdelem);
 258       SFREE(iponoel);SFREE(inoel);
 259   }
 260 
 261   /* if results are requested in too many nodes, it is faster to 
 262      calculate the results in all nodes */
 263 
 264   if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
 265       nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 266   }
 267   
 268   /* subtracting 1 to comply with the C-convention */
 269 
 270   for(i=0;i<nmddof;i++){imddof[i]-=1;}
 271 
 272   RENEW(imddof,ITG,nmddof);
 273   RENEW(imdnode,ITG,nmdnode);
 274   RENEW(imdboun,ITG,nmdboun);
 275   RENEW(imdmpc,ITG,nmdmpc);
 276 
 277   nsectors=1;
 278 
 279   if(!cyclicsymmetry){
 280 
 281       nkg=*nk;
 282       neg=*ne;
 283 
 284       if(fread(&nev,sizeof(ITG),1,f1)!=1){
 285           printf("*ERROR in steadystate reading the number of eigenvalues in the eigenvalue file");
 286           exit(0);
 287       }
 288 
 289       if(nherm==1){
 290           NNEW(d,double,nev);
 291           if(fread(d,sizeof(double),nev,f1)!=nev){
 292               printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
 293               exit(0);
 294           }
 295 
 296           for(i=0;i<nev;i++){
 297               if(d[i]<0){d[i]=0.;}
 298           }
 299       }else{
 300           NNEW(d,double,2*nev);
 301           if(fread(d,sizeof(double),2*nev,f1)!=2*nev){
 302               printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
 303               exit(0);
 304           }
 305       }
 306       
 307       NNEW(ad,double,neq[1]);
 308       NNEW(adb,double,neq[1]);
 309       NNEW(au,double,nzs[2]);
 310       NNEW(aub,double,nzs[1]);
 311       
 312       /* reading the stiffness matrix */
 313       
 314       if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
 315           printf("*ERROR in steadystate reading the diagonal of the stiffness matrix in the eigenvalue file");
 316           exit(0);
 317       }
 318       
 319       if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
 320           printf("*ERROR in steadystate reading the off-diagonals of the stiffness matrix in the eigenvalue file");
 321           exit(0);
 322       }
 323       
 324       /* reading the mass matrix */
 325       
 326       if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 327           printf("*ERROR in steadystate reading the diagonal of the mass matrix in the eigenvalue file");
 328           exit(0);
 329       }
 330       
 331       if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 332           printf("*ERROR in steadystate reading the off-diagonals of the mass matrix in the eigenvalue file");
 333           exit(0);
 334       }
 335       
 336       /* reading the eigenvectors */
 337 
 338       if(nherm==1){
 339           NNEW(z,double,neq[1]*nev);
 340           if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
 341               printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
 342               exit(0);
 343           }
 344       }else{
 345           NNEW(z,double,2*neq[1]*nev);
 346           if(fread(z,sizeof(double),2*neq[1]*nev,f1)!=2*neq[1]*nev){
 347               printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
 348               exit(0);
 349           }
 350       }
 351 
 352       /* reading the orthogonality matrices */
 353 
 354       if(nherm!=1){
 355           NNEW(xmr,double,nev*nev);
 356           NNEW(xmi,double,nev*nev);
 357           if(fread(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
 358               printf("*ERROR in steadystate reading the real orthogonality matrix to the eigenvalue file...");
 359               exit(0);
 360           }
 361           
 362           if(fread(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
 363               printf("*ERROR in steadystate reading the imaginary orthogonality matrix to the eigenvalue file...");
 364               exit(0);
 365           }
 366       }
 367   }
 368   else{
 369       nev=0;
 370       do{
 371           if(fread(&nmd,sizeof(ITG),1,f1)!=1){
 372               break;
 373           }
 374 
 375           if(fread(&nevd,sizeof(ITG),1,f1)!=1){
 376               printf("*ERROR in steadystate reading the number of eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 377               exit(0);
 378           }
 379 
 380           /* reading the eigenvalues (complex for non-Hermitian systems) */
 381 
 382           if(nev==0){
 383               if(nherm==1){NNEW(d,double,nevd);
 384               }else{NNEW(d,double,2*nevd);}
 385               NNEW(nm,ITG,nevd);
 386           }else{
 387               if(nherm!=1){
 388                   printf("*ERROR in steadystate: non-Hermitian systems cannot\n");
 389                   printf("       be combined with multiple modal diameters\n");
 390                   printf("       in cyclic symmetry calculations\n\n");
 391                   FORTRAN(stop,());
 392               }
 393               RENEW(d,double,nev+nevd);
 394               RENEW(nm,ITG,nev+nevd);
 395           }
 396           
 397           if(nherm==1){
 398               if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
 399                   printf("*ERROR in steadystate reading the eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 400                   exit(0);
 401               }
 402               for(i=nev;i<nev+nevd;i++){
 403                   if(d[i]<0){d[i]=0.;}
 404               }
 405           }else{
 406               if(fread(&d[nev],sizeof(double),2*nevd,f1)!=2*nevd){
 407                   printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
 408                   exit(0);
 409               }
 410           }
 411 
 412           for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
 413           
 414           if(nev==0){
 415               NNEW(adb,double,neq[1]);
 416               NNEW(aub,double,nzs[1]);
 417 
 418               if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 419                   printf("*ERROR in steadystate reading the diagonal of the mass matrix in the eigenvalue file");
 420                   exit(0);
 421               }
 422               
 423               if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 424                   printf("*ERROR in steadystate reading the off-diagonals of the mass matrix in the eigenvalue file");
 425                   exit(0);
 426               }
 427           }
 428           
 429           /* reading the eigenvectors */
 430 
 431           if(nev==0){
 432               NNEW(z,double,neq[1]*nevd);
 433           }else{
 434               RENEW(z,double,(long long)neq[1]*(nev+nevd));
 435           }
 436           
 437           if(fread(&z[(long long)neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
 438               printf("*ERROR in steadystate reading the eigenvectors for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 439               exit(0);
 440           }
 441 
 442           /* reading the orthogonality matrices */
 443           
 444           if(nherm!=1){
 445               NNEW(xmr,double,nev*nev);
 446               NNEW(xmi,double,nev*nev);
 447               if(fread(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
 448                   printf("*ERROR in steadystate reading the real orthogonality matrix to the eigenvalue file...");
 449                   exit(0);
 450               }
 451               
 452               if(fread(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
 453                   printf("*ERROR in steadystate reading the imaginary orthogonality matrix to the eigenvalue file...");
 454                   exit(0);
 455               }
 456           }
 457 
 458           nev+=nevd;
 459       }while(1);
 460 
 461       /* determining the maximum amount of sectors */
 462 
 463       for(i=0;i<*mcs;i++){
 464           if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
 465       }
 466 
 467         /* determining the maximum number of sectors to be plotted */
 468 
 469       for(j=0;j<*mcs;j++){
 470           if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
 471       }
 472       nkg=*nk*ngraph;
 473       neg=*ne*ngraph;
 474 
 475       /* allocating field for the expanded structure */
 476 
 477       RENEW(co,double,3**nk*nsectors);
 478 
 479       /* next line is necessary for multiple cyclic symmetry
 480          conditions */
 481 
 482       for(i=3**nk;i<3**nk*nsectors;i++){co[i]=0.;}
 483       if(*ithermal!=0){
 484           RENEW(t0,double,*nk*nsectors);
 485           RENEW(t1old,double,*nk*nsectors);
 486           RENEW(t1,double,*nk*nsectors);
 487           if(*nam>0) RENEW(iamt1,ITG,*nk*nsectors);
 488       }
 489       RENEW(nactdof,ITG,mt**nk*nsectors);
 490       if(*ntrans>0) RENEW(inotr,ITG,2**nk*nsectors);
 491       RENEW(kon,ITG,*nkon*nsectors);
 492       RENEW(ipkon,ITG,*ne*nsectors);
 493       for(i=*ne;i<*ne*nsectors;i++){ipkon[i]=-1;}
 494       RENEW(lakon,char,8**ne*nsectors);
 495       RENEW(ielmat,ITG,mi[2]**ne*nsectors);
 496       if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne*nsectors);
 497 //      RENEW(z,double,(long long)neq[1]*nev*nsectors/2);
 498 
 499       RENEW(nodeboun,ITG,*nboun*nsectors);
 500       RENEW(ndirboun,ITG,*nboun*nsectors);
 501       if(*nam>0) RENEW(iamboun,ITG,*nboun*nsectors);
 502       RENEW(xboun,double,*nboun*nsectors);
 503       RENEW(xbounold,double,*nboun*nsectors);
 504       RENEW(ikboun,ITG,*nboun*nsectors);
 505       RENEW(ilboun,ITG,*nboun*nsectors);
 506 
 507       NNEW(ipompcold,ITG,*nmpc);
 508       NNEW(nodempcold,ITG,3**mpcend);
 509       NNEW(coefmpcold,double,*mpcend);
 510       NNEW(labmpcold,char,20**nmpc);
 511       NNEW(ikmpcold,ITG,*nmpc);
 512       NNEW(ilmpcold,ITG,*nmpc);
 513 
 514       for(i=0;i<*nmpc;i++){ipompcold[i]=ipompc[i];}
 515       for(i=0;i<3**mpcend;i++){nodempcold[i]=nodempc[i];}
 516       for(i=0;i<*mpcend;i++){coefmpcold[i]=coefmpc[i];}
 517       for(i=0;i<20**nmpc;i++){labmpcold[i]=labmpc[i];}
 518       for(i=0;i<*nmpc;i++){ikmpcold[i]=ikmpc[i];}
 519       for(i=0;i<*nmpc;i++){ilmpcold[i]=ilmpc[i];}
 520       nmpcold=*nmpc;
 521       mpcendold=*mpcend;
 522 
 523       RENEW(ipompc,ITG,*nmpc*nsectors);
 524       RENEW(nodempc,ITG,3**mpcend*nsectors);
 525       RENEW(coefmpc,double,*mpcend*nsectors);
 526       RENEW(labmpc,char,20**nmpc*nsectors+1);
 527       RENEW(ikmpc,ITG,*nmpc*nsectors);
 528       RENEW(ilmpc,ITG,*nmpc*nsectors);
 529       RENEW(fmpc,double,*nmpc*nsectors);
 530 
 531       /* reallocating the fields for the nodes in which the
 532          solution has to be calculated */
 533 
 534       RENEW(imddof,ITG,neq[1]/2*nsectors);
 535       RENEW(imdnode,ITG,*nk*nsectors);
 536       RENEW(imdboun,ITG,*nboun*nsectors);
 537       RENEW(imdmpc,ITG,*nmpc*nsectors);
 538 
 539 //izdofNNEW(//      izdof,ITG,1);
 540 
 541       expand(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 542         ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
 543         nforc,nelemload,sideload,xload,nload,nactdof,neq,
 544         nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
 545         alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
 546         t0,ithermal,prestr,iprestr,vold,iperturb,sti,nzs,
 547         adb,aub,filab,eme,plicon,nplicon,plkcon,nplkcon,
 548         xstate,npmat_,matname,mi,ics,cs,mpcend,ncmat_,
 549         nstate_,mcs,nkon,ener,jobnamec,output,set,nset,istartset,
 550         iendset,ialset,nprint,prlab,prset,nener,trab,
 551         inotr,ntrans,ttime,fmpc,&nev,&z,iamboun,xbounold,
 552         &nsectors,nm,icol,irow,nzl,nam,ipompcold,nodempcold,coefmpcold,
 553         labmpcold,&nmpcold,xloadold,iamload,t1old,t1,iamt1,xstiff,&icole,&jqe,
 554         &irowe,isolver,nzse,&adbe,&aube,iexpl,
 555         ibody,xbody,nbody,cocon,ncocon,tieset,&ntie,imddof,&nmddof,
 556         imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,&izdof,&nzdof,
 557         &nherm,xmr,xmi,typeboun,ielprop,prop);
 558 
 559       RENEW(imddof,ITG,nmddof);
 560       RENEW(imdnode,ITG,nmdnode);
 561       RENEW(imdboun,ITG,nmdboun);
 562       RENEW(imdmpc,ITG,nmdmpc);
 563 
 564       SFREE(vold);
 565       NNEW(vold,double,mt**nk);
 566       SFREE(veold);
 567       NNEW(veold,double,mt**nk);
 568       RENEW(eme,double,6*mi[0]**ne);
 569 
 570       if(*nener==1) RENEW(ener,double,mi[0]**ne);
 571   }
 572 
 573   fclose(f1);
 574 
 575   /* allocating space for the friction coefficients */
 576 
 577   if(nherm==1){
 578       NNEW(fric,double,nev);
 579   }else{
 580       NNEW(fric,double,2*nev);
 581   }
 582 
 583   /* check whether there are dashpot elements */
 584 
 585   dashpot=0;
 586   for(i=0;i<*ne;i++){
 587       if(ipkon[i]==-1) continue;
 588       if(strcmp1(&lakon[i*8],"ED")==0){
 589           dashpot=1;break;}
 590   }
 591   if(dashpot){
 592 
 593       if(cyclicsymmetry){
 594           printf("*ERROR in steadystate: dashpots are not allowed in combination with cyclic symmetry\n");
 595           FORTRAN(stop,());
 596       }
 597 
 598       if(nherm!=1){
 599           printf("ERROR in steadystate: dashpots cannot be combined with non-Hermitian systems (in the present version of CalculiX)\n");
 600           FORTRAN(stop,());
 601       }
 602 
 603       /* cc is the reduced damping matrix (damping matrix mapped onto
 604          space spanned by eigenmodes) */
 605 
 606       NNEW(cc,double,nev*nev);
 607 /*      nev2=2*nev;
 608       NNEW(am,double,nev2*nev2);
 609       NNEW(bm,double,nev2);
 610       NNEW(ipiv,ITG,nev2);*/
 611   }
 612 
 613   NNEW(inum,ITG,*nk);
 614   strcpy1(&cflag[0],&filab[4],1);
 615   FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
 616           nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat));
 617 
 618   /* check whether integration point values are requested; if not,
 619      the stress fields do not have to be allocated */
 620 
 621   intpointvar=0;
 622   if(*ithermal<=1){
 623 
 624       /* mechanical */
 625 
 626       if((strcmp1(&filab[174],"S")==0)||
 627          (strcmp1(&filab[261],"E")==0)||
 628          (strcmp1(&filab[348],"RF")==0)||
 629          (strcmp1(&filab[435],"PEEQ")==0)||
 630          (strcmp1(&filab[522],"ENER")==0)||
 631          (strcmp1(&filab[609],"SDV")==0)||
 632          (strcmp1(&filab[1044],"ZZS")==0)||
 633          (strcmp1(&filab[1044],"ERR")==0)||
 634          (strcmp1(&filab[1479],"PHS")==0)||
 635          (strcmp1(&filab[1653],"MAXS")==0)||
 636          (strcmp1(&filab[2175],"CONT")==0)||
 637          (strcmp1(&filab[2262],"CELS")==0)) intpointvar=1;
 638       for(i=0;i<*nprint;i++){
 639           if((strcmp1(&prlab[6*i],"S")==0)||
 640              (strcmp1(&prlab[6*i],"E")==0)||
 641              (strcmp1(&prlab[6*i],"PEEQ")==0)||
 642              (strcmp1(&prlab[6*i],"ENER")==0)||
 643              (strcmp1(&prlab[6*i],"SDV")==0)||
 644              (strcmp1(&prlab[6*i],"RF")==0)) {intpointvar=1;break;}
 645       }
 646   }else{
 647 
 648       /* thermal */
 649 
 650       if((strcmp1(&filab[696],"HFL")==0)||
 651          (strcmp1(&filab[783],"RFL")==0)) intpointvar=1;
 652       for(i=0;i<*nprint;i++){
 653           if((strcmp1(&prlab[6*i],"HFL")==0)||
 654              (strcmp1(&prlab[6*i],"RFL")==0)) {intpointvar=1;break;}
 655       }
 656   }
 657    
 658   if(nfour<=0){
 659 
 660       /* harmonic excitation */
 661 
 662       NNEW(ikactmechr,ITG,neq[1]);
 663       NNEW(ikactmechi,ITG,neq[1]);
 664       nactmechr=0;nactmechi=0;
 665       
 666       /* result fields */
 667       
 668       if(intpointvar==1){
 669           NNEW(fn,double,mt**nk);
 670           NNEW(stnr,double,6**nk);
 671           NNEW(stni,double,6**nk);
 672           NNEW(stx,double,6*mi[0]**ne);
 673           NNEW(eei,double,6*mi[0]**ne);
 674 
 675           if(*ithermal>1) {NNEW(qfn,double,3**nk);
 676               NNEW(qfx,double,3*mi[0]**ne);}
 677           
 678           if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 679           if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 680           if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 681           
 682           if(*nener==1){
 683               NNEW(stiini,double,6*mi[0]**ne);
 684               NNEW(enerini,double,mi[0]**ne);}
 685       }
 686       
 687       /* determining the frequency data points */
 688       
 689       NNEW(freq,double,ndata*(nev+1));
 690       
 691       ndatatot=0.;
 692       freq[0]=fmin;
 693       if(fabs(fmax-fmin)<1.e-10){
 694           ndatatot=1;
 695       }else{
 696 
 697       /* copy the eigenvalues and sort them in ascending order
 698          (important for values from distinct nodal diameters */
 699 
 700           NNEW(e,double,nev);
 701           if(nherm==1){
 702               for(i=0;i<nev;i++){e[i]=sqrt(d[i]);}
 703           }else{
 704 
 705               /* for complex eigenvalues: sorting the real part */
 706 
 707               for(i=0;i<nev;i++){
 708                   e[i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/sqrt(2.);
 709               }
 710           }
 711           FORTRAN(dsort,(e,&idummy,&nev,&iflag));
 712           
 713           for(i=0;i<nev;i++){
 714               if(i!=0){
 715                   if(fabs(e[i]-e[i-1])<1.e-5){continue;}
 716               }
 717               if(e[i]>=fmin){
 718                   if(e[i]<=fmax){
 719                       for(j=1;j<ndata;j++){
 720                           y=-1.+2.*j/((double)(ndata-1));
 721                           if(fabs(y)<1.e-10){freq[ndatatot+j]=
 722                                  (freq[ndatatot]+e[i])/2.;}
 723                           else{
 724                               freq[ndatatot+j]=(freq[ndatatot]+e[i])/2.+
 725                                   (e[i]-freq[ndatatot])*pow(fabs(y),1./bias)*
 726                                   y/(2.*fabs(y));
 727                           }
 728                       }
 729                       ndatatot+=ndata-1;
 730                   }
 731                   else{break;}
 732               }
 733           }
 734 
 735           SFREE(e);
 736 
 737           for(j=1;j<ndata;j++){
 738               y=-1.+2.*j/((double)(ndata-1));
 739               if(fabs(y)<1.e-10){freq[ndatatot+j]=(freq[ndatatot]+fmax)/2.;}
 740               else{
 741                   freq[ndatatot+j]=(freq[ndatatot]+fmax)/2.+
 742                       (fmax-freq[ndatatot])*pow(fabs(y),1./bias)*
 743                       y/(2.*fabs(y));
 744               }
 745           }
 746           ndatatot+=ndata;
 747       }
 748       RENEW(freq,double,ndatatot);
 749       
 750       /*  check for nonzero SPC's */
 751       
 752       iprescribedboundary=0;
 753       for(i=0;i<*nboun;i++){
 754           if(fabs(xboun[i])>1.e-10){
 755               iprescribedboundary=1;
 756               nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 757               break;
 758           }
 759       }
 760 
 761       if((iprescribedboundary)&&(cyclicsymmetry)){
 762           printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
 763           FORTRAN(stop,());
 764       }
 765 
 766       /* calculating the damping coefficients = friction coefficient*2*eigenvalue */
 767       
 768       if(xmodal[10]<0){
 769           for(i=0;i<nev;i++){
 770               if(nherm==1){
 771                   if(sqrt(d[i])>(1.e-10)){
 772                       fric[i]=(alpham+betam*d[i]);
 773                   }
 774                   else {
 775                       printf("*WARNING in steadystate: one of the frequencies is zero\n");
 776                       printf("         no Rayleigh mass damping allowed\n");
 777                       fric[i]=0.;
 778                   }
 779               }else{
 780                   fric[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/
 781                            sqrt(2.);
 782                   fric[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])/
 783                            sqrt(2.);
 784                   if(d[2*i+1]<0.) fric[2*i+1]=-fric[2*i+1];
 785                   fric[2*i]=alpham+betam*fric[2*i];
 786                   fric[2*i+1]=betam*fric[2*i+1];
 787               }
 788           }
 789       }
 790       else{
 791           if(iprescribedboundary){
 792               printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with direct modal damping\n");
 793               FORTRAN(stop,());
 794           }
 795               
 796           /*copy the damping coefficients for every eigenfrequencie from xmodal[11....] */
 797           if(nev<(ITG)xmodal[10]){
 798               imax=nev;
 799               printf("*WARNING in steadystate: too many modal damping coefficients applied\n");
 800               printf("         damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
 801           }
 802           else{
 803               imax=(ITG)xmodal[10];
 804           }
 805           for(i=0; i<imax; i++){
 806               if(nherm==1){
 807                   fric[i]=2.*sqrt(d[i])*xmodal[11+i];    
 808               }else{
 809                   fric[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/
 810                            sqrt(2.);
 811                   fric[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])/
 812                            sqrt(2.);
 813                   if(d[2*i+1]<0.) fric[2*i+1]=-fric[2*i+1];
 814                   fric[2*i]=2.*fric[2*i]*xmodal[11+i];
 815                   fric[2*i+1]=2.*fric[2*i+1]*xmodal[11+i];
 816               } 
 817           }
 818           
 819       }
 820       
 821       /* check whether the loading is real or imaginary */
 822       
 823       NNEW(iphaseforc,ITG,*nforc);
 824       for (i=0;i<*nforc;i++){
 825           if(nodeforc[2*i+1]>=nsectors){
 826               iphaseforc[i]=1;
 827           }
 828       }
 829       
 830       NNEW(iphaseload,ITG,*nload);
 831       for (i=0;i<*nload;i++){
 832           if(nelemload[2*i+1]>=nsectors){
 833               iphaseload[i]=1;
 834           }
 835       }
 836       
 837       if(iprescribedboundary){
 838           NNEW(iphaseboun,ITG,*nboun);
 839           for (i=0;i<*nboun;i++){
 840               if(nodeboun[i]>*nk){
 841                   iphaseboun[i]=1;
 842                   nodeboun[i]=nodeboun[i]-*nk;
 843               }
 844           }
 845       }
 846       
 847       /* allocating actual loading fields */
 848       
 849       NNEW(xforcact,double,*nforc);
 850       NNEW(xforcr,double,*nforc);
 851       NNEW(xforci,double,*nforc);
 852       
 853       NNEW(xloadact,double,2**nload);
 854       NNEW(xloadr,double,2**nload);
 855       NNEW(xloadi,double,2**nload);
 856       
 857       NNEW(xbodyact,double,7**nbody);
 858       NNEW(xbodyr,double,7**nbody);
 859       NNEW(xbodyi,double,7**nbody);
 860       /* copying the rotation axis and/or acceleration vector */
 861       for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
 862       
 863       NNEW(xbounact,double,*nboun);
 864       
 865       if(*ithermal==1) NNEW(t1act,double,*nk);
 866       
 867       /* assigning the body forces to the elements */ 
 868 
 869       if(*nbody>0){
 870           ifreebody=*ne+1;
 871           NNEW(ipobody,ITG,2*ifreebody**nbody);
 872           for(k=1;k<=*nbody;k++){
 873               FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 874                                  iendset,ialset,&inewton,nset,&ifreebody,&k));
 875               RENEW(ipobody,ITG,2*(*ne+ifreebody));
 876           }
 877           RENEW(ipobody,ITG,2*(ifreebody-1));
 878       }
 879       
 880       NNEW(br,double,neq[1]); /* load rhs vector */
 881       NNEW(bi,double,neq[1]); /* load rhs vector */
 882       
 883       if(iprescribedboundary){
 884           NNEW(xbounr,double,*nboun);
 885           NNEW(xbouni,double,*nboun);
 886           
 887           NNEW(fr,double,neq[1]); /* force corresponding to real particular solution */
 888           NNEW(fi,double,neq[1]); /* force corresponding to imaginary particular solution */
 889           
 890           NNEW(ubr,double,neq[1]); /* real particular solution */
 891           NNEW(ubi,double,neq[1]); /* imaginary particular solution */
 892           
 893           NNEW(mubr,double,neq[1]); /* mass times real particular solution */
 894           NNEW(mubi,double,neq[1]); /* mass times imaginary particular solution */
 895       }
 896       
 897       NNEW(bjr,double,nev); /* real response modal decomposition */
 898       NNEW(bji,double,nev); /* imaginary response modal decomposition */
 899       
 900       NNEW(ampli,double,*nam); /* instantaneous amplitude */
 901       
 902       if(nherm==1){
 903           NNEW(aa,double,nev); /* modal coefficients of the real loading */
 904           NNEW(bb,double,nev); /* modal coefficients of the imaginary loading */
 905       }else{
 906           NNEW(aa,double,2*nev); /* modal coefficients of the real loading */
 907           NNEW(bb,double,2*nev); /* modal coefficients of the imaginary loading */
 908       }
 909     
 910       /* result fields */
 911       
 912       NNEW(vr,double,mt**nk);
 913       NNEW(vi,double,mt**nk);
 914       
 915       if(iprescribedboundary){
 916           
 917           /* LU decomposition of the stiffness matrix */
 918           
 919           if(*isolver==0){
 920 #ifdef SPOOLES
 921               spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
 922                       &symmetryflag,&inputformat,&nzs[2]);
 923 #else
 924               printf("*ERROR in steadystate: the SPOOLES library is not linked\n\n");
 925               FORTRAN(stop,());
 926 #endif
 927           }
 928           else if(*isolver==4){
 929 #ifdef SGI
 930               token=1;
 931               sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
 932 #else
 933               printf("*ERROR in steadystate: the SGI library is not linked\n\n");
 934               FORTRAN(stop,());
 935 #endif
 936           }
 937           else if(*isolver==5){
 938 #ifdef TAUCS
 939               tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
 940 #else
 941               printf("*ERROR in steadystate: the TAUCS library is not linked\n\n");
 942               FORTRAN(stop,());
 943 #endif
 944           }
 945           else if(*isolver==7){
 946 #ifdef PARDISO
 947               pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
 948                              &symmetryflag,&inputformat,jq,&nzs[2]);
 949 #else
 950               printf("*ERROR in steadystate: the PARDISO library is not linked\n\n");
 951               FORTRAN(stop,());
 952 #endif
 953           }
 954       }
 955       
 956       for(l=0;l<ndatatot;l=l+*jout){
 957           for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
 958           time=freq[l]/(2.*pi);
 959           timem=-time;
 960           ptime=time;
 961 
 962           /* calculating cc */
 963 
 964           if(dashpot){
 965               NNEW(adc,double,neq[1]);
 966               NNEW(auc,double,nzs[1]);
 967               FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
 968                 ndirboun,xboun,nboun,
 969                 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 970                 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 971                 adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
 972                 ikmpc,ilmpc,ikboun,ilboun,
 973                 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 974                 ielorien,norien,orab,ntmat_,
 975                 t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 976                 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 977                 xstiff,npmat_,&dtime,matname,mi,ncmat_,
 978                 ttime,&time,istep,&iinc,ibody,clearini,&mortar,springarea,
 979                 pslavsurf,pmastsurf,&reltime,&nasym));
 980 
 981               /*  zc = damping matrix * eigenmodes */
 982               
 983               NNEW(zc,double,neq[1]*nev);
 984               for(i=0;i<nev;i++){
 985                   FORTRAN(op,(&neq[1],&z[(long long)i*neq[1]],&zc[i*neq[1]],adc,auc,
 986                               jq,irow));
 987               }
 988               
 989               /* cc is the reduced damping matrix (damping matrix mapped onto
 990                  space spanned by eigenmodes) */
 991               
 992               for(i=0;i<nev*nev;i++){cc[i]=0.;}
 993               for(i=0;i<nev;i++){
 994                   for(j=0;j<=i;j++){
 995                       for(k=0;k<neq[1];k++){
 996                           cc[i*nev+j]+=z[(long long)j*neq[1]+k]*zc[i*neq[1]+k];
 997                       }
 998                   }
 999               }
1000               
1001               /* symmetric part of cc matrix */
1002               
1003               for(i=0;i<nev;i++){
1004                   for(j=i;j<nev;j++){
1005                       cc[i*nev+j]=cc[j*nev+i];
1006                   }
1007               }
1008               SFREE(zc);SFREE(adc);SFREE(auc);
1009           }
1010           
1011           /* calculating the instantaneous loads (forces, surface loading, 
1012              centrifugal and gravity loading or temperature) */
1013           
1014           FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
1015                xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1016                nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,namta,
1017                nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
1018                xbounold,xboun,xbounact,iamboun,nboun,nodeboun,ndirboun,
1019                nodeforc,ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,
1020                ikboun,ilboun,nelemload,sideload,mi,
1021                ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
1022                iendset,ialset,&ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
1023           
1024           /* real part of forces */
1025           
1026           for (i=0;i<*nforc;i++){
1027               xforcr[i]=xforcact[i]*(1-iphaseforc[i]);
1028           }
1029           
1030           for (i=0;i<*nload;i++){
1031               for(j=0;j<2;j++){
1032                   xloadr[2*i+j]=xloadact[2*i+j]*(1-iphaseload[i]);
1033               }
1034           }
1035 
1036           for(i=0;i<*nbody;i++){
1037               for(j=0;j<7;j++){
1038                   xbodyr[7*i+j]=xbodyact[7*i+j];
1039               }
1040               if(ibody[3*i+2]==2){
1041                   xbodyr[7*i]=0.;
1042               }
1043           }
1044           
1045           /* calculating the instantaneous loading vector */
1046           
1047           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1048                        ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcr,
1049                        nforc,nelemload,sideload,xloadr,nload,xbodyr,
1050                        ipobody,nbody,cgr,br,nactdof,&neq[1],nmethod,
1051                        ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1052                        alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1053                        t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1054                        nplicon,plkcon,nplkcon,
1055                        npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1056                        xbodyold,&reltime,veold,matname,mi,ikactmechr,
1057                        &nactmechr,ielprop,prop));
1058           
1059           /* real modal coefficients */
1060           
1061           if(!iprescribedboundary){
1062 
1063               if(nherm==1){
1064                   if(!cyclicsymmetry){
1065                       for(i=0;i<nev;i++){
1066                           i2=(long long)i*neq[1];
1067                           aa[i]=0.;
1068                           if(nactmechr<neq[1]/2){
1069                               for(j=0;j<nactmechr;j++){
1070                                   aa[i]+=z[i2+ikactmechr[j]]*br[ikactmechr[j]];
1071                               }
1072                           }else{
1073                               for(j=0;j<neq[1];j++){
1074                                   aa[i]+=z[i2+j]*br[j];
1075                               }
1076                           }
1077                       }
1078                   }else{
1079                       for(i=0;i<nev;i++){aa[i]=0.;}
1080                       for(j=0;j<nactmechr;j++){
1081                           for(i=0;i<nev;i++){
1082                               FORTRAN(nident,(izdof,&ikactmechr[j],&nzdof,&id));
1083                               if(id!=0){
1084                                   if(izdof[id-1]==ikactmechr[j]){
1085                                       aa[i]+=z[(long long)i*nzdof+id-1]*br[ikactmechr[j]];
1086                                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1087                               }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1088                           }
1089                       }
1090                   }
1091               }else{
1092                   if(!cyclicsymmetry){
1093                       for(i=0;i<nev;i++){
1094                           aa[2*i]=0.;aa[2*i+1]=0.;
1095                           if(nactmechr<neq[1]/2){
1096                               for(j=0;j<nactmechr;j++){
1097                                   aa[2*i]+=z[(long long)2*i*neq[1]+ikactmechr[j]]*br[ikactmechr[j]];
1098                                   aa[2*i+1]+=z[(long long)(2*i+1)*neq[1]+ikactmechr[j]]*br[ikactmechr[j]];
1099                               }
1100                           }else{
1101                               for(j=0;j<neq[1];j++){
1102                                   aa[2*i]+=z[(long long)2*i*neq[1]+j]*br[j];
1103                                   aa[2*i+1]+=z[(long long)(2*i+1)*neq[1]+j]*br[j];
1104                               }
1105                           }
1106                       }
1107                   }else{
1108                       for(i=0;i<nev;i++){aa[2*i]=0.;aa[2*i+1]=0.;}
1109                       for(j=0;j<nactmechr;j++){
1110                           for(i=0;i<nev;i++){
1111                               FORTRAN(nident,(izdof,&ikactmechr[j],&nzdof,&id));
1112                               if(id!=0){
1113                                   if(izdof[id-1]==ikactmechr[j]){
1114                                       aa[i]+=z[(long long)2*i*nzdof+id-1]*br[ikactmechr[j]];
1115                                       aa[i]+=z[(long long)(2*i+1)*nzdof+id-1]*br[ikactmechr[j]];
1116                                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1117                               }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1118                           }
1119                       }
1120                   }
1121               }
1122           
1123           }else{
1124           
1125               /* correction for nonzero SPC's */
1126               
1127               /* next statement makes sure that br is reset to zero at the
1128                  start of rhs.f */
1129               nactmechr=neq[1];
1130 
1131               /* real part of boundary conditions */
1132               
1133               for (i=0;i<*nboun;i++){
1134                   xbounr[i]=xbounact[i]*(1-iphaseboun[i]);
1135               }
1136               
1137               for(j=0;j<neq[1];j++){fr[j]=0.;ubr[j]=0.;}
1138               for(i=0;i<*nboun;i++){
1139                   ic=neq[1]+i;
1140                   for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
1141                       ir=irow[j]-1;
1142                       fr[ir]=fr[ir]-au[j]*xbounr[i];
1143                       ubr[ir]=fr[ir];
1144                   }
1145               }
1146               if(*isolver==0){
1147 #ifdef SPOOLES
1148                   spooles_solve(ubr,&neq[1]);
1149 #endif
1150               }
1151               else if(*isolver==4){
1152 #ifdef SGI
1153                   sgi_solve(ubr,token);
1154 #endif
1155               }
1156               else if(*isolver==5){
1157 #ifdef TAUCS
1158                   tau_solve(ubr,&neq[1]);
1159 #endif
1160               }
1161               else if(*isolver==7){
1162 #ifdef PARDISO
1163                   pardiso_solve(ubr,&neq[1],&symmetryflag);
1164 #endif
1165               }
1166               FORTRAN(op,(&neq[1],ubr,mubr,adb,aub,jq,irow));
1167           }
1168           
1169           /* imaginary part of forces */
1170           
1171           for (i=0;i<*nforc;i++){
1172               xforci[i]=xforcact[i]*iphaseforc[i];
1173           }
1174           
1175           for (i=0;i<*nload;i++){
1176               for(j=0;j<2;j++){
1177                   xloadi[2*i+j]=xloadact[2*i+j]*iphaseload[i];
1178               }
1179           }
1180           
1181           for(i=0;i<*nbody;i++){
1182               for(j=0;j<7;j++){
1183                   xbodyi[7*i+j]=xbodyact[7*i+j];
1184               }
1185               if(ibody[3*i+2]==1){
1186                   xbodyi[7*i]=0.;
1187               }
1188           }
1189           
1190           /* calculating the instantaneous loading vector */
1191           
1192           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1193                        ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforci,
1194                        nforc,nelemload,sideload,xloadi,nload,xbodyi,
1195                        ipobody,nbody,cgr,bi,nactdof,&neq[1],nmethod,
1196                        ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1197                        alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1198                        t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1199                        nplicon,plkcon,nplkcon,
1200                        npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1201                        xbodyold,&reltime,veold,matname,mi,ikactmechi,
1202                        &nactmechi,ielprop,prop));
1203           
1204           /* imaginary modal coefficients */
1205           
1206           if(!iprescribedboundary){
1207 
1208               if(nherm==1){
1209                   if(!cyclicsymmetry){
1210                       for(i=0;i<nev;i++){
1211                           i2=(long long)i*neq[1];
1212                           bb[i]=0.;
1213                           if(nactmechi<neq[1]/2){
1214                               for(j=0;j<nactmechi;j++){
1215                                   bb[i]+=z[i2+ikactmechi[j]]*bi[ikactmechi[j]];
1216                               }
1217                           }else{
1218                               for(j=0;j<neq[1];j++){
1219                                   bb[i]+=z[i2+j]*bi[j];
1220                               }
1221                           }
1222                       }
1223                   }else{
1224                       for(i=0;i<nev;i++){bb[i]=0.;}
1225                       for(j=0;j<nactmechi;j++){
1226                           for(i=0;i<nev;i++){
1227                               FORTRAN(nident,(izdof,&ikactmechi[j],&nzdof,&id));
1228                               if(id!=0){
1229                                   if(izdof[id-1]==ikactmechi[j]){
1230                                       bb[i]+=z[(long long)i*nzdof+id-1]*bi[ikactmechi[j]];
1231                                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1232                               }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1233                           }
1234                       }
1235                   }
1236               }else{
1237                   if(!cyclicsymmetry){
1238                       for(i=0;i<nev;i++){
1239                           bb[2*i]=0.;bb[2*i+1]=0.;
1240                           if(nactmechi<neq[1]/2){
1241                               for(j=0;j<nactmechi;j++){
1242                                   bb[2*i]+=z[(long long)2*i*neq[1]+ikactmechi[j]]*bi[ikactmechi[j]];
1243                                   bb[2*i+1]+=z[(long long)(2*i+1)*neq[1]+ikactmechi[j]]*bi[ikactmechi[j]];
1244                               }
1245                           }else{
1246                               for(j=0;j<neq[1];j++){
1247                                   bb[2*i]+=z[(long long)2*i*neq[1]+j]*bi[j];
1248                                   bb[2*i+1]+=z[(long long)(2*i+1)*neq[1]+j]*bi[j];
1249                               }
1250                           }
1251                       }
1252                   }else{
1253                       for(i=0;i<nev;i++){bb[2*i]=0.;bb[2*i+1]=0.;}
1254                       for(j=0;j<nactmechi;j++){
1255                           for(i=0;i<nev;i++){
1256                               FORTRAN(nident,(izdof,&ikactmechi[j],&nzdof,&id));
1257                               if(id!=0){
1258                                   if(izdof[id-1]==ikactmechi[j]){
1259                                       bb[2*i]+=z[(long long)2*i*nzdof+id-1]*bi[ikactmechi[j]];
1260                                       bb[2*i+1]+=z[(long long)(2*i+1)*nzdof+id-1]*bi[ikactmechi[j]];
1261                                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1262                               }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1263                           }
1264                       }
1265                   }
1266               }
1267 
1268           }else{
1269           
1270           /* correction for nonzero SPC's */
1271               
1272               /* next statement makes sure that bi is reset to zero at the
1273                  start of rhs.f */
1274 
1275               nactmechi=neq[1];
1276               
1277               /* imaginary part of boundary conditions */
1278               
1279               for (i=0;i<*nboun;i++){
1280                   xbouni[i]=xbounact[i]*iphaseboun[i];
1281               }
1282               
1283               for(j=0;j<neq[1];j++){fi[j]=0.;ubi[j]=0.;}
1284               for(i=0;i<*nboun;i++){
1285                   ic=neq[1]+i;
1286                   for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
1287                       ir=irow[j]-1;
1288                       fi[ir]=fi[ir]-au[j]*xbouni[i];
1289                       ubi[ir]=fi[ir];
1290                   }
1291               }
1292               if(*isolver==0){
1293 #ifdef SPOOLES
1294                   spooles_solve(ubi,&neq[1]);
1295 #endif
1296               }
1297               else if(*isolver==4){
1298 #ifdef SGI
1299                   sgi_solve(ubi,token);
1300 #endif
1301               }
1302               else if(*isolver==5){
1303 #ifdef TAUCS
1304                   tau_solve(ubi,&neq[1]);
1305 #endif
1306               }
1307               else if(*isolver==7){
1308 #ifdef PARDISO
1309                   pardiso_solve(ubi,&neq[1],&symmetryflag);
1310 #endif
1311               }
1312               FORTRAN(op,(&neq[1],ubi,mubi,adb,aub,jq,irow));
1313           
1314               /* correction for prescribed boundary conditions */
1315           
1316               for(i=0;i<neq[1];i++){
1317                   br[i]+=freq[l]*(freq[l]*mubr[i]+alpham*mubi[i]+betam*fi[i]);
1318                   bi[i]+=freq[l]*(freq[l]*mubi[i]-alpham*mubr[i]-betam*fr[i]);
1319               }
1320           
1321               /* real and imaginary modal coefficients */
1322 
1323               for(i=0;i<nev;i++){
1324                   aa[i]=0.;
1325                   for(j=0;j<neq[1];j++){
1326                       aa[i]+=z[(long long)i*neq[1]+j]*br[j];
1327                   }
1328               }
1329               
1330               for(i=0;i<nev;i++){
1331                   bb[i]=0.;
1332                   for(j=0;j<neq[1];j++){
1333                       bb[i]+=z[(long long)i*neq[1]+j]*bi[j];
1334                   }
1335               }
1336               
1337           }
1338           
1339           /* calculating the modal coefficients */
1340           
1341           if(nherm==1){
1342               if(dashpot==0){
1343                   for(i=0;i<nev;i++){
1344                       dd=pow(d[i]-pow(freq[l],2),2)+
1345                           pow(fric[i],2)*pow(freq[l],2);
1346                       bjr[i]=(aa[i]*(d[i]-freq[l]*freq[l])+
1347                               bb[i]*fric[i]*freq[l])/dd;
1348                       bji[i]=(bb[i]*(d[i]-freq[l]*freq[l])-
1349                               aa[i]*fric[i]*freq[l])/dd;
1350                   }
1351                   /*    printf("old l=%" ITGFORMAT ",bjr=%f,bji=%f\n",l,bjr[0],bji[0]);*/
1352               }else{
1353                   nev2=2*nev;
1354                   NNEW(am,double,nev2*nev2);
1355                   NNEW(bm,double,nev2);
1356                   NNEW(ipiv,ITG,nev2);
1357                   
1358                   for(i=0;i<nev2;i++){
1359                       for(j=0;j<nev2;j++){
1360                           am[i*nev2+j]=0.;
1361                       }
1362                       bm[i]=0.;
1363                   }
1364                   for(i=0;i<nev;i++){
1365                       am[i*nev2+i]=d[i]-freq[l]*freq[l];
1366                       am[(i+nev)*nev2+i]=-fric[i]*freq[l];
1367                       bm[i]=aa[i];
1368                       am[i*nev2+nev+i]=-am[(i+nev)*nev2+i];
1369                       am[(i+nev)*nev2+nev+i]=am[i*nev2+i];
1370                       bm[nev+i]=bb[i];
1371                       for(j=0;j<nev;j++){
1372                           am[(j+nev)*nev2+i]=am[(j+nev)*nev2+i]
1373                               -cc[i*nev+j]*freq[l];
1374                           am[j*nev2+nev+i]=am[j*nev2+nev+i]
1375                               +cc[i*nev+j]*freq[l];
1376                       }
1377                   }
1378                   
1379                   /* solving the system of equations */
1380                   
1381                   FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
1382                   if(info!=0){
1383                       printf("*ERROR in steadystate: fatal termination of dgesv\n");
1384                       printf("       info=%" ITGFORMAT "\n",info);
1385                       FORTRAN(stop,());
1386                   }
1387                   
1388                   /* storing the solution in bjr and bji */
1389                   
1390                   for(i=0;i<nev;i++){
1391                       bjr[i]=bm[i];
1392                       bji[i]=bm[nev+i];
1393                   }
1394 
1395                   SFREE(am);SFREE(bm);SFREE(ipiv);
1396               }
1397           }else{
1398               nev2=2*nev;
1399               NNEW(am,double,nev2*nev2);
1400               NNEW(bm,double,nev2);
1401               NNEW(ipiv,ITG,nev2);
1402                   
1403               NNEW(ax,double,nev);
1404               NNEW(bx,double,nev);
1405               for(i=0;i<nev;i++){
1406                   ax[i]=-pow(freq[l],2)-freq[l]*fric[2*i+1]+d[2*i];
1407                   bx[i]=-freq[l]*fric[2*i]-d[2*i+1];
1408               }
1409               for(i=0;i<nev;i++){
1410                   for(j=0;j<nev;j++){
1411                       am[j*nev2+i]=xmr[j*nev+i]*ax[j]+xmi[j*nev+i]*bx[j];
1412                       am[(j+nev)*nev2+i]=xmr[j*nev+i]*bx[j]-xmi[j*nev+i]*bx[j];
1413                       am[j*nev2+nev+i]=xmi[j*nev+i]*ax[j]-xmr[j*nev+i]*bx[j];
1414                       am[(j+nev)*nev2+nev+i]=xmi[j*nev+i]*bx[j]+xmr[j*nev+i]*ax[j];
1415                   }
1416                   bm[i]=aa[i]-bb[nev+i];
1417                   bm[nev+i]=bb[i]+aa[nev+i];
1418               }
1419               SFREE(ax);SFREE(bx);
1420                   
1421               /* solving the system of equations */
1422               
1423               FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
1424               if(info!=0){
1425                   printf("*ERROR in steadystate: fatal termination of dgesv\n");
1426                   printf("       info=%" ITGFORMAT "\n",info);
1427                   FORTRAN(stop,());
1428               }
1429               
1430               /* storing the solution in bjr and bji */
1431               
1432               for(i=0;i<nev;i++){
1433                   bjr[i]=bm[i];
1434                   bji[i]=bm[nev+i];
1435               }
1436               
1437               SFREE(am);SFREE(bm);SFREE(ipiv);
1438           }
1439 
1440           /* storing the participation factors */
1441 
1442           if(nherm==1){
1443               NNEW(eig,double,nev);
1444               for(i=0;i<nev;i++){
1445                   eig[i]=sqrt(d[i]);
1446               }
1447           }else{
1448               
1449               NNEW(eig,double,2*nev);
1450               for(i=0;i<nev;i++){
1451                   eig[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])
1452                       /sqrt(2.);
1453                   eig[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])
1454                       /sqrt(2.);
1455                   if(d[2*i+1]<0.) eig[2*i+1]=-eig[2*i+1];
1456               }
1457           }
1458           
1459           mode=0;
1460           FORTRAN(writepf,(eig,bjr,bji,&time,&nev,&mode,&nherm));
1461           SFREE(eig);
1462       
1463           /* calculating the real response */
1464       
1465           if(iprescribedboundary){
1466               if(nmdnode==0){
1467                   memcpy(&br[0],&ubr[0],sizeof(double)*neq[1]);
1468               }else{
1469                   for(i=0;i<nmddof;i++){
1470                       br[imddof[i]]=ubr[imddof[i]];
1471                   }
1472               }
1473           }
1474           else{
1475               if(nmdnode==0){
1476                   DMEMSET(br,0,neq[1],0.);
1477               }else{
1478                   for(i=0;i<nmddof;i++){
1479                       br[imddof[i]]=0.;
1480                   }
1481               }
1482           }
1483           
1484           if(!cyclicsymmetry){
1485               if(nmdnode==0){
1486                   for(i=0;i<neq[1];i++){
1487                       for(j=0;j<nev;j++){
1488                           br[i]+=bjr[j]*z[(long long)j*neq[1]+i];
1489                       }
1490                   }
1491               }else{
1492                   for(i=0;i<nmddof;i++){
1493                       for(j=0;j<nev;j++){
1494                           br[imddof[i]]+=bjr[j]*z[(long long)j*neq[1]+imddof[i]];
1495                       }
1496                   }
1497               }
1498           }else{
1499               for(i=0;i<nmddof;i++){
1500                   FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1501                   if(id!=0){
1502                       if(izdof[id-1]==imddof[i]){
1503                           for(j=0;j<nev;j++){
1504                               br[imddof[i]]+=bjr[j]*z[(long long)j*nzdof+id-1];
1505                           }
1506                       }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1507                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1508               }
1509           }
1510           
1511           if(nmdnode==0){
1512               DMEMSET(vr,0,mt**nk,0.);
1513           }else{
1514               for(jj=0;jj<nmdnode;jj++){
1515                   i=imdnode[jj]-1;
1516                   for(j=1;j<4;j++){
1517                       vr[mt*i+j]=0.;
1518                   }
1519               }
1520           }
1521           
1522           /* calculating the imaginary response */
1523       
1524           if(iprescribedboundary){
1525               if(nmdnode==0){
1526                   memcpy(&bi[0],&ubi[0],sizeof(double)*neq[1]);
1527               }else{
1528                   for(i=0;i<nmddof;i++){
1529                       bi[imddof[i]]=ubi[imddof[i]];
1530                   }
1531               }
1532           }
1533           else{
1534               if(nmdnode==0){
1535                   DMEMSET(bi,0,neq[1],0.);
1536               }else{
1537                   for(i=0;i<nmddof;i++){
1538                       bi[imddof[i]]=0.;
1539                   }
1540               }
1541           }
1542           
1543           if(!cyclicsymmetry){
1544               if(nmdnode==0){
1545                   for(i=0;i<neq[1];i++){
1546                       for(j=0;j<nev;j++){
1547                           bi[i]+=bji[j]*z[(long long)j*neq[1]+i];
1548                       }
1549                   }
1550               }else{
1551                   for(i=0;i<nmddof;i++){
1552                       for(j=0;j<nev;j++){
1553                           bi[imddof[i]]+=bji[j]*z[(long long)j*neq[1]+imddof[i]];
1554                       }
1555                   }
1556               }
1557           }else{
1558               for(i=0;i<nmddof;i++){
1559                   FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1560                   if(id!=0){
1561                       if(izdof[id-1]==imddof[i]){
1562                           for(j=0;j<nev;j++){
1563                               bi[imddof[i]]+=bji[j]*z[(long long)j*nzdof+id-1];
1564                           }
1565                       }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1566                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1567               }
1568           }
1569           
1570           
1571           if(nmdnode==0){
1572               DMEMSET(vi,0,mt**nk,0.);
1573           }else{
1574               for(jj=0;jj<nmdnode;jj++){
1575                   i=imdnode[jj]-1;
1576                   for(j=1;j<4;j++){
1577                       vi[mt*i+j]=0.;
1578                   }
1579               }
1580           }
1581 
1582           /* real response */
1583 
1584           if(iprescribedboundary){
1585       
1586               /* calculating displacements/temperatures */
1587       
1588               FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
1589                     ndirboun,xbounr,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1590                     br,bi,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1591                     imdmpc,&nmdmpc,nmethod,&timem));
1592 
1593               results(co,nk,kon,ipkon,lakon,ne,vr,stnr,inum,
1594                       stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1595                       ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1596                       ithermal,prestr,iprestr,filab,eme,emn,een,
1597                       iperturb,f,fn,nactdof,&iout,qa,
1598                       vold,br,nodeboun,ndirboun,xbounr,nboun,
1599                       ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1600                       veold,accold,&bet,&gam,&dtime,&time,&xnull,
1601                       plicon,nplicon,plkcon,nplkcon,
1602                       xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1603                       &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1604                       enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1605                       set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1606                       qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1607                       ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1608                       xforc,nforc,thicke,shcon,nshcon,
1609                       sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1610                       &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
1611                       islavsurf,ielprop,prop);}
1612           else{
1613       
1614               /* calculating displacements/temperatures */
1615       
1616               FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
1617                     ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1618                     br,bi,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1619                     imdmpc,&nmdmpc,nmethod,&timem));
1620 
1621               results(co,nk,kon,ipkon,lakon,ne,vr,stnr,inum,
1622                       stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1623                       ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1624                       ithermal,prestr,iprestr,filab,eme,emn,een,
1625                       iperturb,f,fn,nactdof,&iout,qa,
1626                       vold,br,nodeboun,ndirboun,xbounact,nboun,
1627                       ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1628                       veold,accold,&bet,&gam,&dtime,&time,&xnull,
1629                       plicon,nplicon,plkcon,nplkcon,
1630                       xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1631                       &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1632                       enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1633                       set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1634                       qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1635                       ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1636                       xforc,nforc,thicke,shcon,nshcon,
1637                       sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1638                       &mortar,islavact,cdn,islavnode,nslavnode,&ntie,
1639                       clearini,islavsurf,ielprop,prop);
1640               
1641               if(nmdnode==0){
1642                   DMEMSET(br,0,neq[1],0.);
1643               }else{
1644                   for(i=0;i<nmddof;i++){
1645                       br[imddof[i]]=0.;
1646                   }
1647               }
1648           }
1649           
1650           (*kode)++;
1651           
1652           mode=-1;
1653           if(strcmp1(&filab[1044],"ZZS")==0){
1654               NNEW(neigh,ITG,40**ne);
1655               NNEW(ipneigh,ITG,*nk);
1656           }
1657 
1658           frd(co,&nkg,kon,ipkon,lakon,&neg,vr,stnr,inum,nmethod,
1659             kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1660             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1661             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1662             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&neg,
1663             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1664             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1665 
1666           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1667 
1668           /* imaginary response */
1669 
1670           if(iprescribedboundary){
1671       
1672               /* calculating displacements/temperatures */
1673       
1674               FORTRAN(dynresults,(nk,vi,ithermal,nactdof,vold,nodeboun,
1675                     ndirboun,xbouni,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1676                     bi,br,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1677                     imdmpc,&nmdmpc,nmethod,&time));
1678 
1679               results(co,nk,kon,ipkon,lakon,ne,vi,stni,inum,
1680                       stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1681                       ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1682                       ithermal,prestr,iprestr,filab,eme,emn,een,
1683                       iperturb,f,fn,nactdof,&iout,qa,
1684                       vold,bi,nodeboun,ndirboun,xbouni,nboun,
1685                       ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1686                       veold,accold,&bet,&gam,&dtime,&time,&xnull,
1687                       plicon,nplicon,plkcon,nplkcon,
1688                       xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1689                       &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1690                       enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1691                       set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1692                       qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1693                       ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1694                       xforc,nforc,thicke,shcon,nshcon,
1695                       sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1696                       &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
1697                       islavsurf,ielprop,prop);}
1698           else{ 
1699       
1700               /* calculating displacements/temperatures */
1701       
1702               FORTRAN(dynresults,(nk,vi,ithermal,nactdof,vold,nodeboun,
1703                     ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1704                     bi,br,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1705                     imdmpc,&nmdmpc,nmethod,&time));
1706 
1707               results(co,nk,kon,ipkon,lakon,ne,vi,stni,inum,
1708                       stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1709                       ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1710                       ithermal,prestr,iprestr,filab,eme,emn,een,
1711                       iperturb,f,fn,nactdof,&iout,qa,
1712                       vold,bi,nodeboun,ndirboun,xbounact,nboun,
1713                       ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1714                       veold,accold,&bet,&gam,&dtime,&time,&xnull,
1715                       plicon,nplicon,plkcon,nplkcon,
1716                       xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1717                       &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1718                       enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1719                       set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1720                       qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1721                       ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1722                       xforc,nforc,thicke,shcon,nshcon,
1723                       sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1724                       &mortar,islavact,cdn,islavnode,nslavnode,&ntie,
1725                       clearini,islavsurf,ielprop,prop);
1726 
1727               if(nmdnode==0){
1728                   DMEMSET(bi,0,neq[1],0.);
1729               }else{
1730                   for(i=0;i<nmddof;i++){
1731                       bi[imddof[i]]=0.;
1732                   }
1733               }
1734           }
1735           
1736           /* calculating the magnitude and phase */
1737           
1738           if(strcmp1(&filab[870],"PU")==0){
1739               
1740               constant=180./pi;
1741               NNEW(va,double,mt**nk);
1742               NNEW(vp,double,mt**nk);
1743               
1744               if(*ithermal<=1){
1745                   if(nmdnode==0){
1746                       for(i=0;i<*nk;i++){
1747                           for(j=1;j<4;j++){
1748                               vreal=vr[mt*i+j];
1749                               va[mt*i+j]=sqrt(vr[mt*i+j]*vr[mt*i+j]+vi[mt*i+j]*vi[mt*i+j]);
1750                               if(fabs(vreal)<1.e-10){
1751                                   if(vi[mt*i+j]>0.){vp[mt*i+j]=90.;}
1752                                   else{vp[mt*i+j]=-90.;}
1753                               }
1754                               else{
1755                                   vp[mt*i+j]=atan(vi[mt*i+j]/vreal)*constant;
1756                                   if(vreal<0.) vp[mt*i+j]+=180.;
1757                               }
1758                           }
1759                       }
1760                   }else{
1761                       for(jj=0;jj<nmdnode;jj++){
1762                           i=imdnode[jj]-1;
1763                           for(j=1;j<4;j++){
1764                               vreal=vr[mt*i+j];
1765                               va[mt*i+j]=sqrt(vr[mt*i+j]*vr[mt*i+j]+vi[mt*i+j]*vi[mt*i+j]);
1766                               if(fabs(vreal)<1.e-10){
1767                                   if(vi[mt*i+j]>0.){vp[mt*i+j]=90.;}
1768                                   else{vp[mt*i+j]=-90.;}
1769                               }
1770                               else{
1771                                   vp[mt*i+j]=atan(vi[mt*i+j]/vreal)*constant;
1772                                   if(vreal<0.) vp[mt*i+j]+=180.;
1773                               }
1774                           }
1775                       }
1776                   }
1777               }
1778               else{
1779                   if(nmdnode==0){
1780                       for(i=0;i<*nk;i++){
1781                           vreal=vr[mt*i];
1782                           va[mt*i]=sqrt(vr[mt*i]*vr[mt*i]+vi[mt*i]*vi[mt*i]);
1783                           if(fabs(vreal)<1.e-10){
1784                               if(vi[mt*i]>0){vp[mt*i]=90.;}
1785                               else{vp[mt*i]=-90.;}
1786                           }
1787                           else{
1788                               vp[mt*i]=atan(vi[mt*i]/vreal)*constant;
1789                               if(vreal<0.) vp[mt*i]+=180.;
1790                           }
1791                       }
1792                   }else{
1793                       for(jj=0;jj<nmdnode;jj++){
1794                           i=imdnode[jj]-1;
1795                           vreal=vr[mt*i];
1796                           va[mt*i]=sqrt(vr[mt*i]*vr[mt*i]+vi[mt*i]*vi[mt*i]);
1797                           if(fabs(vreal)<1.e-10){
1798                               if(vi[mt*i]>0){vp[mt*i]=90.;}
1799                               else{vp[mt*i]=-90.;}
1800                           }
1801                           else{
1802                               vp[mt*i]=atan(vi[mt*i]/vreal)*constant;
1803                               if(vreal<0.) vp[mt*i]+=180.;
1804                           }
1805                       }
1806                   }
1807               }
1808           }
1809           
1810           if(strcmp1(&filab[1479],"PHS")==0){
1811               
1812               constant=180./pi;
1813               NNEW(stna,double,6**nk);
1814               NNEW(stnp,double,6**nk);
1815               
1816               if(*ithermal<=1){
1817                   if(nmdnode==0){
1818                       for(i=0;i<*nk;i++){
1819                           for(j=0;j<6;j++){
1820                               vreal=stnr[6*i+j];
1821                               stna[6*i+j]=sqrt(stnr[6*i+j]*stnr[6*i+j]+stni[6*i+j]*stni[6*i+j]);
1822                               if(fabs(vreal)<1.e-10){
1823                                   if(stni[6*i+j]>0.){stnp[6*i+j]=90.;}
1824                                   else{stnp[6*i+j]=-90.;}
1825                               }
1826                               else{
1827                                   stnp[6*i+j]=atan(stni[6*i+j]/vreal)*constant;
1828                                   if(vreal<0.) stnp[6*i+j]+=180.;
1829                               }
1830                           }
1831                       }
1832                   }else{
1833                       for(jj=0;jj<nmdnode;jj++){
1834                           i=imdnode[jj]-1;
1835                           for(j=0;j<6;j++){
1836                               vreal=stnr[6*i+j];
1837                               stna[6*i+j]=sqrt(stnr[6*i+j]*stnr[6*i+j]+stni[6*i+j]*stni[6*i+j]);
1838                               if(fabs(vreal)<1.e-10){
1839                                   if(stni[6*i+j]>0.){stnp[6*i+j]=90.;}
1840                                   else{stnp[6*i+j]=-90.;}
1841                               }
1842                               else{
1843                                   stnp[6*i+j]=atan(stni[6*i+j]/vreal)*constant;
1844                                   if(vreal<0.) stnp[6*i+j]+=180.;
1845                               }
1846                           }
1847                       }
1848                   }
1849               }
1850           }
1851  
1852           (*kode)++;
1853           mode=0;
1854           
1855           if(strcmp1(&filab[1044],"ZZS")==0){
1856               NNEW(neigh,ITG,40**ne);
1857               NNEW(ipneigh,ITG,*nk);
1858           }
1859 
1860           frd(co,&nkg,kon,ipkon,lakon,&neg,vi,stni,inum,nmethod,
1861             kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1862             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1863             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1864             mi,stx,va,vp,stna,stnp,vmax,stnmax,&ngraph,veold,ener,&neg,
1865             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1866             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1867 
1868           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1869 
1870           SFREE(va);SFREE(vp);SFREE(stna);SFREE(stnp);
1871           
1872       }
1873       
1874       /* restoring the imaginary loading */
1875       
1876       SFREE(iphaseforc);SFREE(xforcr);SFREE(xforci);
1877 
1878       SFREE(iphaseload);SFREE(xloadr);SFREE(xloadi);
1879       
1880       SFREE(xbodyr);SFREE(xbodyi);
1881       
1882       if(iprescribedboundary){
1883           for (i=0;i<*nboun;i++){
1884               if(iphaseboun[i]==1){
1885                   nodeboun[i]=nodeboun[i]+*nk;
1886               }
1887           }
1888           SFREE(iphaseboun);
1889       }
1890       
1891       /* updating the loading at the end of the step; 
1892          important in case the amplitude at the end of the step
1893          is not equal to one */
1894       
1895       for(k=0;k<*nboun;++k){xboun[k]=xbounact[k];}
1896       for(k=0;k<*nforc;++k){xforc[k]=xforcact[k];}
1897       for(k=0;k<2**nload;++k){xload[k]=xloadact[k];}
1898       for(k=0;k<7**nbody;k=k+7){xbody[k]=xbodyact[k];}
1899       if(*ithermal==1){
1900           for(k=0;k<*nk;++k){t1[k]=t1act[k];}
1901       }
1902       
1903       SFREE(br);SFREE(bi);SFREE(bjr);SFREE(bji),SFREE(freq);
1904       SFREE(xforcact);SFREE(xloadact);SFREE(xbounact);SFREE(aa);SFREE(bb);
1905       SFREE(ampli);SFREE(xbodyact);SFREE(vr);SFREE(vi);if(*nbody>0) SFREE(ipobody);
1906       
1907       if(*ithermal==1) SFREE(t1act);
1908       
1909       if(iprescribedboundary){
1910           if(*isolver==0){
1911 #ifdef SPOOLES
1912               spooles_cleanup();
1913 #endif
1914           }
1915           else if(*isolver==4){
1916 #ifdef SGI
1917               sgi_cleanup(token);
1918 #endif
1919           }
1920           else if(*isolver==5){
1921 #ifdef TAUCS
1922               tau_cleanup();
1923 #endif
1924           }
1925           else if(*isolver==7){
1926 #ifdef PARDISO
1927               pardiso_cleanup(&neq[1],&symmetryflag);
1928 #endif
1929           }
1930           SFREE(xbounr);SFREE(xbouni);SFREE(fr);SFREE(fi);SFREE(ubr);SFREE(ubi);
1931           SFREE(mubr);SFREE(mubi);
1932       }
1933 
1934       SFREE(ikactmechr);SFREE(ikactmechi);
1935       
1936       if(intpointvar==1){
1937           SFREE(fn);
1938           SFREE(stnr);SFREE(stni);SFREE(stx);SFREE(eei);
1939 
1940           if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
1941           
1942           if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
1943           if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1944           if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
1945           
1946           if(*nener==1){SFREE(stiini);SFREE(enerini);}
1947       }
1948 
1949   }else{
1950 
1951       /* steady state response to a nonharmonic periodic loading */
1952 
1953       NNEW(ikactmech,ITG,neq[1]);
1954       nactmech=0;
1955       
1956       NNEW(xforcact,double,nfour**nforc);
1957       NNEW(xloadact,double,nfour*2**nload);
1958       NNEW(xbodyact,double,nfour*7**nbody);
1959       NNEW(xbounact,double,nfour**nboun);
1960       NNEW(xbounacttime,double,nfour**nboun);
1961       if(*ithermal==1) NNEW(t1act,double,*nk);
1962 
1963       NNEW(r,double,nfour);
1964       NNEW(wsave,double,2*nfour);
1965       NNEW(isave,ITG,15);
1966       
1967       /*  check for nonzero SPC's */
1968       
1969       iprescribedboundary=0;
1970       for(i=0;i<*nboun;i++){
1971           if(fabs(xboun[i])>1.e-10){
1972               iprescribedboundary=1;
1973               break;
1974           }
1975       }
1976 
1977       if((iprescribedboundary)&&(cyclicsymmetry)){
1978           printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
1979           FORTRAN(stop,());
1980       }
1981 
1982       /* calculating the damping coefficients = friction coefficient*2*eigenvalue */
1983       
1984       if(xmodal[10]<0){
1985           for(i=0;i<nev;i++){
1986               if(sqrt(d[i])>(1.e-10)){
1987                   fric[i]=(alpham+betam*d[i]);
1988               }
1989               else {
1990                   printf("*WARNING in steadystate: one of the frequencies is zero\n");
1991                   printf("         no Rayleigh mass damping allowed\n");
1992                   fric[i]=0.;
1993               }
1994           }
1995       }
1996       else{
1997           if(iprescribedboundary){
1998               printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with direct modal damping\n");
1999               FORTRAN(stop,());
2000           }
2001               
2002           /*copy the damping coefficients for every eigenfrequencie from xmodal[11....] */
2003           if(nev<(ITG)xmodal[10]){
2004               imax=nev;
2005               printf("*WARNING in steadystate: too many modal damping coefficients applied\n");
2006               printf("         damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
2007           }
2008           else{
2009               imax=(ITG)xmodal[10];
2010           }
2011           for(i=0; i<imax; i++){
2012               fric[i]=2.*sqrt(d[i])*xmodal[11+i];     
2013           }
2014           
2015       }
2016 
2017       /* determining the load time history */
2018       
2019       NNEW(ampli,double,*nam); /* instantaneous amplitude */
2020 
2021       for(l=0;l<nfour;l++){
2022 
2023           time=tmin+(tmax-tmin)*(double)l/(double)nfour;
2024               
2025           FORTRAN(tempload,(xforcold,xforc,&xforcact[l**nforc],iamforc,nforc,
2026             xloadold,xload,&xloadact[l*2**nload],iamload,nload,ibody,xbody,
2027             nbody,xbodyold,&xbodyact[l*7**nbody],t1old,t1,t1act,
2028             iamt1,nk,amta,namta,nam,ampli,&time,&reltime,ttime,&dtime,
2029             ithermal,nmethod,xbounold,xboun,&xbounact[l**nboun],iamboun,nboun,
2030             nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,co,vold,itg,&ntg,
2031             amname,ikboun,ilboun,nelemload,sideload,mi,
2032             ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
2033             iendset,ialset,&ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
2034           
2035       }
2036 
2037       SFREE(ampli);
2038 
2039       for(i=0;i<l**nboun;i++){xbounacttime[i]=xbounact[i];}
2040 
2041       /* determining the load frequency history:
2042          frequency transform of the load time history */
2043 
2044       for(i=0;i<*nforc;i++){
2045           for(l=0;l<nfour;l++){
2046               r[l]=xforcact[l**nforc+i];
2047           }
2048           FORTRAN(drffti,(&nfour,wsave,isave));
2049           FORTRAN(drfftf,(&nfour,r,wsave,isave));
2050           for(l=0;l<nfour;l++){
2051               xforcact[l**nforc+i]=r[l]/nfour*2.;
2052           }
2053           xforcact[i]=xforcact[i]/2.;
2054       }
2055 
2056       for(i=0;i<*nload;i++){
2057           for(l=0;l<nfour;l++){
2058               r[l]=xloadact[l*2**nload+2*i];
2059           }
2060           FORTRAN(drffti,(&nfour,wsave,isave));
2061           FORTRAN(drfftf,(&nfour,r,wsave,isave));
2062           for(l=0;l<nfour;l++){
2063               xloadact[l*2**nload+2*i]=r[l]/nfour*2.;
2064           }
2065           xloadact[2*i]=xloadact[2*i]/2.;
2066       }
2067 
2068       for(i=0;i<*nbody;i++){
2069           for(l=0;l<nfour;l++){
2070               r[l]=xbodyact[l**nbody+7*i];
2071           }
2072           FORTRAN(drffti,(&nfour,wsave,isave));
2073           FORTRAN(drfftf,(&nfour,r,wsave,isave));
2074           for(l=0;l<nfour;l++){
2075               xbodyact[l**nbody+7*i]=r[l]/nfour*2.;
2076           }
2077           xbodyact[7*i]=xbodyact[7*i]/2.;
2078       }
2079 
2080       if(iprescribedboundary){
2081           for(i=0;i<*nboun;i++){
2082               for(l=0;l<nfour;l++){
2083                   r[l]=xbounact[l**nboun+i];
2084               }
2085               FORTRAN(drffti,(&nfour,wsave,isave));
2086               FORTRAN(drfftf,(&nfour,r,wsave,isave));
2087               for(l=0;l<nfour;l++){
2088                   xbounact[l**nboun+i]=r[l]/nfour*2.;
2089               }
2090               xbounact[i]=xbounact[i]/2.;
2091           }
2092           
2093           /* LU decomposition of the stiffness matrix */
2094           
2095           if(*isolver==0){
2096 #ifdef SPOOLES
2097               spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
2098                         &symmetryflag,&inputformat,&nzs[2]);
2099 #else
2100               printf("*ERROR in steadystate: the SPOOLES library is not linked\n\n");
2101               FORTRAN(stop,());
2102 #endif
2103           }
2104           else if(*isolver==4){
2105 #ifdef SGI
2106               token=1;
2107               sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
2108 #else
2109               printf("*ERROR in steadystate: the SGI library is not linked\n\n");
2110               FORTRAN(stop,());
2111 #endif
2112           }
2113           else if(*isolver==5){
2114 #ifdef TAUCS
2115               tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
2116 #else
2117               printf("*ERROR in steadystate: the TAUCS library is not linked\n\n");
2118               FORTRAN(stop,());
2119 #endif
2120           }
2121           else if(*isolver==7){
2122 #ifdef PARDISO
2123               pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
2124                              &symmetryflag,&inputformat,jq,&nzs[2]);
2125 #else
2126               printf("*ERROR in steadystate: the PARDISO library is not linked\n\n");
2127               FORTRAN(stop,());
2128 #endif
2129           }
2130 
2131       }
2132 
2133       SFREE(r);SFREE(wsave);SFREE(isave);
2134 
2135       /* determining the frequency data points */
2136       
2137       NNEW(freqnh,double,ndata*(nev+1));
2138       
2139       ndatatot=0.;
2140       freqnh[0]=fmin;
2141       if(fabs(fmax-fmin)<1.e-10){
2142           ndatatot=1;
2143       }else{
2144 
2145       /* copy the eigenvalues and sort them in ascending order
2146          (important for values from distinct nodal diameters */
2147 
2148           NNEW(e,double,nev);
2149           for(i=0;i<nev;i++){e[i]=sqrt(d[i]);}
2150           FORTRAN(dsort,(e,&idummy,&nev,&iflag));
2151 
2152           for(i=0;i<nev;i++){
2153               if(i!=0){
2154                   if(fabs(e[i]-e[i-1])<1.e-5){continue;}
2155               }
2156               if(e[i]>=fmin){
2157                   if(e[i]<=fmax){
2158                       for(j=1;j<ndata;j++){
2159                           y=-1.+2.*j/((double)(ndata-1));
2160                           if(fabs(y)<1.e-10){freqnh[ndatatot+j]=
2161                                  (freqnh[ndatatot]+e[i])/2.;}
2162                           else{
2163                               freqnh[ndatatot+j]=(freqnh[ndatatot]+e[i])/2.+
2164                                   (e[i]-freqnh[ndatatot])*pow(fabs(y),1./bias)
2165                                   *y/(2.*fabs(y));
2166                           }
2167                       }
2168                       ndatatot+=ndata-1;
2169                   }
2170                   else{break;}
2171               }
2172           }
2173           SFREE(e);
2174           for(j=1;j<ndata;j++){
2175               y=-1.+2.*j/((double)(ndata-1));
2176               if(fabs(y)<1.e-10){freqnh[ndatatot+j]=
2177                          (freqnh[ndatatot]+fmax)/2.;}
2178               else{
2179                   freqnh[ndatatot+j]=(freqnh[ndatatot]+fmax)/2.+
2180                       (fmax-freqnh[ndatatot])*pow(fabs(y),1./bias)*
2181                       y/(2.*fabs(y));
2182               }
2183           }
2184           ndatatot+=ndata;
2185       }
2186       RENEW(freqnh,double,ndatatot);
2187 
2188       for(ii=0;ii<ndatatot;ii++){
2189           for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
2190 
2191           sprintf(description,"%12f",freqnh[ii]/(2.*pi));
2192           
2193           NNEW(xforcr,double,*nforc);
2194           NNEW(xloadr,double,2**nload);
2195           NNEW(xbodyr,double,7**nbody);
2196           for(k=0;k<7**nbody;k++){xbodyr[k]=xbody[k];}
2197           if(iprescribedboundary){
2198               NNEW(xbounr,double,*nboun);
2199               NNEW(fr,double,neq[1]); /* force corresponding to real particular solution */
2200               NNEW(ubr,double,neq[1]); /* real particular solution */
2201               NNEW(mubr,double,neq[1]); /* mass times real particular solution */
2202           }
2203           
2204           /* assigning the body forces to the elements */ 
2205 
2206           if(*nbody>0){
2207               ifreebody=*ne+1;
2208               NNEW(ipobody,ITG,2*ifreebody**nbody);
2209               for(k=1;k<=*nbody;k++){
2210                   FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
2211                             iendset,ialset,&inewton,nset,&ifreebody,&k));
2212                   RENEW(ipobody,ITG,2*(*ne+ifreebody));
2213               }
2214               RENEW(ipobody,ITG,2*(ifreebody-1));
2215           }
2216           
2217           NNEW(br,double,neq[1]); /* load rhs vector (real part) */
2218           NNEW(bi,double,neq[1]); /* load rhs vector (imaginary part) */
2219           NNEW(btot,double,nfour*neq[1]);
2220           NNEW(bp,double,nfour*neq[1]);
2221           
2222           NNEW(bjr,double,nev); /* real response modal decomposition */
2223           NNEW(bji,double,nev); /* imaginary response modal decomposition */
2224           
2225           NNEW(aa,double,nev); /* modal coefficients of the real loading */
2226           NNEW(bb,double,nev); /* modal coefficients of the imaginary loading */
2227           
2228           /* loop over all Fourier frequencies */
2229           
2230           NNEW(freq,double,nfour);
2231           
2232           for(l=0;l<nfour;l++){
2233               
2234               /* frequency */
2235               
2236               freq[l]=freqnh[ii]*floor((l+1.)/2.+0.1);
2237 
2238           /* calculating cc */
2239 
2240               if(dashpot){
2241                   NNEW(adc,double,neq[1]);
2242                   NNEW(auc,double,nzs[1]);
2243                   FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
2244                     ndirboun,xboun,nboun,
2245                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
2246                     nforc,nelemload,sideload,xload,nload,xbody,ipobody,
2247                     nbody,cgr,adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
2248                     ikmpc,ilmpc,ikboun,ilboun,
2249                     elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2250                     ielorien,norien,orab,ntmat_,
2251                     t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
2252                     nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
2253                     xstiff,npmat_,&dtime,matname,mi,ncmat_,
2254                     ttime,&freq[l],istep,&iinc,ibody,clearini,&mortar,springarea,
2255                     pslavsurf,pmastsurf,&reltime,&nasym));
2256                   
2257                   /*  zc = damping matrix * eigenmodes */
2258                   
2259                   NNEW(zc,double,neq[1]*nev);
2260                   for(i=0;i<nev;i++){
2261                       FORTRAN(op,(&neq[1],&z[(long long)i*neq[1]],&zc[i*neq[1]],
2262                                   adc,auc,jq,irow));
2263                   }
2264                   
2265                   /* cc is the reduced damping matrix (damping matrix mapped onto
2266                      space spanned by eigenmodes) */
2267                   
2268                   for(i=0;i<nev*nev;i++){cc[i]=0.;}
2269                   for(i=0;i<nev;i++){
2270                       for(j=0;j<=i;j++){
2271                           for(k=0;k<neq[1];k++){
2272                               cc[i*nev+j]+=z[(long long)j*neq[1]+k]*zc[i*neq[1]+k];
2273                           }
2274                       }
2275                   }
2276                   
2277                   /* symmetric part of cc matrix */
2278                   
2279                   for(i=0;i<nev;i++){
2280                       for(j=i;j<nev;j++){
2281                           cc[i*nev+j]=cc[j*nev+i];
2282                       }
2283                   }
2284                   SFREE(zc);SFREE(adc);SFREE(auc);
2285               }
2286               
2287               /* loading for this frequency */
2288               
2289               for(i=0;i<*nforc;i++){
2290                   xforcr[i]=xforcact[l**nforc+i];
2291               }
2292               
2293               for(i=0;i<*nload;i++){
2294                   xloadr[2*i]=xloadact[l*2**nload+2*i];
2295               }
2296               
2297               for(i=0;i<*nbody;i++){
2298                   xbodyr[7*i]=xbodyact[l**nbody+7*i];
2299               }
2300               
2301               /* calculating the instantaneous loading vector */
2302               
2303               FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
2304                 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcr,
2305                 nforc,nelemload,sideload,xloadr,nload,xbodyr,
2306                 ipobody,nbody,cgr,br,nactdof,&neq[1],nmethod,
2307                 ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
2308                 alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
2309                 t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
2310                 nplicon,plkcon,nplkcon,
2311                 npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
2312                 xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
2313                 ielprop,prop));
2314           
2315               /* real modal coefficients */
2316           
2317               if(!iprescribedboundary){
2318                   if(!cyclicsymmetry){
2319                       for(i=0;i<nev;i++){
2320                           i2=(long long)i*neq[1];
2321                           aa[i]=0.;
2322                           if(nactmech<neq[1]/2){
2323                               for(j=0;j<nactmech;j++){
2324                                   aa[i]+=z[i2+ikactmech[j]]*br[ikactmech[j]];
2325                               }
2326                           }else{
2327                               for(j=0;j<neq[1];j++){
2328                                   aa[i]+=z[i2+j]*br[j];
2329                               }
2330                           }
2331                       }
2332                   }else{
2333                       for(i=0;i<nev;i++){aa[i]=0.;}
2334                       for(j=0;j<nactmech;j++){
2335                           for(i=0;i<nev;i++){
2336                               FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
2337                               if(id!=0){
2338                                   if(izdof[id-1]==ikactmech[j]){
2339                                       aa[i]+=z[(long long)i*nzdof+id-1]*br[ikactmech[j]];
2340                                   }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2341                               }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2342                           }
2343                       }
2344                   }
2345           
2346                   /* imaginary modal coefficients */
2347                   
2348                   for(i=0;i<nev;i++){
2349                       bb[i]=0.;
2350                   }
2351 
2352               }else{
2353 
2354                   /* prescribed boundary conditions */
2355               
2356                   /* next statement makes sure that br is reset to zero at the
2357                      start of rhs.f */
2358 
2359                   nactmech=neq[1];
2360               
2361                   for(i=0;i<neq[1];i++){bi[i]=0.;}
2362                   
2363                   for(i=0;i<*nboun;i++){
2364                       xbounr[i]=xbounact[l**nboun+i];
2365                   }
2366                   
2367                   for(j=0;j<neq[1];j++){fr[j]=0.;ubr[j]=0.;}
2368                   for(i=0;i<*nboun;i++){
2369                       ic=neq[1]+i;
2370                       for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
2371                           ir=irow[j]-1;
2372                           fr[ir]=fr[ir]-au[j]*xbounr[i];
2373                           ubr[ir]=fr[ir];
2374                       }
2375                   }
2376                   if(*isolver==0){
2377 #ifdef SPOOLES
2378                       spooles_solve(ubr,&neq[1]);
2379 #endif
2380                   }
2381                   else if(*isolver==4){
2382 #ifdef SGI
2383                       sgi_solve(ubr,token);
2384 #endif
2385                   }
2386                   else if(*isolver==5){
2387 #ifdef TAUCS
2388                       tau_solve(ubr,&neq[1]);
2389 #endif
2390                   }
2391                   else if(*isolver==7){
2392 #ifdef PARDISO
2393                       pardiso_solve(ubr,&neq[1],&symmetryflag);
2394 #endif
2395                   }
2396                   FORTRAN(op,(&neq[1],ubr,mubr,adb,aub,jq,irow));
2397                   
2398                   for(i=0;i<neq[1];i++){
2399                       br[i]+=freq[l]*(freq[l]*mubr[i]);
2400                       bi[i]+=freq[l]*(-alpham*mubr[i]-betam*fr[i]);
2401                   }
2402               
2403                   /* real and imaginary modal coefficients */
2404                   
2405                   for(i=0;i<nev;i++){
2406                       aa[i]=0.;
2407                       for(j=0;j<neq[1];j++){
2408                           aa[i]+=z[(long long)i*neq[1]+j]*br[j];
2409                       }
2410                   }
2411                   
2412                   for(i=0;i<nev;i++){
2413                       bb[i]=0.;
2414                       for(j=0;j<neq[1];j++){
2415                           bb[i]+=z[(long long)i*neq[1]+j]*bi[j];
2416                       }
2417                   }
2418               }
2419               
2420               /* calculating the modal coefficients */
2421               
2422               if(dashpot==0){
2423                   for(i=0;i<nev;i++){
2424                       dd=pow(d[i]-pow(freq[l],2),2)+
2425                           pow(fric[i],2)*pow(freq[l],2);
2426                       bjr[i]=(aa[i]*(d[i]-freq[l]*freq[l])+
2427                               bb[i]*fric[i]*freq[l])/dd;
2428                       bji[i]=(bb[i]*(d[i]-freq[l]*freq[l])-
2429                               aa[i]*fric[i]*freq[l])/dd;
2430                   }
2431               }else{
2432                   nev2=2*nev;
2433                   NNEW(am,double,nev2*nev2);
2434                   NNEW(bm,double,nev2);
2435                   NNEW(ipiv,ITG,nev2);
2436                   
2437                   for(i=0;i<nev2;i++){
2438                       for(j=0;j<nev2;j++){
2439                           am[i*nev2+j]=0.;
2440                       }
2441                       bm[i]=0.;
2442                   }
2443                   for(i=0;i<nev;i++){
2444                       am[i*nev2+i]=d[i]-freq[l]*freq[l];
2445                       am[(i+nev)*nev2+i]=-fric[i]*freq[l];
2446                       bm[i]=aa[i];
2447                       am[i*nev2+nev+i]=-am[(i+nev)*nev2+i];
2448                       am[(i+nev)*nev2+nev+i]=am[i*nev2+i];
2449                       bm[nev+i]=bb[i];
2450                       for(j=0;j<nev;j++){
2451                           am[(j+nev)*nev2+i]=am[(j+nev)*nev2+i]
2452                               -cc[i*nev+j]*freq[l];
2453                           am[j*nev2+nev+i]=am[j*nev2+nev+i]
2454                               +cc[i*nev+j]*freq[l];
2455                       }
2456                   }
2457                   
2458                   /* solving the system of equations */
2459                   
2460                   FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
2461                   if(info!=0){
2462                       printf("*ERROR in steadystate: fatal termination of dgesv\n");
2463                       printf("       info=%" ITGFORMAT "\n",info);
2464 /*                FORTRAN(stop,());*/
2465                   }
2466                   
2467                   /* storing the solution in bjr and bji */
2468                   
2469                   for(i=0;i<nev;i++){
2470                       bjr[i]=bm[i];
2471                       bji[i]=bm[nev+i];
2472                   }
2473                   
2474                   SFREE(am);SFREE(bm);SFREE(ipiv);
2475               }
2476       
2477               /* calculating the real response */
2478               
2479               if(iprescribedboundary){
2480                   if(nmdnode==0){
2481                       memcpy(&br[0],&ubr[0],sizeof(double)*neq[1]);
2482                   }else{
2483                       for(i=0;i<nmddof;i++){
2484                           br[imddof[i]]=ubr[imddof[i]];
2485                       }
2486                   }
2487               }
2488               else{
2489                   if(nmdnode==0){
2490                       DMEMSET(br,0,neq[1],0.);
2491                   }else{
2492                       for(i=0;i<nmddof;i++){
2493                           br[imddof[i]]=0.;
2494                       }
2495                   }
2496               }
2497               
2498               if(!cyclicsymmetry){
2499                   if(nmdnode==0){
2500                       for(i=0;i<neq[1];i++){
2501                           for(j=0;j<nev;j++){
2502                               br[i]+=bjr[j]*z[(long long)j*neq[1]+i];
2503                           }
2504                       }
2505                   }else{
2506                       for(i=0;i<nmddof;i++){
2507                           for(j=0;j<nev;j++){
2508                               br[imddof[i]]+=bjr[j]*z[(long long)j*neq[1]+imddof[i]];
2509                           }
2510                       }
2511                   }
2512               }else{
2513                   for(i=0;i<nmddof;i++){
2514                       FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
2515                       if(id!=0){
2516                           if(izdof[id-1]==imddof[i]){
2517                               for(j=0;j<nev;j++){
2518                                   br[imddof[i]]+=bjr[j]*z[(long long)j*nzdof+id-1];
2519                               }
2520                           }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2521                       }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2522                   }
2523               }
2524           
2525               /* calculating the imaginary response */
2526               
2527               if(nmdnode==0){
2528                   DMEMSET(bi,0,neq[1],0.);
2529               }else{
2530                   for(i=0;i<nmddof;i++){
2531                       bi[imddof[i]]=0.;
2532                   }
2533               }
2534           
2535               if(!cyclicsymmetry){
2536                   if(nmdnode==0){
2537                       for(i=0;i<neq[1];i++){
2538                           for(j=0;j<nev;j++){
2539                               bi[i]+=bji[j]*z[(long long)j*neq[1]+i];
2540                           }
2541                       }
2542                   }else{
2543                       for(i=0;i<nmddof;i++){
2544                           for(j=0;j<nev;j++){
2545                               bi[imddof[i]]+=bji[j]*z[(long long)j*neq[1]+imddof[i]];
2546                           }
2547                       }
2548                   }
2549               }else{
2550                   for(i=0;i<nmddof;i++){
2551                       FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
2552                       if(id!=0){
2553                           if(izdof[id-1]==imddof[i]){
2554                               for(j=0;j<nev;j++){
2555                                   bi[imddof[i]]+=bji[j]*z[(long long)j*nzdof+id-1];
2556                               }
2557                           }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2558                       }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2559                   }
2560               }
2561               
2562               if(nmdnode==0){
2563               
2564                   /* magnitude and phase of the response */
2565 
2566                   for(i=0;i<neq[1];i++){
2567                       breal=br[i];
2568                       br[i]=sqrt(br[i]*br[i]+bi[i]*bi[i]);
2569                       if(fabs(breal)<1.e-10){
2570                           if(bi[i]>0.){bi[i]=pi/2.;}
2571                           else{bi[i]=-pi/2.;}
2572                       }
2573                       else{
2574                           bi[i]=atan(bi[i]/breal);
2575                           if(breal<0.){bi[i]+=pi;}
2576                       }
2577                   }
2578                   
2579                   /* correction for the sinus terms */
2580                   
2581                   if((l!=0)&&(2*(ITG)floor(l/2.+0.1)==l)){
2582                       for(i=0;i<neq[1];i++){
2583 //                        bi[i]-=pi/2.;}
2584                           bi[i]+=pi/2.;}
2585                   }
2586                   
2587                   /* contribution to the time response */
2588                   
2589                   for(j=0;j<nfour;j++){
2590                       time=tmin+2.*pi/freqnh[ii]*(double)j/(double)nfour;
2591                       for(i=0;i<neq[1];i++){
2592                           btot[j*neq[1]+i]+=br[i]*cos(freq[l]*time+bi[i]);
2593                           bp[j*neq[1]+i]-=freq[l]*br[i]*sin(freq[l]*time+bi[i]);
2594                       }
2595                   }
2596               }else{
2597               
2598                   /* magnitude and phase of the response */
2599 
2600                   for(jj=0;jj<nmddof;jj++){
2601                       i=imddof[jj];
2602                       breal=br[i];
2603                       br[i]=sqrt(br[i]*br[i]+bi[i]*bi[i]);
2604                       if(fabs(breal)<1.e-10){
2605                           if(bi[i]>0.){bi[i]=pi/2.;}
2606                           else{bi[i]=-pi/2.;}
2607                       }
2608                       else{
2609                           bi[i]=atan(bi[i]/breal);
2610                           if(breal<0.){bi[i]+=pi;}
2611                       }
2612                   }
2613                   
2614                   /* correction for the sinus terms */
2615                   
2616                   if((l!=0)&&(2*(ITG)floor(l/2.+0.1)==l)){
2617                       for(jj=0;jj<nmddof;jj++){
2618                           i=imddof[jj];
2619 //                        bi[i]-=pi/2.;}
2620                           bi[i]+=pi/2.;}
2621                   }
2622                   
2623                   /* contribution to the time response */
2624                   
2625                   for(j=0;j<nfour;j++){
2626                       time=tmin+2.*pi/freqnh[ii]*(double)j/(double)nfour;
2627                       for(jj=0;jj<nmddof;jj++){
2628                           i=imddof[jj];
2629                           btot[j*neq[1]+i]+=br[i]*cos(freq[l]*time+bi[i]);
2630                           bp[j*neq[1]+i]-=freq[l]*br[i]*sin(freq[l]*time+bi[i]);
2631                       }
2632                   }
2633               }
2634 
2635               /* resetting the part of br occupied by the variables to be printed
2636                  to zero */
2637 
2638               if(!iprescribedboundary){
2639                   if(nmdnode==0){
2640                       DMEMSET(br,0,neq[1],0.);
2641                   }else{
2642                       for(i=0;i<nmddof;i++){
2643                           br[imddof[i]]=0.;
2644                       }
2645                   }
2646               }
2647               
2648           }
2649           
2650           SFREE(xforcr);SFREE(xloadr);SFREE(xbodyr);SFREE(br);SFREE(bi);SFREE(freq);
2651           SFREE(bjr);SFREE(bji);SFREE(aa);SFREE(bb);
2652 
2653           if(*nbody>0) SFREE(ipobody);
2654           if(iprescribedboundary) {SFREE(xbounr);SFREE(fr);SFREE(ubr);SFREE(mubr);}
2655           
2656           
2657           /* result fields */
2658           
2659           NNEW(vr,double,mt**nk);
2660 
2661           if(intpointvar==1){
2662               NNEW(fn,double,mt**nk);
2663               NNEW(stn,double,6**nk);
2664               NNEW(stx,double,6*mi[0]**ne);
2665 
2666               if(*ithermal>1) {
2667                   NNEW(qfn,double,3**nk);
2668                   NNEW(qfx,double,3*mi[0]**ne);}
2669               
2670               if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
2671               if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2672               if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
2673               
2674               NNEW(eei,double,6*mi[0]**ne);
2675               if(*nener==1){
2676                   NNEW(stiini,double,6*mi[0]**ne);
2677                   NNEW(enerini,double,mi[0]**ne);}
2678           }
2679           
2680           /* storing the results */
2681           
2682           for(l=0;l<nfour;l++){
2683               time=tmin+2.*pi/freqnh[ii]*(double)l/(double)nfour;
2684               ptime=time;
2685               
2686               if(nmdnode==0){
2687                   DMEMSET(vr,0,mt**nk,0.);
2688               }else{
2689                   for(jj=0;jj<nmdnode;jj++){
2690                       i=imdnode[jj]-1;
2691                       for(j=1;j<4;j++){
2692                           vr[mt*i+j]=0.;
2693                       }
2694                   }
2695               }
2696       
2697               /* calculating displacements/temperatures */
2698 
2699               *nmethod=4;
2700               FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
2701                     ndirboun,&xbounacttime[l**nboun],nboun,ipompc,nodempc,
2702                     coefmpc,labmpc,nmpc,&btot[l*neq[1]],&bp[l*neq[1]],veold,&dtime,mi,
2703                     imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,nmethod,&time));
2704               *nmethod=5;
2705 
2706               results(co,nk,kon,ipkon,lakon,ne,vr,stn,inum,
2707                       stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
2708                       ielmat,ielorien,norien,orab,ntmat_,t0,t1,
2709                       ithermal,prestr,iprestr,filab,eme,emn,een,
2710                       iperturb,f,fn,nactdof,&iout,qa,
2711                       vold,&btot[l*neq[1]],nodeboun,ndirboun,
2712                       &xbounacttime[l**nboun],nboun,
2713                       ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
2714                       veold,accold,&bet,&gam,&dtime,&time,&xnull,
2715                       plicon,nplicon,plkcon,nplkcon,
2716                       xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
2717                       &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
2718                       enern,emeini,xstaten,eei,enerini,cocon,ncocon,
2719                       set,nset,istartset,iendset,ialset,nprint,prlab,prset,
2720                       qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,ikmpc,
2721                       ilmpc,istep,&iinc,springarea,&reltime,&ne0,xforc,nforc,
2722                       thicke,shcon,nshcon,
2723                       sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2724                       &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
2725                       islavsurf,ielprop,prop);
2726           
2727               (*kode)++;
2728               mode=-1;
2729               
2730               if(strcmp1(&filab[1044],"ZZS")==0){
2731                   NNEW(neigh,ITG,40**ne);
2732                   NNEW(ipneigh,ITG,*nk);
2733               }
2734 
2735               frd(co,&nkg,kon,ipkon,lakon,&neg,vr,stn,inum,nmethod,
2736                   kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2737                   nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2738                   ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2739                   mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&neg,
2740                   cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2741                   thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
2742 
2743               if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2744 
2745           }
2746           
2747           SFREE(vr);SFREE(btot);SFREE(bp);
2748 
2749           if(intpointvar==1){
2750               SFREE(fn);SFREE(stn);SFREE(stx);SFREE(eei);
2751               if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
2752               
2753               if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
2754               if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2755               if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
2756               
2757               if(*nener==1){SFREE(stiini);SFREE(enerini);}
2758           }
2759           
2760       }
2761       SFREE(xforcact);SFREE(xloadact);SFREE(xbodyact);SFREE(xbounact);
2762       SFREE(xbounacttime);SFREE(freqnh);
2763       if(*ithermal==1) SFREE(t1act);
2764       if(iprescribedboundary){
2765           if(*isolver==0){
2766 #ifdef SPOOLES
2767               spooles_cleanup();
2768 #endif
2769           }
2770           else if(*isolver==4){
2771 #ifdef SGI
2772               sgi_cleanup(token);
2773 #endif
2774           }
2775           else if(*isolver==5){
2776 #ifdef TAUCS
2777               tau_cleanup();
2778 #endif
2779           }
2780           else if(*isolver==7){
2781 #ifdef PARDISO
2782               pardiso_cleanup(&neq[1],&symmetryflag);
2783 #endif
2784           }
2785       }
2786 
2787       SFREE(ikactmech);
2788 
2789   }
2790 
2791   SFREE(adb);SFREE(aub);SFREE(z);SFREE(d);SFREE(inum);
2792 
2793   if(!cyclicsymmetry){
2794       SFREE(ad);SFREE(au);
2795   }else{
2796       SFREE(izdof);SFREE(nm);
2797 
2798       *nk/=nsectors;
2799       *ne/=nsectors;
2800       *nboun/=nsectors;
2801       neq[1]=neq[1]*2/nsectors;
2802 
2803       RENEW(co,double,3**nk);
2804       if(*ithermal!=0){
2805           RENEW(t0,double,*nk);
2806           RENEW(t1old,double,*nk);
2807           RENEW(t1,double,*nk);
2808           if(*nam>0) RENEW(iamt1,ITG,*nk);
2809       }
2810       RENEW(nactdof,ITG,mt**nk);
2811       if(*ntrans>0) RENEW(inotr,ITG,2**nk);
2812       RENEW(kon,ITG,*nkon);
2813       RENEW(ipkon,ITG,*ne);
2814       RENEW(lakon,char,8**ne);
2815       RENEW(ielmat,ITG,mi[2]**ne);
2816       if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne);
2817       RENEW(nodeboun,ITG,*nboun);
2818       RENEW(ndirboun,ITG,*nboun);
2819       if(*nam>0) RENEW(iamboun,ITG,*nboun);
2820       RENEW(xboun,double,*nboun);
2821       RENEW(xbounold,double,*nboun);
2822       RENEW(ikboun,ITG,*nboun);
2823       RENEW(ilboun,ITG,*nboun);
2824 
2825       /* recovering the original multiple point constraints */
2826 
2827       RENEW(ipompc,ITG,*nmpc);
2828       RENEW(nodempc,ITG,3**mpcend);
2829       RENEW(coefmpc,double,*mpcend);
2830       RENEW(labmpc,char,20**nmpc+1);
2831       RENEW(ikmpc,ITG,*nmpc);
2832       RENEW(ilmpc,ITG,*nmpc);
2833       RENEW(fmpc,double,*nmpc);
2834 
2835       *nmpc=nmpcold;
2836       *mpcend=mpcendold;
2837       for(i=0;i<*nmpc;i++){ipompc[i]=ipompcold[i];}
2838       for(i=0;i<3**mpcend;i++){nodempc[i]=nodempcold[i];}
2839       for(i=0;i<*mpcend;i++){coefmpc[i]=coefmpcold[i];}
2840       for(i=0;i<20**nmpc;i++){labmpc[i]=labmpcold[i];}
2841       for(i=0;i<*nmpc;i++){ikmpc[i]=ikmpcold[i];}
2842       for(i=0;i<*nmpc;i++){ilmpc[i]=ilmpcold[i];}
2843       SFREE(ipompcold);SFREE(nodempcold);SFREE(coefmpcold);
2844       SFREE(labmpcold);SFREE(ikmpcold);SFREE(ilmpcold);
2845 
2846       RENEW(vold,double,mt**nk);
2847       RENEW(veold,double,mt**nk);
2848       RENEW(eme,double,6*mi[0]**ne);
2849       if(*nener==1)RENEW(ener,double,mi[0]**ne);
2850 
2851 /* distributed loads */
2852 
2853       for(i=0;i<*nload;i++){
2854           if(nelemload[2*i]<=*ne*nsectors){
2855               nelemload[2*i]-=*ne*nelemload[2*i+1];
2856           }else{
2857               nelemload[2*i]-=*ne*(nsectors+nelemload[2*i+1]-1);
2858           }
2859       }
2860 
2861   /*  sorting the elements with distributed loads */
2862 
2863       if(*nload>0){
2864           if(*nam>0){
2865               FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
2866           }else{
2867               FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
2868           }
2869       }
2870       
2871 /* point loads */
2872       
2873       for(i=0;i<*nforc;i++){
2874           if(nodeforc[2*i]<=*nk*nsectors){
2875               nodeforc[2*i]-=*nk*nodeforc[2*i+1];
2876           }else{
2877               nodeforc[2*i]-=*nk*(nsectors+nodeforc[2*i+1]-1);
2878           }
2879       }
2880   }
2881 
2882   SFREE(xstiff);SFREE(fric);
2883 
2884   if(dashpot){SFREE(cc);}
2885 
2886   if(nherm!=1){SFREE(xmr);SFREE(xmi);}
2887 
2888   SFREE(imddof);SFREE(imdnode);SFREE(imdboun);SFREE(imdmpc);SFREE(imdelem);
2889 
2890   *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
2891   *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
2892   *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;*veoldp=veold;
2893   *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
2894   *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
2895   *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2896   *fmpcp=fmpc;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
2897 
2898 //  (*ttime)+=(*tper);
2899 
2900   return;
2901 }

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