root/src/dyna.c

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

DEFINITIONS

This source file includes following definitions.
  1. dyna

   1 /*     CalculiX - A 3-dimensional finite element program                   */
   2 /*              Copyright (C) 1998-2015 Guido Dhondt                          */
   3 /*     This program is free software; you can redistribute it and/or     */
   4 /*     modify it under the terms of the GNU General Public License as    */
   5 /*     published by the Free Software Foundation(version 2);    */
   6 /*                    */
   7 
   8 /*     This program is distributed in the hope that it will be useful,   */
   9 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  10 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  11 /*     GNU General Public License for more details.                      */
  12 
  13 /*     You should have received a copy of the GNU General Public License */
  14 /*     along with this program; if not, write to the Free Software       */
  15 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  16 
  17 #include <stdio.h>
  18 #include <math.h>
  19 #include <stdlib.h>
  20 #include <string.h>
  21 #include "CalculiX.h"
  22 
  23 #ifdef SPOOLES
  24    #include "spooles.h"
  25 #endif
  26 #ifdef SGI
  27    #include "sgi.h"
  28 #endif
  29 #ifdef TAUCS
  30    #include "tau.h"
  31 #endif
  32 #ifdef PARDISO
  33    #include "pardiso.h"
  34 #endif
  35 
  36 void dyna(double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, 
  37                ITG **nodebounp, ITG **ndirbounp, double **xbounp, ITG *nboun,
  38                ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp,
  39                ITG *nmpc, ITG *nodeforc,ITG *ndirforc,double *xforc, 
  40                ITG *nforc,ITG *nelemload, char *sideload,double *xload,
  41                ITG *nload, 
  42                ITG **nactdofp,ITG *neq, ITG *nzl,ITG *icol, ITG *irow, 
  43                ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG **ikbounp, 
  44                ITG **ilbounp,double *elcon, ITG *nelcon, double *rhcon, 
  45                ITG *nrhcon,double *cocon, ITG *ncocon,
  46                double *alcon, ITG *nalcon, double *alzero, 
  47                ITG **ielmatp,ITG **ielorienp, ITG *norien, double *orab, 
  48                ITG *ntmat_,double **t0p, 
  49                double **t1p,ITG *ithermal,double *prestr, ITG *iprestr, 
  50                double **voldp,ITG *iperturb, double **stip, ITG *nzs, 
  51                double *timepar, double *xmodal,
  52                double **veoldp, char *amname, double *amta,
  53                ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload,
  54                ITG **iamt1p,ITG *jout,
  55                ITG *kode, char *filab,double **emep, double *xforcold, 
  56                double *xloadold,
  57                double **t1oldp, ITG **iambounp, double **xbounoldp, ITG *iexpl,
  58                double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
  59                double **xstatep, ITG *npmat_, char *matname, ITG *mi,
  60                ITG *ncmat_, ITG *nstate_, double **enerp, char *jobnamec,
  61                double *ttime, char *set, ITG *nset, ITG *istartset,
  62                ITG *iendset, ITG **ialsetp, ITG *nprint, char *prlab,
  63                char *prset, ITG *nener, double *trab, 
  64                ITG **inotrp, ITG *ntrans, double **fmpcp, char *cbody, ITG *ibody,
  65                double *xbody, ITG *nbody, double *xbodyold, ITG *istep,
  66                ITG *isolver,ITG *jq, char *output, ITG *mcs, ITG *nkon,
  67                ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset,
  68                ITG *idrct, ITG *jmax,
  69                double *ctrl, ITG *itpamp, double *tietol,ITG *nalset,
  70                ITG *ikforc,ITG *ilforc,double *thicke,ITG *nslavs,ITG *nmat,
  71                char *typeboun,ITG *ielprop,double *prop){
  72 
  73   char fneig[132]="",description[13]="            ",*lakon=NULL,*labmpc=NULL,
  74     *labmpcold=NULL,lakonl[9]="        \0",*tchar1=NULL,*tchar2=NULL,
  75     *tchar3=NULL,cflag[1]=" ";
  76 
  77   ITG nev,i,j,k,idof,*inum=NULL,*ipobody=NULL,inewton=0,id,
  78     iinc=0,jprint=0,l,iout,ielas=0,icmd,iprescribedboundary,init,ifreebody,
  79     mode=-1,noddiam=-1,*kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,
  80     *inotr=NULL,*nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,
  81     *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
  82     *ilmpc=NULL,nsectors,nmpcold,mpcendold,*ipompcold=NULL,*nodempcold=NULL,
  83     *ikmpcold=NULL,*ilmpcold=NULL,kflag=2,nmd,nevd,*nm=NULL,*iamt1=NULL,
  84     *itg=NULL,ntg=0,symmetryflag=0,inputformat=0,dashpot,lrw,liw,iddebdf=0,
  85     *iwork=NULL,ngraph=1,nkg,neg,ncont,ne0,nkon0, *itietri=NULL,
  86     *koncont=NULL,konl[20],imat,nope,kodem,indexe,j1,jdof,
  87     *ipneigh=NULL,*neigh=NULL,inext,itp=0,*islavact=NULL,
  88     ismallsliding=0,isteadystate,mpcfree,im,cyclicsymmetry,
  89     memmpc_,imax,iener=0,*icole=NULL,*irowe=NULL,*jqe=NULL,nzse[3],
  90     nalset_=*nalset,*ialset=*ialsetp,*istartset_=NULL,*iendset_=NULL,
  91     *itiefac=NULL,*islavsurf=NULL,*islavnode=NULL,mt=mi[1]+1,
  92     *imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,mortar=0,*imastop=NULL,
  93     *iponoels=NULL,*inoels=NULL,*imddof=NULL,nmddof,
  94     *ikactcont=NULL,nactcont,nactcont_=100,*ikactmech=NULL,nactmech,
  95     iabsload=0,*ipe=NULL,*ime=NULL,iprev=1,inonlinmpc=0,ielem,
  96     *imdnode=NULL,nmdnode,*imdboun=NULL,nmdboun,*imdmpc=NULL,
  97     nmdmpc,intpointvar,kmin,kmax,i1,ifacecount,*izdof=NULL,
  98     nzdof,iload,iforc,*iponoel=NULL,*inoel=NULL,*imdelem=NULL,nmdelem,
  99     irenewxstate,nasym=0,*nshcon=NULL,nherm,icfd=0,*inomat=NULL;
 100 
 101   long long i2;
 102 
 103   double *d=NULL, *z=NULL, *b=NULL, *zeta=NULL,*stiini=NULL,
 104     *cd=NULL, *cv=NULL, *xforcact=NULL, *xloadact=NULL,*cc=NULL,
 105     *t1act=NULL, *ampli=NULL, *aa=NULL, *bb=NULL, *aanew=NULL, *bj=NULL, 
 106     *v=NULL,*aamech=NULL,*emn=NULL,*cdn=NULL,ptime,
 107     *stn=NULL, *stx=NULL, *een=NULL, *adb=NULL,*xstiff=NULL,*bjp=NULL,
 108     *aub=NULL, *temp_array1=NULL, *temp_array2=NULL, *aux=NULL,
 109     *f=NULL, *fn=NULL, *xbounact=NULL,*epn=NULL,*xstateini=NULL,
 110     *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
 111     *qfx=NULL, *xbodyact=NULL, *cgr=NULL, *au=NULL, *vbounact=NULL,
 112     *abounact=NULL,dtime,reltime,*t0=NULL,*t1=NULL,*t1old=NULL,
 113     physcon[1],zetaj,dj,ddj,h1,h2,h3,h4,h5,h6,sum,aai,bbi,tstart,tend,
 114     qa[3],cam[5],accold[1],bet,gam,*ad=NULL,sigma=0.,alpham,betam,
 115     *bact=NULL,*bmin=NULL,*co=NULL,*xboun=NULL,*xbounold=NULL,*vold=NULL,
 116     *eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,*coefmpcold,*veold=NULL,
 117     *xini=NULL,*rwork=NULL,*adc=NULL,*auc=NULL,*zc=NULL, *rpar=NULL,
 118     *cg=NULL,*straight=NULL,xl[27],voldl[mt*9],elas[21],fnl[27],t0l,t1l,
 119     elconloc[21],veoldl[mt*9],setnull,deltmx,bbmax,dd,dtheta,dthetaref,
 120     theta,*vini=NULL,dthetaold,*bcont=NULL,*vr=NULL,*vi=NULL,*bcontini=NULL,
 121     *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,precision,resultmaxprev,
 122     resultmax,func,funcp,fexp,fexm,fcos,fsin,sump,*bp=NULL,h14,senergy=0.0,
 123     *bv=NULL,*cstr=NULL,*aube=NULL,*adbe=NULL,*sti=*stip,time0=0.0,
 124     time=0.0,*xforcdiff=NULL,*xloaddiff=NULL,*xbodydiff=NULL,*t1diff=NULL,
 125     *xboundiff=NULL,*bprev=NULL,*bdiff=NULL,*areaslav=NULL,
 126     *springarea=NULL, *bold=NULL,*eenmax=NULL,*fnr=NULL,*fni=NULL,
 127     *xmastnor=NULL,*emeini=NULL,*xstate=NULL,*clearini=NULL,
 128     *shcon=NULL,*xmr=NULL,*xmi=NULL,*xnoels=NULL,*pslavsurf=NULL,
 129       *pmastsurf=NULL,*cdnr=NULL,*cdni=NULL,*tinc,*tper,*tmin,*tmax;
 130 
 131   FILE *f1;
 132 
 133   /* dummy variables for nonlinmpc */
 134 
 135   ITG *iaux=NULL,maxlenmpc,icascade=0,newstep=0,iit=-1,idiscon;
 136 
 137 #ifdef SGI
 138   ITG token;
 139 #endif
 140 
 141   /*    if iabsload=0: aamech is modified by the present incremental 
 142                        contribution of b
 143            iabsload=1: the last incremental contribution is
 144                        subtracted before the new one is added to b;
 145                        this latter incremental contribution is used
 146                        to update aamech
 147            iabsload=2: aamech is determined by the absolute
 148                        contribution of b (no incremental procedure
 149                        for the load; this is necessary if
 150                        - nonlinear MPC's are applied or
 151                        - user dloads are applied */
 152 
 153   co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
 154   ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
 155   ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;xstate=*xstatep;
 156   xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
 157   vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
 158   coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 159   fmpc=*fmpcp;veold=*veoldp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
 160 
 161   tinc=&timepar[0];
 162   tper=&timepar[1];
 163   tmin=&timepar[2];
 164   tmax=&timepar[3];
 165 
 166   if(ithermal[0]<=1){
 167       kmin=1;kmax=3;
 168   }else if(ithermal[0]==2){
 169       kmin=0;kmax=mi[1];if(kmax>2)kmax=2;
 170   }else{
 171       kmin=0;kmax=3;
 172   }
 173 
 174 //xstiffNNEW(//  xstiff,double,(long long)27*mi[0]**ne);
 175 
 176   dtime=*tinc;
 177 
 178   alpham=xmodal[0];
 179   betam=xmodal[1];
 180 
 181   dd=ctrl[16];deltmx=ctrl[26];
 182 
 183   /* determining nzl */
 184 
 185   *nzl=0;
 186   for(i=neq[1];i>0;i--){
 187       if(icol[i-1]>0){
 188           *nzl=i;
 189           break;
 190       }
 191   }
 192 
 193   /* opening the eigenvalue file and checking for cyclic symmetry */
 194 
 195   strcpy(fneig,jobnamec);
 196   strcat(fneig,".eig");
 197 
 198   if((f1=fopen(fneig,"rb"))==NULL){
 199     printf("*ERROR in dyna: cannot open eigenvalue file for reading");
 200     exit(0);
 201   }
 202 
 203   printf(" *INFO  in dyna: if there are problems reading the .eig file this may be due to:\n");
 204   printf("        1) the nonexistence of the .eig file\n");
 205   printf("        2) other boundary conditions than in the input deck\n");
 206   printf("           which created the .eig file\n\n");
 207 
 208   if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
 209       printf("*ERROR in dyna reading the cyclic symmetry flag in the eigenvalue file");
 210       exit(0);
 211   }
 212 
 213   if(fread(&nherm,sizeof(ITG),1,f1)!=1){
 214       printf("*ERROR in dyna reading the Hermitian flag in the eigenvalue file");
 215       exit(0);
 216   }
 217 
 218   if(nherm!=1){
 219       printf("*ERROR in dyna: the eigenvectors in the .eig-file result\n");
 220       printf("       from a non-Hermitian eigenvalue problem. The modal\n");
 221       printf("       dynamic procedure cannot handle that yet\n\n");
 222       FORTRAN(stop,());
 223   }
 224 
 225   /* creating imddof containing the degrees of freedom
 226      retained by the user and imdnode containing the nodes */
 227 
 228   nmddof=0;nmdnode=0;nmdboun=0;nmdmpc=0;nmdelem=0;
 229 
 230   NNEW(imddof,ITG,*nk*3);
 231   NNEW(imdnode,ITG,*nk);
 232   NNEW(imdboun,ITG,*nboun);
 233   NNEW(imdmpc,ITG,*nmpc);
 234   FORTRAN(createmddof,(imddof,&nmddof,istartset,iendset,
 235                        ialset,nactdof,ithermal,mi,imdnode,&nmdnode,
 236                        ikmpc,ilmpc,ipompc,nodempc,nmpc,
 237                        imdmpc,&nmdmpc,imdboun,&nmdboun,ikboun,
 238                        nboun,nset,ntie,tieset,set,lakon,kon,ipkon,labmpc,
 239                        ilboun,filab,prlab,prset,nprint,ne,&cyclicsymmetry));
 240 
 241   /* if results are requested in too many nodes, it is faster to 
 242      calculate the results in all nodes */
 243 
 244   if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
 245       nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 246   }
 247 
 248   if(nmdnode!=0){
 249       if(!cyclicsymmetry){
 250           for(i=0;i<*nload;i++){
 251               iload=i+1;
 252               FORTRAN(addimdnodedload,(nelemload,sideload,ipkon,kon,lakon,
 253                       &iload,imdnode,&nmdnode,ikmpc,ilmpc,ipompc,nodempc,nmpc,
 254                       imddof,&nmddof,nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 255                       ikboun,nboun,ilboun,ithermal));
 256           }
 257           for(i=0;i<*nforc;i++){
 258               iforc=i+1;
 259               FORTRAN(addimdnodecload,(nodeforc,&iforc,imdnode,&nmdnode,xforc,
 260                    ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
 261                    nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 262                    ikboun,nboun,ilboun,ithermal));
 263           }
 264       }
 265           
 266       /* determining the elements belonging to a given node */
 267       
 268       NNEW(iponoel,ITG,*nk);
 269       NNEW(inoel,ITG,2**nkon);
 270       FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
 271       NNEW(imdelem,ITG,*ne);
 272 
 273       /* storing the elements in which integration point results
 274          are needed; storing the nodes which are needed to
 275          calculate these results */
 276 
 277       FORTRAN(createmdelem,(imdnode,&nmdnode,xforc,
 278                    ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
 279                    nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
 280                    ikboun,nboun,ilboun,ithermal,imdelem,&nmdelem,
 281                    iponoel,inoel,prlab,prset,nprint,lakon,set,nset,
 282                    ialset,ipkon,kon,istartset,iendset,nforc,
 283                    ikforc,ilforc));
 284 
 285       RENEW(imdelem,ITG,nmdelem);
 286       SFREE(iponoel);SFREE(inoel);
 287   }
 288 
 289   /* if results are requested in too many nodes, it is faster to 
 290      calculate the results in all nodes */
 291 
 292   if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
 293       nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 294   }
 295   
 296   /* subtracting 1 to comply with the C-convention */
 297 
 298   for(i=0;i<nmddof;i++){imddof[i]-=1;}
 299   RENEW(imddof,ITG,nmddof);
 300   RENEW(imdnode,ITG,nmdnode);
 301   RENEW(imdboun,ITG,nmdboun);
 302   RENEW(imdmpc,ITG,nmdmpc);
 303 
 304   nsectors=1;
 305 
 306   /* reading the eigenvalues / eigenmodes */
 307 
 308   if(!cyclicsymmetry){
 309 
 310       nkg=*nk;
 311       neg=*ne;
 312 
 313       if(fread(&nev,sizeof(ITG),1,f1)!=1){
 314           printf("*ERROR in dyna reading the number of eigenvalues in the eigenvalue file");
 315           exit(0);
 316       }
 317       
 318       NNEW(d,double,nev);
 319       
 320       if(fread(d,sizeof(double),nev,f1)!=nev){
 321           printf("*ERROR in dyna reading the eigenvalues in the eigenvalue file");
 322           exit(0);
 323       }
 324 
 325       for(i=0;i<nev;i++){
 326           if(d[i]>0){d[i]=sqrt(d[i]);}else{d[i]=0.;}
 327       }
 328       
 329       NNEW(ad,double,neq[1]);
 330       NNEW(adb,double,neq[1]);
 331       NNEW(au,double,nzs[2]);
 332       NNEW(aub,double,nzs[1]);
 333       
 334       if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
 335           printf("*ERROR in dyna reading the diagonal of the stiffness matrix in the eigenvalue file");
 336           exit(0);
 337       }
 338       
 339       if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
 340           printf("*ERROR in dyna reading the off-diagonals of the stiffness matrix in the eigenvalue file");
 341           exit(0);
 342       }
 343       
 344       if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 345           printf("*ERROR in dyna reading the diagonal of the mass matrix in the eigenvalue file");
 346           exit(0);
 347       }
 348       
 349       if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 350           printf("*ERROR in dyna reading the off-diagonals of the mass matrix in the  eigenvalue file");
 351           exit(0);
 352       }
 353       
 354       NNEW(z,double,neq[1]*nev);
 355       
 356       if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
 357           printf("*ERROR in dyna reading the eigenvectors in the eigenvalue file");
 358           exit(0);
 359       }
 360   }
 361   else{
 362       nev=0;
 363       do{
 364           if(fread(&nmd,sizeof(ITG),1,f1)!=1){
 365               break;
 366           }
 367           if(fread(&nevd,sizeof(ITG),1,f1)!=1){
 368               printf("*ERROR in dyna reading the number of eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 369               exit(0);
 370           }
 371           if(nev==0){
 372               NNEW(d,double,nevd);
 373               NNEW(nm,ITG,nevd);
 374           }else{
 375               RENEW(d,double,nev+nevd);
 376               RENEW(nm,ITG,nev+nevd);
 377           }
 378           
 379           if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
 380               printf("*ERROR in dyna reading the eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 381               exit(0);
 382           }
 383 
 384           for(i=nev;i<nev+nevd;i++){
 385               if(d[i]>0){d[i]=sqrt(d[i]);}else{d[i]=0.;}
 386           }
 387 
 388           for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
 389           
 390           if(nev==0){
 391               NNEW(adb,double,neq[1]);
 392               NNEW(aub,double,nzs[1]);
 393 
 394               if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
 395                   printf("*ERROR in dyna reading the diagonal of the mass matrix in the eigenvalue file");
 396                   exit(0);
 397               }
 398               
 399               if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 400                   printf("*ERROR in dyna reading the off-diagonals of the mass matrix in the eigenvalue file");
 401                   exit(0);
 402               }
 403           }
 404           
 405           if(nev==0){
 406               NNEW(z,double,neq[1]*nevd);
 407           }else{
 408               RENEW(z,double,(long long)neq[1]*(nev+nevd));
 409           }
 410           
 411           if(fread(&z[(long long)neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
 412               printf("*ERROR in dyna reading the eigenvectors for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
 413               exit(0);
 414           }
 415           nev+=nevd;
 416       }while(1);
 417 
 418       /* determining the maximum amount of segments */
 419 
 420       for(i=0;i<*mcs;i++){
 421 //        if(cs[17*i]>nsectors) nsectors=cs[17*i];
 422           if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
 423       }
 424 
 425         /* determining the maximum number of sectors to be plotted */
 426 
 427       for(j=0;j<*mcs;j++){
 428           if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
 429       }
 430       nkg=*nk*ngraph;
 431       neg=*ne*ngraph;
 432 
 433       /* allocating field for the expanded structure */
 434 
 435       RENEW(co,double,3**nk*nsectors);
 436 
 437       /* next line is necessary for multiple cyclic symmetry
 438          conditions */
 439 
 440       for(i=3**nk;i<3**nk*nsectors;i++){co[i]=0.;}
 441       if(*ithermal!=0){
 442           RENEW(t0,double,*nk*nsectors);
 443           RENEW(t1old,double,*nk*nsectors);
 444           RENEW(t1,double,*nk*nsectors);
 445           if(*nam>0) RENEW(iamt1,ITG,*nk*nsectors);
 446       }
 447       RENEW(nactdof,ITG,mt**nk*nsectors);
 448       if(*ntrans>0) RENEW(inotr,ITG,2**nk*nsectors);
 449       RENEW(kon,ITG,*nkon*nsectors);
 450       RENEW(ipkon,ITG,*ne*nsectors);
 451       for(i=*ne;i<*ne*nsectors;i++){ipkon[i]=-1;}
 452       RENEW(lakon,char,8**ne*nsectors);
 453       RENEW(ielmat,ITG,mi[2]**ne*nsectors);
 454       if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne*nsectors);
 455 //      RENEW(z,double,(long long)neq[1]*nev*nsectors/2);
 456 
 457       RENEW(nodeboun,ITG,*nboun*nsectors);
 458       RENEW(ndirboun,ITG,*nboun*nsectors);
 459       if(*nam>0) RENEW(iamboun,ITG,*nboun*nsectors);
 460       RENEW(xboun,double,*nboun*nsectors);
 461       RENEW(xbounold,double,*nboun*nsectors);
 462       RENEW(ikboun,ITG,*nboun*nsectors);
 463       RENEW(ilboun,ITG,*nboun*nsectors);
 464 
 465       NNEW(ipompcold,ITG,*nmpc);
 466       NNEW(nodempcold,ITG,3**mpcend);
 467       NNEW(coefmpcold,double,*mpcend);
 468       NNEW(labmpcold,char,20**nmpc);
 469       NNEW(ikmpcold,ITG,*nmpc);
 470       NNEW(ilmpcold,ITG,*nmpc);
 471 
 472       for(i=0;i<*nmpc;i++){ipompcold[i]=ipompc[i];}
 473       for(i=0;i<3**mpcend;i++){nodempcold[i]=nodempc[i];}
 474       for(i=0;i<*mpcend;i++){coefmpcold[i]=coefmpc[i];}
 475       for(i=0;i<20**nmpc;i++){labmpcold[i]=labmpc[i];}
 476       for(i=0;i<*nmpc;i++){ikmpcold[i]=ikmpc[i];}
 477       for(i=0;i<*nmpc;i++){ilmpcold[i]=ilmpc[i];}
 478       nmpcold=*nmpc;
 479       mpcendold=*mpcend;
 480 
 481       RENEW(ipompc,ITG,*nmpc*nsectors);
 482       RENEW(nodempc,ITG,3**mpcend*nsectors);
 483       RENEW(coefmpc,double,*mpcend*nsectors);
 484       RENEW(labmpc,char,20**nmpc*nsectors+1);
 485       RENEW(ikmpc,ITG,*nmpc*nsectors);
 486       RENEW(ilmpc,ITG,*nmpc*nsectors);
 487       RENEW(fmpc,double,*nmpc*nsectors);
 488 
 489       /* determining the space needed to expand the
 490          contact surfaces */
 491 
 492       NNEW(tchar1,char,81);
 493       NNEW(tchar2,char,81);
 494       NNEW(tchar3,char,81);
 495       for(i=0; i<*ntie; i++){
 496         if(tieset[i*(81*3)+80]=='C'){
 497 
 498           //a contact constraint was found, so increase nalset
 499 
 500           memcpy(tchar2,&tieset[i*(81*3)+81],81);
 501           tchar2[80]='\0';
 502           memcpy(tchar3,&tieset[i*(81*3)+81+81],81);
 503           tchar3[80]='\0';
 504           for(j=0; j<*nset; j++){
 505             memcpy(tchar1,&set[j*81],81);
 506             tchar1[80]='\0';
 507             if(strcmp(tchar1,tchar2)==0){
 508 
 509               //dependent nodal surface was found
 510 
 511               (*nalset)+=(iendset[j]-istartset[j]+1)*(nsectors);
 512             }
 513             else if(strcmp(tchar1,tchar3)==0){
 514 
 515               //independent element face surface was found
 516 
 517               (*nalset)+=(iendset[j]-istartset[j]+1)*(nsectors);
 518             }
 519           }
 520         }
 521       }
 522       SFREE(tchar1);
 523       SFREE(tchar2);
 524       SFREE(tchar3);
 525 
 526       RENEW(ialset,ITG,*nalset);
 527 
 528       /* save the information in istarset and isendset */
 529 
 530       NNEW(istartset_,ITG,*nset);
 531       NNEW(iendset_,ITG,*nset);
 532       for(j=0; j<*nset; j++){
 533         istartset_[j]=istartset[j];
 534         iendset_[j]=iendset[j];
 535       }
 536 
 537       /* reallocating the fields for the nodes in which the
 538          solution has to be calculated */
 539 
 540       RENEW(imddof,ITG,neq[1]/2*nsectors);
 541       RENEW(imdnode,ITG,*nk*nsectors);
 542       RENEW(imdboun,ITG,*nboun*nsectors);
 543       RENEW(imdmpc,ITG,*nmpc*nsectors);
 544 
 545 //izdofNNEW(//      izdof,ITG,1);
 546       
 547       expand(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 548         ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
 549         nforc,nelemload,sideload,xload,nload,nactdof,neq,
 550         nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
 551         alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
 552         t0,ithermal,prestr,iprestr,vold,iperturb,sti,nzs,
 553         adb,aub,filab,eme,plicon,nplicon,plkcon,nplkcon,
 554         xstate,npmat_,matname,mi,ics,cs,mpcend,ncmat_,
 555         nstate_,mcs,nkon,ener,jobnamec,output,set,nset,istartset,
 556         iendset,ialset,nprint,prlab,prset,nener,trab,
 557         inotr,ntrans,ttime,fmpc,&nev,&z,iamboun,xbounold,
 558         &nsectors,nm,icol,irow,nzl,nam,ipompcold,nodempcold,coefmpcold,
 559         labmpcold,&nmpcold,xloadold,iamload,t1old,t1,iamt1,xstiff,
 560         &icole,&jqe,&irowe,isolver,nzse,&adbe,&aube,iexpl,
 561         ibody,xbody,nbody,cocon,ncocon,tieset,ntie,imddof,&nmddof,
 562         imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,&izdof,&nzdof,
 563         &nherm,xmr,xmi,typeboun,ielprop,prop);
 564 
 565       RENEW(imddof,ITG,nmddof);
 566       RENEW(imdnode,ITG,nmdnode);
 567       RENEW(imdboun,ITG,nmdboun);
 568       RENEW(imdmpc,ITG,nmdmpc);
 569 
 570       SFREE(vold);
 571       NNEW(vold,double,mt**nk);
 572       SFREE(veold);
 573       NNEW(veold,double,mt**nk);
 574       RENEW(eme,double,6*mi[0]**ne);
 575       RENEW(sti,double,6*mi[0]**ne);
 576 
 577 //      RENEW(xstiff,double,(long long)27*mi[0]**ne);
 578       if(*nener==1) RENEW(ener,double,mi[0]**ne*2);
 579   }
 580 
 581   fclose(f1);
 582           
 583   /* checking for steadystate calculations */
 584 
 585   if(*tper<0){
 586       precision=-*tper;
 587       *tper=1.e10;
 588       isteadystate=1;
 589   }else{
 590       isteadystate=0;
 591   }
 592 
 593   /* checking for nonlinear MPC's */
 594 
 595   for(i=0;i<*nmpc;i++){
 596       if((strcmp1(&labmpc[20*i]," ")!=0)&&
 597          (strcmp1(&labmpc[20*i],"CONTACT")!=0)&&
 598          (strcmp1(&labmpc[20*i],"CYCLIC")!=0)&&
 599          (strcmp1(&labmpc[20*i],"SUBCYCLIC")!=0)){
 600           inonlinmpc=1;
 601           iabsload=2;
 602           break;
 603       }
 604   }
 605 
 606 
 607   /* normalizing the time */
 608 
 609   FORTRAN(checktime,(itpamp,namta,tinc,ttime,amta,tmin,&inext,&itp,istep));
 610   dtheta=(*tinc)/(*tper);
 611   dthetaref=dtheta;
 612   dthetaold=dtheta;
 613 
 614   *tmin=*tmin/(*tper);
 615   *tmax=*tmax/(*tper);
 616   theta=0.;
 617 
 618   /* check for rigid body modes 
 619      if there is a jump of 1.e4 in two subsequent eigenvalues
 620      all eigenvalues preceding the jump are considered to
 621      be rigid body modes and their frequency is set to zero */
 622 
 623   setnull=1.;
 624   for(i=nev-2;i>-1;i--){
 625       if(fabs(d[i])<0.0001*fabs(d[i+1])) setnull=0.;
 626       d[i]*=setnull;
 627   }
 628 
 629   /* check whether there are dashpot elements */
 630 
 631   dashpot=0;
 632   for(i=0;i<*ne;i++){
 633       if(ipkon[i]<0) continue;
 634       if(strcmp1(&lakon[i*8],"ED")==0){
 635           dashpot=1;break;}
 636   }
 637 
 638   if(dashpot){
 639 
 640       if(cyclicsymmetry){
 641           printf("*ERROR in dyna: dashpots are not allowed in combination with cyclic symmetry\n");
 642           FORTRAN(stop,());
 643       }
 644 
 645       liw=51;
 646       NNEW(iwork,ITG,liw);
 647       lrw=130+42*nev;
 648       NNEW(rwork,double,lrw);
 649       NNEW(xini,double,2*nev);
 650       NNEW(adc,double,neq[1]);
 651       NNEW(auc,double,nzs[1]);
 652       FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 653               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 654               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 655               adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
 656               ikmpc,ilmpc,ikboun,ilboun,
 657               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 658               ielorien,norien,orab,ntmat_,
 659               t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 660               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 661               xstiff,npmat_,&dtime,matname,mi,ncmat_,
 662               ttime,&time0,istep,&iinc,ibody,clearini,&mortar,springarea,
 663               pslavsurf,pmastsurf,&reltime,&nasym));
 664 
 665       /*  zc = damping matrix * eigenmodes */
 666 
 667       NNEW(zc,double,neq[1]*nev);
 668       for(i=0;i<nev;i++){
 669           FORTRAN(op,(&neq[1],&z[i*neq[1]],&zc[i*neq[1]],adc,auc,
 670           jq,irow));
 671       }
 672 
 673       /* cc is the reduced damping matrix (damping matrix mapped onto
 674          space spanned by eigenmodes) */
 675 
 676       NNEW(cc,double,nev*nev);
 677       for(i=0;i<nev;i++){
 678           for(j=0;j<=i;j++){
 679               for(k=0;k<neq[1];k++){
 680                   cc[i*nev+j]+=z[j*neq[1]+k]*zc[i*neq[1]+k];
 681               }
 682           }
 683       }
 684 
 685       /* symmetric part of cc matrix */
 686 
 687       for(i=0;i<nev;i++){
 688           for(j=i;j<nev;j++){
 689               cc[i*nev+j]=cc[j*nev+i];
 690           }
 691       }
 692       SFREE(zc);
 693   }
 694 
 695   /* contact conditions */
 696 
 697   if(*nslavs==0){irenewxstate=1;}else{irenewxstate=0;}
 698   inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
 699           lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
 700           &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
 701           &mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,&ifacecount,
 702           nmpc,&mpcfree,&memmpc_,
 703           &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 704           iperturb,ikboun,nboun,co,istep,&xnoels);
 705 
 706   if(ncont!=0){
 707 
 708       if(mortar>0){
 709           printf("*ERROR in dyna: modal dynamics cannot be combined with\n");
 710           printf("       face-to-face penalty contact\n\n");
 711           FORTRAN(stop,());
 712       }
 713           
 714       if(dashpot){
 715           printf("*ERROR in dyna: contact is not allowed in combination with dashpots\n");
 716           FORTRAN(stop,());
 717       }
 718       RENEW(ipkon,ITG,*ne+*nslavs);
 719       RENEW(lakon,char,8*(*ne+*nslavs));
 720       if(*nener==1){
 721         RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
 722       }
 723 
 724       /* 11 instead of 10: last position is reserved for how
 725          many dependent nodes are paired to this face */
 726     
 727       RENEW(kon,ITG,*nkon+11**nslavs);
 728       if(*norien>0){
 729           RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
 730           for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
 731       }
 732       RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
 733       for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
 734       NNEW(cg,double,3*ncont);
 735       NNEW(straight,double,16*ncont);
 736 
 737       /* internal state variables for contact */
 738 
 739       if((irenewxstate==1)&&(*nslavs!=0)&&(*nstate_>0)){
 740           RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
 741           for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
 742               xstate[k]=0.;
 743           }
 744       }
 745       if(*nstate_>0){
 746           NNEW(xstateini,double,*nstate_*mi[0]*(*ne+*nslavs));
 747           for(k=0;k<*nstate_*mi[0]*(*ne+*nslavs);++k){
 748               xstateini[k]=xstate[k];
 749           }
 750       }
 751 
 752       NNEW(xmastnor,double,3*nmastnode[*ntie]);
 753       NNEW(areaslav,double,ifacecount);
 754       NNEW(springarea,double,2**nslavs);
 755       NNEW(vini,double,mt**nk);
 756       NNEW(bcontini,double,neq[1]);
 757       NNEW(bcont,double,neq[1]);
 758       NNEW(ikactcont,ITG,nactcont_);
 759   }
 760 
 761   /* storing the element and topology information before introducing 
 762      contact elements */
 763 
 764   ne0=*ne;nkon0=*nkon;
 765 
 766   NNEW(zeta,double,nev);
 767   NNEW(cstr,double,6);
 768 
 769   /* calculating the damping coefficients*/
 770 
 771   if(xmodal[10]<0){
 772       for(i=0;i<nev;i++){
 773         if(fabs(d[i])>(1.e-10)){
 774           zeta[i]=(alpham+betam*d[i]*d[i])/(2.*d[i]);
 775         }
 776         else {
 777           printf("*WARNING in dyna: one of the frequencies is zero\n");
 778           printf("         no Rayleigh mass damping allowed\n");
 779           zeta[i]=0.;
 780         }
 781       }
 782   }
 783   else{
 784 
 785     /*copy the damping coefficients for every eigenfrequency from xmodal[11....] */
 786 
 787     if(nev<(ITG)xmodal[10]){
 788       imax=nev;
 789       printf("*WARNING in dyna: too many modal damping coefficients applied\n");
 790       printf("         damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
 791     }
 792     else{
 793       imax=(ITG)xmodal[10];
 794     }
 795     for(i=0; i<imax; i++){
 796       zeta[i]=xmodal[11+i];     
 797     }
 798     
 799   }
 800 
 801   /* modal decomposition of the initial conditions */
 802   /* for cyclic symmetric structures the initial conditions
 803      are assumed to be zero */
 804   
 805   NNEW(cd,double,nev);
 806   NNEW(cv,double,nev);
 807 
 808   if(!cyclicsymmetry){
 809     NNEW(temp_array1,double,neq[1]);
 810     NNEW(temp_array2,double,neq[1]);
 811     for(i=0;i<neq[1];i++){temp_array1[i]=0;temp_array2[i]=0;}
 812     
 813     /* displacement initial conditions */
 814     
 815     for(i=0;i<*nk;i++){
 816       for(j=0;j<mt;j++){
 817         if(nactdof[mt*i+j]!=0){
 818           idof=nactdof[mt*i+j]-1;
 819           temp_array1[idof]=vold[mt*i+j];
 820         }
 821       }
 822     }
 823     
 824     FORTRAN(op,(&neq[1],temp_array1,temp_array2,adb,aub,jq,irow));
 825     
 826     for(i=0;i<neq[1];i++){
 827       for(k=0;k<nev;k++){
 828         cd[k]+=z[k*neq[1]+i]*temp_array2[i];
 829       }
 830     }
 831     
 832     /* velocity initial conditions */
 833     
 834     for(i=0;i<neq[1];i++){temp_array1[i]=0;temp_array2[i]=0;}
 835     for(i=0;i<*nk;i++){
 836       for(j=0;j<mt;j++){
 837         if(nactdof[mt*i+j]!=0){
 838           idof=nactdof[mt*i+j]-1;
 839           temp_array1[idof]=veold[mt*i+j];
 840         }
 841       }
 842     }
 843     
 844     FORTRAN(op,(&neq[1],temp_array1,temp_array2,adb,aub,jq,irow));
 845     
 846     for(i=0;i<neq[1];i++){
 847       for(k=0;k<nev;k++){
 848         cv[k]+=z[k*neq[1]+i]*temp_array2[i];
 849       }
 850     }
 851     
 852     SFREE(temp_array1);SFREE(temp_array2);
 853                       
 854   }
 855   NNEW(xforcact,double,*nforc);
 856   NNEW(xforcdiff,double,*nforc);
 857   NNEW(xloadact,double,2**nload);
 858   NNEW(xloaddiff,double,2**nload);
 859   NNEW(xbodyact,double,7**nbody);
 860   NNEW(xbodydiff,double,7**nbody);
 861 
 862   /* copying the rotation axis and/or acceleration vector */
 863 
 864   for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
 865   for(k=0;k<7**nbody;k++){xbodydiff[k]=xbody[k];}
 866   NNEW(xbounact,double,*nboun);
 867   NNEW(xboundiff,double,*nboun);
 868   if(*ithermal==1) {NNEW(t1act,double,*nk);
 869       NNEW(t1diff,double,*nk);}
 870   
 871   /* assigning the body forces to the elements */ 
 872 
 873   if(*nbody>0){
 874       ifreebody=*ne+1;
 875       NNEW(ipobody,ITG,2**ne);
 876       for(k=1;k<=*nbody;k++){
 877           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 878                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 879           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 880       }
 881       RENEW(ipobody,ITG,2*(ifreebody-1));
 882   }
 883                
 884   NNEW(b,double,neq[1]); /* load rhs vector and displacement solution vector */
 885   NNEW(bp,double,neq[1]); /* velocity solution vector */
 886   NNEW(bj,double,nev); /* response modal decomposition */
 887   NNEW(bjp,double,nev); /* derivative of the response modal decomposition */
 888   NNEW(ampli,double,*nam); /* instantaneous amplitude */
 889 
 890   /* constant coefficient of the linear amplitude function */
 891 
 892   NNEW(aa,double,nev); 
 893   NNEW(aanew,double,nev);
 894   NNEW(aamech,double,nev);
 895 
 896   /* linear coefficient of the linear amplitude function */
 897 
 898   NNEW(bb,double,nev);
 899   
 900   NNEW(v,double,mt**nk);
 901   NNEW(fn,double,mt**nk);
 902   NNEW(stn,double,6**nk);
 903   NNEW(inum,ITG,*nk);
 904   strcpy1(&cflag[0],&filab[4],1);
 905   FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
 906           nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat));
 907 
 908   if(*ithermal>1) {NNEW(qfn,double,3**nk);
 909       NNEW(qfx,double,3*mi[0]**ne);}
 910 
 911   if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 912   if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 913   if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
 914   if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 915 
 916   NNEW(eei,double,6*mi[0]**ne);
 917   if(*nener==1){
 918     NNEW(stiini,double,6*mi[0]**ne);
 919     NNEW(enerini,double,mi[0]**ne);}
 920 
 921   /*  check for nonzero SPC's */
 922 
 923   iprescribedboundary=0;
 924   for(i=0;i<*nboun;i++){
 925       if(fabs(xboun[i])>1.e-10){
 926           iprescribedboundary=1;
 927           nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
 928           break;
 929       }
 930   }
 931 
 932 /*     calculating the instantaneous loads (forces, surface loading, 
 933        centrifugal and gravity loading or temperature) at time 0 
 934        setting iabsload to 2 if user subroutine dload is used */
 935 
 936 /*  FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,
 937      xload,xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,
 938      xbodyact,t1old,t1,t1act,iamt1,nk,
 939      amta,namta,nam,ampli,&time0,&reltime,ttime,&dtime,ithermal,nmethod,
 940      xbounold,xboun,xbounact,iamboun,nboun,
 941      nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 942      co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi));*/
 943       
 944   FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
 945               xloadold,xload,xloadact,iamload,nload,ibody,xbody,
 946               nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
 947               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
 948               nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
 949               ndirboun,nodeforc,
 950               ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
 951               nelemload,sideload,mi,
 952               xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
 953               &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
 954               fn));
 955 
 956       if(iabsload==2) NNEW(bold,double,neq[1]);
 957 
 958   /*  calculating the instantaneous loading vector at time 0 */
 959 
 960   NNEW(ikactmech,ITG,neq[1]);
 961   nactmech=0;
 962   FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
 963        ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
 964        nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,nbody,
 965        cgr,b,nactdof,&neq[1],nmethod,
 966        ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,alcon,
 967        nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,t0,t1act,
 968        ithermal,iprestr,vold,iperturb,iexpl,plicon,
 969        nplicon,plkcon,nplkcon,npmat_,ttime,&time0,istep,&iinc,&dtime,
 970        physcon,ibody,xbodyold,&reltime,veold,matname,mi,ikactmech,
 971        &nactmech,ielprop,prop));
 972   
 973   /*  correction for nonzero SPC's */
 974   
 975   if(iprescribedboundary){
 976 
 977       if(cyclicsymmetry){
 978           printf("*ERROR in dyna: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
 979           FORTRAN(stop,());
 980       }
 981 
 982       if(*idrct!=1){
 983           printf("*ERROR in dyna: variable increment length is not allwed in combination with prescribed boundaries\n");
 984           FORTRAN(stop,());
 985       }
 986       
 987       /* LU decomposition of the stiffness matrix */
 988       
 989       if(*isolver==0){
 990 #ifdef SPOOLES
 991           spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
 992                          &symmetryflag,&inputformat,&nzs[2]);
 993 #else
 994           printf("*ERROR in dyna: the SPOOLES library is not linked\n\n");
 995           FORTRAN(stop,());
 996 #endif
 997       }
 998       else if(*isolver==4){
 999 #ifdef SGI
1000           token=1;
1001           sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
1002 #else
1003           printf("*ERROR in dyna: the SGI library is not linked\n\n");
1004           FORTRAN(stop,());
1005 #endif
1006       }
1007       else if(*isolver==5){
1008 #ifdef TAUCS
1009           tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
1010 #else
1011           printf("*ERROR in dyna: the TAUCS library is not linked\n\n");
1012           FORTRAN(stop,());
1013 #endif
1014       }
1015       else if(*isolver==7){
1016 #ifdef PARDISO
1017           pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
1018                          &symmetryflag,&inputformat,jq,&nzs[2]);
1019 #else
1020           printf("*ERROR in dyna: the PARDISO library is not linked\n\n");
1021           FORTRAN(stop,());
1022 #endif
1023       }
1024 
1025       NNEW(bact,double,neq[1]);
1026       NNEW(bmin,double,neq[1]);
1027       NNEW(bv,double,neq[1]);
1028       NNEW(bprev,double,neq[1]);
1029       NNEW(bdiff,double,neq[1]);
1030 
1031       init=1;
1032       dynboun(amta,namta,nam,ampli,&time0,ttime,&dtime,xbounold,xboun,
1033               xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1034               aub,icol,irow,neq,nzs,&sigma,b,isolver,
1035               &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1036               bprev,bdiff,&nactmech,&iabsload,&iprev);
1037       init=0;
1038   }
1039 
1040 /* creating contact elements and calculating the contact forces
1041    (normal and shear) */
1042 
1043   if(ncont!=0){
1044       DMEMSET(bcont,0,neq[1],0.);
1045       contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1046               ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
1047               co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
1048               &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
1049               &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,iperturb,
1050               ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
1051               itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
1052               imastnode,nmastnode,xmastnor,filab,mcs,ics,
1053               &nasym,xnoels,&mortar,pslavsurf,pmastsurf,clearini,&theta);
1054 
1055       RENEW(ikactcont,ITG,nactcont_);
1056       DMEMSET(ikactcont,0,nactcont_,0.);
1057       nactcont=0;
1058      
1059       for(i=ne0;i<*ne;i++){
1060           indexe=ipkon[i];
1061           imat=ielmat[mi[2]*i];
1062           kodem=nelcon[2*imat-2];
1063           for(j=0;j<8;j++){lakonl[j]=lakon[8*i+j];}
1064           nope=atoi(&lakonl[7])+1;
1065           for(j=0;j<nope;j++){
1066               konl[j]=kon[indexe+j];
1067               for(j1=0;j1<3;j1++){
1068                   xl[j*3+j1]=co[3*(konl[j]-1)+j1];
1069                   voldl[mt*j+j1+1]=vold[mt*(konl[j]-1)+j1+1];
1070                   veoldl[mt*j+j1+1]=veold[mt*(konl[j]-1)+j1+1];
1071               }
1072           }
1073           konl[nope]=kon[indexe+nope];
1074 
1075           FORTRAN(springforc_n2f,(xl,konl,voldl,&imat,elcon,nelcon,elas,
1076               fnl,ncmat_,ntmat_,&nope,lakonl,&t1l,&kodem,elconloc,
1077               plicon,nplicon,npmat_,&senergy,&iener,cstr,mi,
1078               &springarea[2*(konl[nope]-1)],nmethod,&ne0,nstate_,
1079               xstateini,xstate,&reltime,&ielas));
1080 
1081           storecontactdof(&nope,nactdof,&mt,konl,&ikactcont,&nactcont,
1082                           &nactcont_,bcont,fnl,ikmpc,nmpc,ilmpc,ipompc,nodempc, 
1083                           coefmpc);
1084 
1085       }
1086       if(nactcont>100){nactcont_=nactcont;}else{nactcont_=100;}
1087       RENEW(ikactcont,ITG,nactcont_);
1088 
1089   }
1090 
1091   iit=1;
1092 
1093   /* load at the start of a new step:
1094      mechanical loading without contact */
1095 
1096   if(!cyclicsymmetry){
1097       for(i=0;i<nev;i++){
1098           i2=(long long)i*neq[1];
1099           aamech[i]=0.;
1100           if(nactmech<neq[1]/2){
1101               for(j=0;j<nactmech;j++){
1102                   aamech[i]+=z[i2+ikactmech[j]]*b[ikactmech[j]];
1103               }
1104           }else{
1105               for(j=0;j<neq[1];j++){
1106                   aamech[i]+=z[i2+j]*b[j];
1107               }
1108           }
1109           aanew[i]=aamech[i];
1110           if(ncont!=0){
1111               for(j=0;j<nactcont;j++){
1112                   aanew[i]+=z[i2+ikactcont[j]]*bcont[ikactcont[j]];
1113               }
1114           }
1115       }
1116   }else{
1117       for(i=0;i<nev;i++){aamech[i]=0.;}
1118       for(j=0;j<nactmech;j++){
1119           FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
1120           if(id!=0){
1121               if(izdof[id-1]==ikactmech[j]){
1122                   for(i=0;i<nev;i++){
1123                       aamech[i]+=z[(long long)i*nzdof+id-1]*b[ikactmech[j]];
1124                   }
1125               }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1126           }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1127       }
1128       memcpy(&aanew[0],&aamech[0],sizeof(double)*nev);
1129       if(ncont!=0){
1130           for(j=0;j<nactcont;j++){
1131               FORTRAN(nident,(izdof,&ikactcont[j],&nzdof,&id));
1132               if(id!=0){
1133                   if(izdof[id-1]==ikactcont[j]){
1134                       for(i=0;i<nev;i++){
1135                           aanew[i]+=z[(long long)i*nzdof+id-1]*bcont[ikactcont[j]];
1136                       }
1137                   }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1138               }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1139           }
1140       }
1141   }
1142 
1143   /* check whether integration point values are requested; if not,
1144      the stress fields do not have to be allocated */
1145 
1146   intpointvar=0;
1147   if(*ithermal<=1){
1148 
1149       /* mechanical */
1150 
1151       if((strcmp1(&filab[174],"S")==0)||
1152          (strcmp1(&filab[261],"E")==0)||
1153          (strcmp1(&filab[348],"RF")==0)||
1154          (strcmp1(&filab[435],"PEEQ")==0)||
1155          (strcmp1(&filab[522],"ENER")==0)||
1156          (strcmp1(&filab[609],"SDV")==0)||
1157          (strcmp1(&filab[1044],"ZZS")==0)||
1158          (strcmp1(&filab[1044],"ERR")==0)||
1159          (strcmp1(&filab[1479],"PHS")==0)||
1160          (strcmp1(&filab[1653],"MAXS")==0)||
1161          (strcmp1(&filab[2175],"CONT")==0)||
1162          (strcmp1(&filab[2262],"CELS")==0)) intpointvar=1;
1163       for(i=0;i<*nprint;i++){
1164           if((strcmp1(&prlab[6*i],"S")==0)||
1165              (strcmp1(&prlab[6*i],"E")==0)||
1166              (strcmp1(&prlab[6*i],"PEEQ")==0)||
1167              (strcmp1(&prlab[6*i],"ENER")==0)||
1168              (strcmp1(&prlab[6*i],"SDV")==0)||
1169              (strcmp1(&prlab[6*i],"CDIS")==0)||
1170              (strcmp1(&prlab[6*i],"CSTR")==0)||
1171              (strcmp1(&prlab[6*i],"CELS")==0)||
1172              (strcmp1(&prlab[6*i],"RF")==0)) {intpointvar=1;break;}
1173       }
1174   }else{
1175 
1176       /* thermal */
1177 
1178       if((strcmp1(&filab[696],"HFL")==0)||
1179          (strcmp1(&filab[783],"RFL")==0)) intpointvar=1;
1180       for(i=0;i<*nprint;i++){
1181           if((strcmp1(&prlab[6*i],"HFL")==0)||
1182              (strcmp1(&prlab[6*i],"RFL")==0)) {intpointvar=1;break;}
1183       }
1184   }
1185 
1186   if((intpointvar==1)) NNEW(stx,double,6*mi[0]**ne);
1187 
1188   /* major loop */
1189 
1190   resultmaxprev=0.;
1191   resultmax=0.;
1192 
1193   while(1.-theta>1.e-6){
1194     
1195     time0=time;
1196     
1197 //    printf("\nnew increment\n");
1198     
1199     if(*nener==1){
1200       memcpy(&enerini[0],&ener[0],sizeof(double)*mi[0]*ne0);
1201       if(*ithermal!=2){
1202           memcpy(&stiini[0],&sti[0],sizeof(double)*6*mi[0]*ne0);
1203       }
1204     }
1205 
1206     if(ncont!=0){
1207         if(nmdnode!=0){
1208             for(i=0;i<nmdnode;i++){
1209                 i1=mt*(imdnode[i]-1);
1210                 for(j=kmin;j<=kmax;j++){
1211                     vini[i1+j]=vold[i1+j];
1212                 }
1213             }
1214         }else{
1215             memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
1216         }
1217         if(*nstate_>0){
1218             for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
1219                 xstateini[k]=xstate[k];
1220             }
1221         }       
1222     }
1223     iinc++;
1224     jprint++;
1225 
1226     if(dashpot)RENEW(rpar,double,4+nev*(3+nev));
1227       
1228     /* check for max. # of increments */
1229     
1230     if(iinc>*jmax){
1231       printf(" *ERROR in dyna: max. # of increments reached\n\n");
1232       FORTRAN(stop,());
1233     }
1234     
1235       if(iinc>1){
1236           memcpy(&cd[0],&bj[0],sizeof(double)*nev);
1237           memcpy(&cv[0],&bjp[0],sizeof(double)*nev);
1238       }
1239 
1240       
1241       if((*idrct!=1)&&(iinc!=1)){
1242         
1243         /* increasing the increment size */
1244         
1245         dthetaold=dtheta;
1246         dtheta=dthetaref*dd;
1247         
1248         /* check increment length whether
1249            - it does not exceed tmax
1250            - the step length is not exceeded
1251            - a time point is not exceeded  */
1252         
1253         dthetaref=dtheta;
1254         checkinclength(&time0,ttime,&theta,&dtheta,idrct,tper,tmax,
1255                        tmin,ctrl, amta,namta,itpamp,&inext,&dthetaref,&itp,
1256                        &jprint,jout);
1257     }
1258       
1259       reltime=theta+dtheta;
1260       time=reltime**tper;
1261       dtime=dtheta**tper;
1262       
1263       /* calculating the instantaneous loads (forces, surface loading, 
1264          centrifugal and gravity loading or temperature) */
1265       
1266       FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
1267               xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1268               nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
1269               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
1270               nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
1271               ndirboun,nodeforc,
1272               ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
1273               nelemload,sideload,mi,
1274               xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
1275               &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
1276               fn));
1277               
1278       /* calculating the instantaneous loading vector */
1279               
1280       if(iabsload!=2){
1281           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1282                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcdiff,
1283                     nforc,nelemload,sideload,xloaddiff,nload,xbodydiff,
1284                     ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1285                     ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1286                     alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1287                     t0,t1diff,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1288                     nplicon,plkcon,nplkcon,
1289                     npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1290                     xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1291                     ielprop,prop));
1292       }else{
1293           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1294                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1295                     nforc,nelemload,sideload,xloadact,nload,xbodyact,
1296                     ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1297                     ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1298                     alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1299                     t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1300                     nplicon,plkcon,nplkcon,
1301                     npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1302                     xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1303                     ielprop,prop));
1304       }
1305               
1306       /* correction for nonzero SPC's */
1307       
1308       if(iprescribedboundary){
1309           dynboun(amta,namta,nam,ampli,&time,ttime,&dtime,
1310                       xbounold,xboun,
1311                       xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1312                       aub,icol,irow,neq,nzs,&sigma,b,isolver,
1313                       &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1314                       bprev,bdiff,&nactmech,&iabsload,&iprev);
1315       }
1316 
1317       if(*idrct==0){
1318           bbmax=0.;
1319           if(iabsload!=2){
1320               if(nactmech<neq[1]/2){
1321                   for(i=0;i<nactmech;i++){
1322                       if(fabs(b[ikactmech[i]])>bbmax) bbmax=fabs(b[ikactmech[i]]);
1323                   }
1324               }else{
1325                   for(i=0;i<neq[1];i++){
1326                       if(fabs(b[i])>bbmax) bbmax=fabs(b[i]);
1327                   }
1328               }
1329           }else{
1330 
1331               /* bbmax is to be calculated from the difference of b and bold */
1332 
1333               if(nactmech<neq[1]/2){
1334                   for(i=0;i<nactmech;i++){
1335                       if(fabs(b[ikactmech[i]]-bold[ikactmech[i]])>bbmax) 
1336                           bbmax=fabs(b[ikactmech[i]]-bold[ikactmech[i]]);
1337                   }
1338               }else{
1339                   for(i=0;i<neq[1];i++){
1340                       if(fabs(b[i])>bbmax) bbmax=fabs(b[i]-bold[i]);
1341                   }
1342               }
1343               
1344               /* copy b into bold */
1345 
1346               if(nactmech<neq[1]/2){
1347                   for(i=0;i<nactmech;i++){
1348                       bold[ikactmech[i]]=b[ikactmech[i]];
1349                   }
1350               }else{
1351                   memcpy(&bold[0],&b[0],sizeof(double)*neq[1]);
1352               }
1353           }
1354 
1355           /* check for size of mechanical force */
1356 
1357           if((bbmax>deltmx)&&(((itp==1)&&(dtheta>*tmin))||(itp==0))){
1358 
1359               /* force increase too big: increment size is decreased */
1360 
1361               if(iabsload==0) iabsload=1;
1362               dtheta=dtheta*deltmx/bbmax;
1363               dthetaref=dtheta;
1364               if(itp==1){
1365                   inext--;
1366                   itp=0;
1367               }
1368               
1369               /* check whether the new increment size is not too small */
1370               
1371               if(dtheta<*tmin){
1372                   dtheta=*tmin;
1373               }
1374 
1375               reltime=theta+dtheta;
1376               time=reltime**tper;
1377               dtime=dtheta**tper;
1378               
1379               /* calculating the instantaneous loads (forces, surface loading, 
1380                  centrifugal and gravity loading or temperature) */
1381               
1382               FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
1383                 xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1384                 nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
1385                 namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
1386                 nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
1387                 ndirboun,nodeforc,
1388                 ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
1389                 nelemload,sideload,mi,
1390                 xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
1391                 &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
1392                 fn));
1393               
1394               /* calculating the instantaneous loading vector */
1395               
1396               if(iabsload!=2){
1397                   FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1398                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcdiff,
1399                     nforc,nelemload,sideload,xloaddiff,nload,xbodydiff,
1400                     ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1401                     ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1402                     alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1403                     t0,t1diff,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1404                     nplicon,plkcon,nplkcon,
1405                     npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1406                     xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1407                     ielprop,prop));
1408               }else{
1409                   FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1410                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1411                     nforc,nelemload,sideload,xloadact,nload,xbodyact,
1412                     ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1413                     ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1414                     alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1415                     t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1416                     nplicon,plkcon,nplkcon,
1417                     npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1418                     xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1419                     ielprop,prop));
1420               }
1421               
1422               /* correction for nonzero SPC's */
1423               
1424               if(iprescribedboundary){
1425                   dynboun(amta,namta,nam,ampli,&time,ttime,&dtime,
1426                       xbounold,xboun,
1427                       xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1428                       aub,icol,irow,neq,nzs,&sigma,b,isolver,
1429                       &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1430                       bprev,bdiff,&nactmech,&iabsload,&iprev);
1431               }
1432               if(iabsload==1) iabsload=0;
1433           }
1434           
1435           if(ncont!=0){
1436               for(i=0;i<nactcont;i++){
1437                   jdof=ikactcont[i];
1438                   bcontini[jdof]=bcont[jdof];
1439               }
1440           }
1441 
1442       }
1443 
1444       /* step length is OK for mechanical load
1445          calculating equation for linearized loading */
1446 
1447       /* load: actual mechanical load +
1448                contact from last increment */
1449 
1450       if(!cyclicsymmetry){
1451           for(i=0;i<nev;i++){
1452               i2=(long long)i*neq[1];
1453               aa[i]=aanew[i];
1454               if(iabsload==2){aamech[i]=0.;}
1455               if(nactmech<neq[1]/2){
1456                   for(j=0;j<nactmech;j++){
1457                       aamech[i]+=z[i2+ikactmech[j]]*b[ikactmech[j]];
1458                   }
1459               }else{
1460                   for(j=0;j<neq[1];j++){
1461                       aamech[i]+=z[i2+j]*b[j];
1462                   }
1463               }
1464               
1465               aanew[i]=aamech[i];
1466               if(ncont!=0){
1467                   for(j=0;j<nactcont;j++){
1468                       aanew[i]+=z[i2+ikactcont[j]]*bcont[ikactcont[j]];
1469                   }
1470               }
1471               
1472               bb[i]=(aanew[i]-aa[i])/dtime;
1473               aa[i]=aanew[i]-bb[i]*time;
1474           }
1475       }else{
1476           for(i=0;i<nev;i++){
1477               memcpy(&aa[0],&aanew[0],sizeof(double)*nev);
1478               if(iabsload==2){aamech[i]=0.;}
1479           }
1480           for(j=0;j<nactmech;j++){
1481               FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
1482               if(id!=0){
1483                   if(izdof[id-1]==ikactmech[j]){
1484                       for(i=0;i<nev;i++){
1485                           aamech[i]+=z[(long long)i*nzdof+id-1]*b[ikactmech[j]];
1486                       }
1487                   }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1488               }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1489           }
1490           memcpy(&aanew[0],&aamech[0],sizeof(double)*nev);
1491           if(ncont!=0){
1492               for(j=0;j<nactcont;j++){
1493                   FORTRAN(nident,(izdof,&ikactcont[j],&nzdof,&id));
1494                   if(id!=0){
1495                       if(izdof[id-1]==ikactcont[j]){
1496                           for(i=0;i<nev;i++){
1497                               aanew[i]+=z[(long long)i*nzdof+id-1]*bcont[ikactcont[j]];
1498                           }
1499                       }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1500                   }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1501               }
1502           }
1503           for(i=0;i<nev;i++){
1504                 bb[i]=(aanew[i]-aa[i])/(dtime);
1505                 aa[i]=aanew[i]-bb[i]*time;
1506           }     
1507       }
1508       
1509       /* calculating the response due to unchanged contact force during
1510          the increment */
1511 
1512       if(dashpot){
1513           FORTRAN(subspace,(d,aa,bb,cc,&alpham,&betam,&nev,xini,
1514                             cd,cv,&time,rwork,&lrw,&iinc,jout,rpar,bj,
1515                             iwork,&liw,&iddebdf,bjp));
1516           if(iddebdf==2){
1517               liw=56+2*nev;
1518               RENEW(iwork,ITG,liw);
1519               for(i=0;i<liw;i++){iwork[i]=0;}
1520               lrw=250+20*nev+4*nev*nev;
1521               RENEW(rwork,double,lrw);
1522               for(i=0;i<lrw;i++){rwork[i]=0.;}
1523               iddebdf=1;
1524               FORTRAN(subspace,(d,aa,bb,cc,&alpham,&betam,&nev,xini,
1525                                 cd,cv,&time,rwork,&lrw,&iinc,jout,rpar,bj,
1526                                 iwork,&liw,&iddebdf,bjp));
1527           }
1528       }
1529       else{
1530           for(l=0;l<nev;l++){
1531               zetaj=zeta[l];
1532               dj=d[l];
1533               
1534               /* zero eigenfrequency: rigid body mode */
1535               
1536               if(fabs(d[l])<=1.e-10){
1537                   aai=aa[l];
1538                   bbi=bb[l];
1539                   tstart=time0;
1540                   tend=time;
1541                   sum=tend*(aai*time+
1542                             tend*((bbi*time-aai)/2.-bbi*tend/3.))-
1543                       tstart*(aai*time+
1544                               tstart*((bbi*time-aai)/2.-bbi*tstart/3.));
1545                   sump=tend*(aai+bbi*tend/2.)-tstart*(aai+bbi*tstart/2.);
1546                   bj[l]=sum+cd[l]+dtime*cv[l];
1547                   bjp[l]=sump+cv[l];
1548               }
1549               
1550               /*   subcritical damping */
1551               
1552               else if(zetaj<1.-1.e-6){
1553                   ddj=dj*sqrt(1.-zetaj*zetaj);
1554                   h1=zetaj*dj;
1555                   h2=h1*h1+ddj*ddj;
1556                   h3=h1*h1-ddj*ddj;
1557                   h4=2.*h1*ddj/h2;
1558                   h14=h1/ddj;
1559                   tstart=0.;                  
1560                   FORTRAN(fsub,(&time,&dtime,&aa[l],&bb[l],&ddj,
1561                                 &h1,&h2,&h3,&h4,&func,&funcp));
1562                   sum=func;sump=funcp;
1563                   FORTRAN(fsub,(&time,&tstart,&aa[l],&bb[l],&ddj,
1564                                 &h1,&h2,&h3,&h4,&func,&funcp));
1565                   sum-=func;sump-=funcp;
1566                   fexp=exp(-h1*dtime);
1567                   fsin=sin(ddj*dtime);
1568                   fcos=cos(ddj*dtime);
1569 
1570                   bj[l]=sum/ddj+fexp*(fcos+zetaj/sqrt(1.-zetaj*zetaj)*fsin)*cd[l]+
1571                       fexp*fsin*cv[l]/ddj;
1572                   bjp[l]=sump/ddj+fexp*((-h1+ddj*h14)*fcos+(-ddj-h1*h14)*fsin)*cd[l]
1573                     +fexp*(-h1*fsin+ddj*fcos)*cv[l]/ddj;
1574 
1575               }
1576               
1577               /*      supercritical damping */
1578               
1579               else if(zetaj>1.+1.e-6){
1580                   ddj=dj*sqrt(zetaj*zetaj-1.);
1581                   h1=ddj-zetaj*dj;
1582                   h2=ddj+zetaj*dj;
1583                   h3=1./h1;
1584                   h4=1./h2;
1585                   h5=h3*h3;
1586                   h6=h4*h4;
1587                   tstart=0.;                  
1588                   FORTRAN(fsuper,(&time,&dtime,&aa[l],&bb[l],
1589                                   &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
1590                   sum=func;sump=funcp;
1591                   FORTRAN(fsuper,(&time,&tstart,&aa[l],&bb[l],
1592                                   &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
1593                   sum-=func;sump-=funcp;
1594                   
1595                   fexm=exp(h1*dtime);
1596                   fexp=exp(-h2*dtime);
1597                   h14=zetaj*dj/ddj;
1598                   bj[l]=sum/(2.*ddj)+(fexm+fexp)*cd[l]/2.+zetaj*(fexm-fexp)/(2.*
1599                       sqrt(zetaj*zetaj-1.))*cd[l]+(fexm-fexp)*cv[l]/(2.*ddj);
1600                   bjp[l]=sump/(2.*ddj)+(h1*fexm-h2*fexp)*cd[l]/2.
1601                     +(h14*cd[l]+cv[l]/ddj)*(h1*fexm+h2*fexp)/2.;
1602               }
1603               
1604               /* critical damping */
1605               
1606               else{
1607                   h1=zetaj*dj;
1608                   h2=1./h1;
1609                   h3=h2*h2;
1610                   h4=h2*h3;
1611                   tstart=0.;
1612                   FORTRAN(fcrit,(&time,&dtime,&aa[l],&bb[l],&zetaj,&dj,
1613                                  &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
1614                   sum=func;sump=funcp;
1615                   FORTRAN(fcrit,(&time,&tstart,&aa[l],&bb[l],&zetaj,&dj,
1616                                  &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
1617                   sum-=func;sump-=funcp;
1618                   fexp=exp(-h1*dtime);
1619                   bj[l]=sum+fexp*((1.+h1*dtime)*cd[l]+dtime*cv[l]);
1620                   bjp[l]=sump+fexp*(-h1*h1*dtime*cd[l]+
1621                                     (1.-h1*dtime)*cv[l]);
1622               }
1623           }
1624       }
1625       
1626       /* composing the response */
1627       
1628       if(iprescribedboundary){
1629           if(nmdnode==0){
1630               memcpy(&b[0],&bmin[0],sizeof(double)*neq[1]);
1631               memcpy(&bp[0],&bv[0],sizeof(double)*neq[1]);
1632           }else{
1633               for(i=0;i<nmddof;i++){
1634                   b[imddof[i]]=bmin[imddof[i]];
1635                   bp[imddof[i]]=bv[imddof[i]];
1636               }
1637           }
1638       }
1639       else{
1640           if(nmdnode==0){
1641               DMEMSET(b,0,neq[1],0.);
1642               DMEMSET(bp,0,neq[1],0.);
1643           }else{
1644               for(i=0;i<nmddof;i++){
1645                   b[imddof[i]]=0.;
1646                   bp[imddof[i]]=0.;
1647               }
1648           }
1649       }
1650       
1651       if(!cyclicsymmetry){
1652           if(nmdnode==0){
1653               for(i=0;i<neq[1];i++){
1654                   for(j=0;j<nev;j++){
1655                       b[i]+=bj[j]*z[(long long)j*neq[1]+i];
1656                       bp[i]+=bjp[j]*z[(long long)j*neq[1]+i];
1657                   }
1658               }
1659           }else{
1660               for(i=0;i<nmddof;i++){
1661                   for(j=0;j<nev;j++){
1662                       b[imddof[i]]+=bj[j]*z[(long long)j*neq[1]+imddof[i]];
1663                       bp[imddof[i]]+=bjp[j]*z[(long long)j*neq[1]+imddof[i]];
1664                   }
1665               }
1666           }
1667       }else{
1668           for(i=0;i<nmddof;i++){
1669               FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1670               if(id!=0){
1671                   if(izdof[id-1]==imddof[i]){
1672                       for(j=0;j<nev;j++){
1673                           b[imddof[i]]+=bj[j]*z[(long long)j*nzdof+id-1];
1674                           bp[imddof[i]]+=bjp[j]*z[(long long)j*nzdof+id-1];
1675                       }
1676                   }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1677               }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1678           }
1679       }
1680       
1681       /* update nonlinear MPC-coefficients (e.g. for rigid
1682          body MPC's */
1683 
1684       if(inonlinmpc==1){
1685           FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1686                          nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1687                          &maxlenmpc,ikmpc,ilmpc,&icascade,
1688                          kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1689                          &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1690       }
1691       
1692       /* calculating displacements/temperatures */
1693       
1694       FORTRAN(dynresults,(nk,v,ithermal,nactdof,vold,nodeboun,
1695               ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1696               b,bp,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1697               imdmpc,&nmdmpc,nmethod,&time));
1698       
1699       /* creating contact elements and calculating the contact forces
1700          based on the displacements at the end of the present increment */
1701       
1702       if(ncont!=0){
1703           dynacont(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
1704               ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
1705               nforc,nelemload,sideload,xload,nload,nactdof,neq,nzl,icol,
1706               irow,
1707               nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
1708               nrhcon,cocon,ncocon,alcon,nalcon,alzero,ielmat,ielorien,
1709               norien,orab,ntmat_,t0,t1,ithermal,prestr,iprestr,
1710               vold,iperturb,sti,nzs,tinc,tper,xmodal,veold,amname,amta,
1711               namta,nam,iamforc,iamload,iamt1,jout,filab,eme,xforcold,
1712               xloadold,t1old,iamboun,xbounold,iexpl,plicon,nplicon,plkcon,
1713               nplkcon,xstate,npmat_,matname,mi,ncmat_,nstate_,ener,jobnamec,
1714               ttime,set,nset,istartset,iendset,ialset,nprint,prlab,
1715               prset,nener,trab,inotr,ntrans,fmpc,cbody,ibody,xbody,nbody,
1716               xbodyold,istep,isolver,jq,output,mcs,nkon,mpcend,ics,cs,ntie,
1717               tieset,idrct,jmax,tmin,tmax,ctrl,itpamp,tietol,&iit,
1718               &ncont,&ne0,&reltime,&dtime,bcontini,bj,aux,iaux,bcont,
1719               &nev,v,&nkon0,&deltmx,&dtheta,&theta,&iprescribedboundary,
1720               &mpcfree,&memmpc_,itietri,koncont,cg,straight,&iinc,
1721               vini,aa,bb,aanew,d,z,zeta,b,&time0,&time,ipobody,
1722               xforcact,xloadact,t1act,xbounact,xbodyact,cd,cv,ampli,
1723               &dthetaref,bjp,bp,cstr,imddof,&nmddof, 
1724               &ikactcont,&nactcont,&nactcont_,aamech,bprev,&iprev,&inonlinmpc,
1725               &ikactmech,&nactmech,imdnode,&nmdnode,imdboun,&nmdboun,
1726               imdmpc,&nmdmpc,&itp,&inext,imastop,
1727               nslavnode,islavnode,islavsurf,itiefac,areaslav,iponoels,
1728               inoels,springarea,izdof,&nzdof,fn,imastnode,nmastnode,xmastnor,
1729               xstateini,nslavs,&cyclicsymmetry,xnoels,&ielas,ielprop,prop);
1730       }   
1731           
1732       theta+=dtheta;
1733 //      (*ttime)+=dtime;
1734       
1735       /* check whether a time point was reached */
1736       
1737       if((*itpamp>0)&&(*idrct==0)){
1738           if(itp==1){
1739               jprint=*jout;
1740           }else{
1741               jprint=*jout+1;
1742           }
1743       }
1744       
1745       /* check whether output is needed */
1746       
1747       if((*jout==jprint)||(1.-theta<=1.e-6)){
1748           iout=2;
1749           jprint=0;
1750       }else if((*nener==1)){
1751           iout=-2;
1752       }else{
1753           iout=0;
1754       }
1755       
1756       if((iout==2)||(iout==-2)){
1757 
1758         /* deactivating the elements for which the stresses are not
1759            needed */
1760 
1761         if(nmdnode>0){
1762             if((intpointvar==1)){
1763                 for(k=0;k<ne0;k++){
1764                     if(ipkon[k]<-1){
1765                         printf("*ERROR in dyna: contact remeshing of quadratic elements is not allowed\n\n");
1766                         FORTRAN(stop,());
1767                     }else if(ipkon[k]!=-1){
1768                         ipkon[k]=-ipkon[k]-2;
1769                     }
1770                 }
1771                 for(k=0;k<nmdelem;k++){
1772                     ielem=imdelem[k]-1;
1773                     ipkon[ielem]=-2-ipkon[ielem];
1774                 }
1775             }
1776         }
1777 
1778         results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1779                 stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1780                 ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1781                 ithermal,prestr,iprestr,filab,eme,emn,een,
1782                 iperturb,f,fn,nactdof,&iout,qa,
1783                 vold,b,nodeboun,ndirboun,xbounact,nboun,
1784                 ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1785                 veold,accold,&bet,&gam,&dtime,&time,ttime,
1786                 plicon,nplicon,plkcon,nplkcon,
1787                 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1788                 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1789                 enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1790                 set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1791                 qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,ikmpc,
1792                 ilmpc,istep,&iinc,springarea,&reltime,&ne0,xforc,nforc,
1793                 thicke,shcon,nshcon,
1794                 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1795                 &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1796                 islavsurf,ielprop,prop);
1797 
1798         /* restoring */
1799 
1800         if(nmdnode>0){
1801             if((intpointvar==1)){
1802                 for(k=0;k<ne0;k++){
1803                     if(ipkon[k]<-1){ipkon[k]=-2-ipkon[k];}
1804                 }
1805             }
1806         }
1807         
1808         if((*ithermal!=2)&&(intpointvar==1)){
1809           for(k=0;k<6*mi[0]*ne0;++k){
1810             sti[k]=stx[k];
1811           }
1812         }
1813       }
1814       if(iout==2){
1815         (*kode)++;
1816         if(strcmp1(&filab[1044],"ZZS")==0){
1817             NNEW(neigh,ITG,40**ne);
1818             NNEW(ipneigh,ITG,*nk);
1819         }
1820 
1821         ptime=*ttime+time;
1822         frd(co,&nkg,kon,ipkon,lakon,&neg,v,stn,inum,nmethod,
1823             kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1824             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1825             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1826             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1827             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1828             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1829         
1830         if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1831       }
1832       
1833       if(isteadystate==1){
1834           
1835           /* calculate maximum displacement/temperature */
1836           
1837           resultmax=0.;
1838           if(*ithermal<2){
1839               for(i=1;i<mt**nk;i=i+mt){
1840                   if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1841               for(i=2;i<mt**nk;i=i+mt){
1842                   if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1843               for(i=3;i<mt**nk;i=i+mt){
1844                   if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1845           }else if(*ithermal==2){
1846               for(i=0;i<mt**nk;i=i+mt){
1847                   if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1848           }else{
1849               printf("*ERROR in dyna: coupled temperature-displacement calculations are not allowed\n");
1850           }
1851           if(fabs((resultmax-resultmaxprev)/resultmax)<precision){
1852               break;
1853           }else{resultmaxprev=resultmax;}
1854       }
1855      
1856   }
1857       
1858   if((intpointvar==1)) SFREE(stx);
1859 
1860   /* calculating the displacements and velocities in all nodes as 
1861      initial condition for the next step; only needed if
1862      - nonzero initial conditions are allowed (-> no cyclic symmetry)
1863      - the output was restricted (-> nmdnode nonzero) */
1864 
1865   if((nmdnode!=0)&&(!cyclicsymmetry)){
1866 
1867       /* determining the solution in the independent nodes */
1868 
1869       DMEMSET(b,0,neq[1],0.);
1870       DMEMSET(bp,0,neq[1],0.);
1871 
1872       for(i=0;i<neq[1];i++){
1873           for(j=0;j<nev;j++){
1874               b[i]+=bj[j]*z[(long long)j*neq[1]+i];
1875               bp[i]+=bjp[j]*z[(long long)j*neq[1]+i];
1876           }
1877       }
1878       
1879       /* update nonlinear MPC-coefficients (e.g. for rigid
1880          body MPC's */
1881 
1882       if(inonlinmpc==1){
1883           FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1884                          nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1885                          &maxlenmpc,ikmpc,ilmpc,&icascade,
1886                          kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1887                          &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1888       }
1889       
1890       /* calculating displacements/temperatures */
1891       
1892       nmdnode=0;
1893       FORTRAN(dynresults,(nk,v,ithermal,nactdof,vold,nodeboun,
1894               ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1895               b,bp,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1896               imdmpc,&nmdmpc,nmethod,&time));
1897   }
1898   
1899   SFREE(eei);
1900   SFREE(vbounact);
1901   SFREE(abounact);
1902 
1903   if(*nener==1){SFREE(stiini);SFREE(enerini);}
1904 
1905   if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
1906   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1907   if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
1908   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
1909   if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
1910 
1911   /* updating the loading at the end of the step; 
1912      important in case the amplitude at the end of the step
1913      is not equal to one */
1914 
1915   for(k=0;k<*nboun;++k){xboun[k]=xbounact[k];}
1916   for(k=0;k<*nforc;++k){xforc[k]=xforcact[k];}
1917   for(k=0;k<2**nload;++k){xload[k]=xloadact[k];}
1918   for(k=0;k<7**nbody;k=k+7){xbody[k]=xbodyact[k];}
1919   if(*ithermal==1){
1920     for(k=0;k<*nk;++k){t1[k]=t1act[k];}
1921   }
1922   
1923   SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(adb);SFREE(d);
1924   SFREE(aub);SFREE(z);SFREE(b);SFREE(zeta);SFREE(bj);SFREE(cd);SFREE(cv);
1925   SFREE(xforcact);SFREE(xloadact);SFREE(xbounact);SFREE(aa);SFREE(bb);SFREE(aanew);
1926   SFREE(ampli);SFREE(xbodyact);SFREE(bjp);SFREE(bp);SFREE(aamech);SFREE(ikactmech);
1927   SFREE(xforcdiff);SFREE(xloaddiff);SFREE(xboundiff),SFREE(xbodydiff);
1928 
1929   if(*ithermal==1) {SFREE(t1act);SFREE(t1diff);}
1930 
1931   if(iprescribedboundary){
1932       if(*isolver==0){
1933 #ifdef SPOOLES
1934           spooles_cleanup();
1935 #endif
1936       }
1937       else if(*isolver==4){
1938 #ifdef SGI
1939           sgi_cleanup(token);
1940 #endif
1941       }
1942       else if(*isolver==5){
1943 #ifdef TAUCS
1944           tau_cleanup();
1945 #endif
1946       }
1947       else if(*isolver==7){
1948 #ifdef PARDISO
1949           pardiso_cleanup(&neq[1],&symmetryflag);
1950 #endif
1951       }
1952       SFREE(bact);SFREE(bmin);SFREE(bv);SFREE(bprev);SFREE(bdiff);
1953   }
1954 
1955   /* deleting the contact information */
1956 
1957 //  *ne=ne0; *nkon=nkon0;
1958   if(ncont!=0){
1959       *ne=ne0; *nkon=nkon0;
1960       if(*nener==1){
1961         RENEW(ener,double,mi[0]**ne*2);
1962       }
1963       RENEW(ipkon,ITG,*ne);
1964       RENEW(lakon,char,8**ne);
1965       RENEW(kon,ITG,*nkon);
1966       if(*norien>0){
1967           RENEW(ielorien,ITG,mi[2]**ne);
1968       }
1969       RENEW(ielmat,ITG,mi[2]**ne);
1970       SFREE(cg);SFREE(straight);
1971 
1972       SFREE(vini);SFREE(bcont);SFREE(bcontini);SFREE(ikactcont);
1973 
1974       SFREE(imastop);SFREE(itiefac);SFREE(islavsurf);SFREE(islavnode);
1975       SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
1976       SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
1977       SFREE(springarea);SFREE(xmastnor);
1978 
1979       SFREE(areaslav);
1980 
1981       if(*nstate_>0){SFREE(xstateini);}
1982 
1983   }
1984 
1985   if(!cyclicsymmetry){
1986       SFREE(ad);SFREE(au);
1987   }else{
1988       SFREE(adbe); SFREE(aube);SFREE(icole); SFREE(irowe); SFREE(jqe);SFREE(izdof);
1989       SFREE(nm);
1990 
1991       *nk/=nsectors;
1992       *ne/=nsectors;
1993       *nkon/=nsectors;
1994       *nboun/=nsectors;
1995       neq[1]=neq[1]*2/nsectors;
1996 
1997       RENEW(ialset,ITG,nalset_);
1998 
1999       /* restore the infomration in istartset and iendset */
2000 
2001       for(j=0; j<*nset; j++){
2002         istartset[j]=istartset_[j];
2003         iendset[j]=iendset_[j];
2004       }
2005       SFREE(istartset_);
2006       SFREE(iendset_);
2007 
2008       RENEW(co,double,3**nk);
2009       if((*ithermal!=0)&&(*nam>0)) RENEW(iamt1,ITG,*nk);
2010       RENEW(nactdof,ITG,mt**nk);
2011       if(*ntrans>0) RENEW(inotr,ITG,2**nk);
2012       RENEW(kon,ITG,*nkon);
2013       RENEW(ipkon,ITG,*ne);
2014       RENEW(lakon,char,8**ne);
2015       RENEW(ielmat,ITG,mi[2]**ne);
2016       if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne);
2017       RENEW(nodeboun,ITG,*nboun);
2018       RENEW(ndirboun,ITG,*nboun);
2019       if(*nam>0) RENEW(iamboun,ITG,*nboun);
2020       RENEW(xboun,double,*nboun);
2021       RENEW(xbounold,double,*nboun);
2022       RENEW(ikboun,ITG,*nboun);
2023       RENEW(ilboun,ITG,*nboun);
2024 
2025       /* recovering the original multiple point constraints */
2026 
2027       RENEW(ipompc,ITG,*nmpc);
2028       RENEW(nodempc,ITG,3**mpcend);
2029       RENEW(coefmpc,double,*mpcend);
2030       RENEW(labmpc,char,20**nmpc+1);
2031       RENEW(ikmpc,ITG,*nmpc);
2032       RENEW(ilmpc,ITG,*nmpc);
2033       RENEW(fmpc,double,*nmpc);
2034 
2035       *nmpc=nmpcold;
2036       *mpcend=mpcendold;
2037       for(i=0;i<*nmpc;i++){ipompc[i]=ipompcold[i];}
2038       for(i=0;i<3**mpcend;i++){nodempc[i]=nodempcold[i];}
2039       for(i=0;i<*mpcend;i++){coefmpc[i]=coefmpcold[i];}
2040       for(i=0;i<20**nmpc;i++){labmpc[i]=labmpcold[i];}
2041       for(i=0;i<*nmpc;i++){ikmpc[i]=ikmpcold[i];}
2042       for(i=0;i<*nmpc;i++){ilmpc[i]=ilmpcold[i];}
2043       SFREE(ipompcold);SFREE(nodempcold);SFREE(coefmpcold);
2044       SFREE(labmpcold);SFREE(ikmpcold);SFREE(ilmpcold);
2045 
2046       RENEW(vold,double,mt**nk);
2047       RENEW(veold,double,mt**nk);
2048       RENEW(eme,double,6*mi[0]**ne);
2049       RENEW(sti,double,6*mi[0]**ne);
2050       if(*nener==1)RENEW(ener,double,mi[0]**ne*2);
2051 
2052 /* distributed loads */
2053 
2054       for(i=0;i<*nload;i++){
2055           if(nelemload[2*i]<nsectors){
2056               nelemload[2*i]-=*ne*nelemload[2*i+1];
2057           }else{
2058               nelemload[2*i]-=*ne*(nelemload[2*i+1]-nsectors);
2059           }
2060       }
2061 
2062   /*  sorting the elements with distributed loads */
2063 
2064       if(*nload>0){
2065           if(*nam>0){
2066               FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
2067           }else{
2068               FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
2069           }
2070       }
2071       
2072 /* point loads */
2073       
2074       for(i=0;i<*nforc;i++){
2075           if(nodeforc[2*i+1]<nsectors){
2076               nodeforc[2*i]-=*nk*nodeforc[2*i+1];
2077           }else{
2078               nodeforc[2*i]-=*nk*(nodeforc[2*i+1]-nsectors);
2079           }
2080       }
2081   }
2082 
2083   SFREE(xstiff);if(*nbody>0) SFREE(ipobody);
2084 
2085   if(dashpot){
2086       SFREE(xini);SFREE(rwork);SFREE(adc);SFREE(auc);SFREE(cc);
2087       SFREE(rpar);SFREE(iwork);}
2088 
2089   SFREE(cstr);
2090 
2091   SFREE(imddof);SFREE(imdnode);SFREE(imdboun);SFREE(imdmpc);SFREE(imdelem);
2092 
2093   if(iabsload==2) SFREE(bold);
2094 
2095   *ialsetp=ialset;
2096   *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
2097   *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
2098   *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;
2099   *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
2100   *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
2101   *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2102   *fmpcp=fmpc;*veoldp=veold;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
2103   *stip=sti;*xstatep=xstate;
2104 
2105   (*ttime)+=(*tper);
2106 
2107   return;
2108 }
2109 

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