root/src/nonlingeo.c

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

DEFINITIONS

This source file includes following definitions.
  1. nonlingeo

   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 "CalculiX.h"
  22 #ifdef SPOOLES
  23    #include "spooles.h"
  24 #endif
  25 #ifdef SGI
  26    #include "sgi.h"
  27 #endif
  28 #ifdef TAUCS
  29    #include "tau.h"
  30 #endif
  31 #ifdef PARDISO
  32    #include "pardiso.h"
  33 #endif
  34 
  35 
  36 void nonlingeo(double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp,
  37              ITG *ne, 
  38              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  39              ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp,
  40              ITG *nmpc, 
  41              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  42              ITG *nelemload, char *sideload, double *xload,ITG *nload, 
  43              ITG *nactdof, 
  44              ITG **icolp, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, 
  45              ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG *ikboun, 
  46              ITG *ilboun,
  47              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  48              double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
  49              ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
  50              double *t0, double *t1, double *t1old, 
  51              ITG *ithermal,double *prestr, ITG *iprestr, 
  52              double **voldp,ITG *iperturb, double *sti, ITG *nzs,  
  53              ITG *kode, char *filab, 
  54              ITG *idrct, ITG *jmax, ITG *jout, double *timepar,
  55              double *eme,
  56              double *xbounold, double *xforcold, double *xloadold,
  57              double *veold, double *accold,
  58              char *amname, double *amta, ITG *namta, ITG *nam,
  59              ITG *iamforc, ITG *iamload,
  60              ITG *iamt1, double *alpha, ITG *iexpl,
  61              ITG *iamboun, double *plicon, ITG *nplicon, double *plkcon,
  62              ITG *nplkcon,
  63              double **xstatep, ITG *npmat_, ITG *istep, double *ttime,
  64              char *matname, double *qaold, ITG *mi,
  65              ITG *isolver, ITG *ncmat_, ITG *nstate_, ITG *iumat,
  66              double *cs, ITG *mcs, ITG *nkon, double **enerp, ITG *mpcinfo,
  67              char *output,
  68              double *shcon, ITG *nshcon, double *cocon, ITG *ncocon,
  69              double *physcon, ITG *nflow, double *ctrl, 
  70              char *set, ITG *nset, ITG *istartset,
  71              ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  72              char *prset, ITG *nener,ITG *ikforc,ITG *ilforc, double *trab, 
  73              ITG *inotr, ITG *ntrans, double **fmpcp, char *cbody,
  74              ITG *ibody, double *xbody, ITG *nbody, double *xbodyold,
  75              ITG *ielprop, double *prop, ITG *ntie, char *tieset,
  76              ITG *itpamp, ITG *iviewfile, char *jobnamec, double *tietol,
  77              ITG *nslavs, double *thicke, ITG *ics, 
  78              ITG *nintpoint,ITG *mortar,ITG *ifacecount,char *typeboun,
  79              ITG **islavsurfp,double **pslavsurfp,double **clearinip,
  80              ITG *nmat,double *xmodal,ITG *iaxial,ITG *inext){
  81 
  82   char description[13]="            ",*lakon=NULL,jobnamef[396]="",
  83       *sideface=NULL,*labmpc=NULL,fnffrd[132]="",*lakonf=NULL; 
  84  
  85   ITG *inum=NULL,k,iout=0,icntrl,iinc=0,jprint=0,iit=-1,jnz=0,
  86       icutb=0,istab=0,ifreebody,uncoupled,n1,n2,itruecontact,
  87       iperturb_sav[2],ilin,*icol=NULL,*irow=NULL,ielas=0,icmd=0,
  88       memmpc_,mpcfree,icascade,maxlenmpc,*nodempc=NULL,*iaux=NULL,
  89       *nodempcref=NULL,memmpcref_,mpcfreeref,*itg=NULL,*ineighe=NULL,
  90       *ieg=NULL,ntg=0,ntr,*kontri=NULL,*nloadtr=NULL,idamping=0,
  91       *ipiv=NULL,ntri,newstep,mode=-1,noddiam=-1,nasym=0,im,
  92       ntrit,*inocs=NULL,inewton=0,*ipobody=NULL,*nacteq=NULL,
  93       *nactdog=NULL,nteq,network,*itietri=NULL,*koncont=NULL,
  94       ncont,ne0,nkon0,*ipkon=NULL,*kon=NULL,*ielorien=NULL,
  95       *ielmat=NULL,itp=0,symmetryflag=0,inputformat=0,
  96       *iruc=NULL,iitterm=0,iturbulent,ngraph=1,ismallsliding=0,
  97       *ipompc=NULL,*ikmpc=NULL,*ilmpc=NULL,i0ref,irref,icref,
  98       *itiefac=NULL,*islavsurf=NULL,*islavnode=NULL,*imastnode=NULL,
  99       *nslavnode=NULL,*nmastnode=NULL,*imastop=NULL,imat,
 100       *iponoels=NULL,*inoels=NULL,*islavsurfold=NULL,maxlenmpcref,
 101       *islavact=NULL,mt=mi[1]+1,*nactdofinv=NULL,*ipe=NULL, 
 102       *ime=NULL,*ikactmech=NULL,nactmech,inode,idir,neold,
 103       iemchange=0,nzsrad,*mast1rad=NULL,*irowrad=NULL,*icolrad=NULL,
 104       *jqrad=NULL,*ipointerrad=NULL,*integerglob=NULL,negpres=0,
 105       mass[2]={0,0}, stiffness=1, buckling=0, rhsi=1, intscheme=0,idiscon=0,
 106       coriolis=0,*ipneigh=NULL,*neigh=NULL,maxprevcontel,nslavs_prev_step,
 107       *nelemface=NULL,*ipoface=NULL,*nodface=NULL,*ifreestream=NULL,iex,
 108       *isolidsurf=NULL,*neighsolidsurf=NULL,*iponoel=NULL,*inoel=NULL,
 109       nef=0,nface,nfreestream,nsolidsurf,i,indexe,icfd=0,id,*neij=NULL,
 110       node,networknode,iflagact=0,*nodorig=NULL,*ipivr=NULL,iglob=0,
 111       *inomat=NULL,*ipnei=NULL,ntrimax,*nx=NULL,*ny=NULL,*nz=NULL,
 112       *neifa=NULL,*neiel=NULL,*ielfa=NULL,*ifaext=NULL,nflnei,nfaext,
 113       idampingwithoutcontact=0,*nactdoh=NULL,*nactdohinv=NULL,*ipkonf=NULL,
 114       *ielmatf=NULL,*ielorienf=NULL;
 115 
 116   double *stn=NULL,*v=NULL,*een=NULL,cam[5],*epn=NULL,*cg=NULL,
 117          *cdn=NULL,*vel=NULL,*vfa=NULL,*pslavsurfold=NULL,
 118          *f=NULL,*fn=NULL,qa[3]={0.,0.,-1.},qam[2]={0.,0.},dtheta,theta,
 119          err,ram[8]={0.,0.,0.,0.,0.,0.,0.,0.},*areaslav=NULL,
 120          *springarea=NULL,ram1[8]={0.,0.,0.,0.,0.,0.,0.,0.},
 121          ram2[8]={0.,0.,0.,0.,0.,0.,0.,0.},deltmx,ptime,
 122          uam[2]={0.,0.},*vini=NULL,*ac=NULL,qa0,qau,ea,*straight=NULL,
 123          *t1act=NULL,qamold[2],*xbounact=NULL,*bc=NULL,
 124          *xforcact=NULL,*xloadact=NULL,*fext=NULL,*clearini=NULL,
 125          reltime,time,bet=0.,gam=0.,*aux2=NULL,dtime,*fini=NULL,
 126          *fextini=NULL,*veini=NULL,*accini=NULL,*xstateini=NULL,
 127          *ampli=NULL,scal1,*eei=NULL,*t1ini=NULL,pressureratio,
 128          *xbounini=NULL,dev,*xstiff=NULL,*stx=NULL,*stiini=NULL,
 129          *enern=NULL,*coefmpc=NULL,*aux=NULL,*xstaten=NULL,
 130          *coefmpcref=NULL,*enerini=NULL,*emn=NULL,alpham,betam,
 131          *tarea=NULL,*tenv=NULL,*erad=NULL,*fnr=NULL,*fni=NULL,
 132          *adview=NULL,*auview=NULL,*qfx=NULL,*cvini=NULL,*cv=NULL,
 133          *qfn=NULL,*co=NULL,*vold=NULL,*fenv=NULL,sigma=0.,
 134          *xbodyact=NULL,*cgr=NULL,dthetaref, *vcontu=NULL,*vr=NULL,*vi=NULL,
 135          *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,*fmpc=NULL,*ener=NULL,
 136          *f_cm=NULL, *f_cs=NULL,*adc=NULL,*auc=NULL,
 137          *xstate=NULL,*eenmax=NULL,*adrad=NULL,*aurad=NULL,*bcr=NULL,
 138          *xmastnor=NULL,*emeini=NULL,*tinc,*tper,*tmin,*tmax,*tincf,
 139          *doubleglob=NULL,*xnoels=NULL,*au=NULL,
 140          *ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
 141          *x=NULL,*y=NULL,*z=NULL,*xo=NULL,
 142          *yo=NULL,*zo=NULL,*cdnr=NULL,*cdni=NULL;
 143 
 144 #ifdef SGI
 145   ITG token;
 146 #endif
 147   
 148   icol=*icolp;irow=*irowp;co=*cop;vold=*voldp;
 149   ipkon=*ipkonp;lakon=*lakonp;kon=*konp;ielorien=*ielorienp;
 150   ielmat=*ielmatp;ener=*enerp;xstate=*xstatep;
 151   
 152   ipompc=*ipompcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 153   fmpc=*fmpcp;nodempc=*nodempcp;coefmpc=*coefmpcp;
 154 
 155   islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
 156 
 157   tinc=&timepar[0];
 158   tper=&timepar[1];
 159   tmin=&timepar[2];
 160   tmax=&timepar[3];
 161   tincf=&timepar[4];
 162 
 163   if(*ithermal==4){
 164       uncoupled=1;
 165       *ithermal=3;
 166   }else{
 167       uncoupled=0;
 168   }
 169   
 170   if(*mortar!=1){
 171       maxprevcontel=*nslavs;
 172   }else if(*mortar==1){
 173       maxprevcontel=*nintpoint;
 174       if(*nstate_!=0){
 175           if(maxprevcontel!=0){
 176               NNEW(islavsurfold,ITG,2**ifacecount+2);
 177               NNEW(pslavsurfold,double,3**nintpoint);
 178               memcpy(&islavsurfold[0],&islavsurf[0],
 179                      sizeof(ITG)*(2**ifacecount+2));
 180               memcpy(&pslavsurfold[0],&pslavsurf[0],
 181                      sizeof(double)*(3**nintpoint));
 182           }
 183       }
 184   }
 185   nslavs_prev_step=*nslavs;
 186 
 187   /* turbulence model 
 188      iturbulent==0: laminar
 189      iturbulent==1: k-epsilon
 190      iturbulent==2: q-omega
 191      iturbulent==3: SST */
 192   
 193   iturbulent=(ITG)physcon[8];
 194   
 195   for(k=0;k<3;k++){
 196       strcpy1(&jobnamef[k*132],&jobnamec[k*132],132);
 197   }
 198   
 199   qa0=ctrl[20];qau=ctrl[21];ea=ctrl[23];deltmx=ctrl[26];
 200   i0ref=ctrl[0];irref=ctrl[1];icref=ctrl[3];
 201   
 202   memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
 203   maxlenmpc=mpcinfo[3];
 204 
 205   alpham=xmodal[0];
 206   betam=xmodal[1];
 207 
 208   /* check whether, for a dynamic calculation, damping is involved */
 209 
 210   if(*nmethod==4){
 211       if(*iexpl<=1){
 212           if((fabs(alpham)>1.e-30)||(fabs(betam)>1.e-30)){
 213               idamping=1;
 214               idampingwithoutcontact=1;
 215           }else{
 216               for(i=0;i<*ne;i++){
 217                   if(ipkon[i]<0) continue;
 218                   if(strcmp1(&lakon[i*8],"ED")==0){
 219                       idamping=1;idampingwithoutcontact=1;break;
 220                   }
 221               }
 222           }
 223       }
 224   }
 225 
 226   if((icascade==2)&&(*iexpl>=2)){
 227       printf(" *ERROR in nonlingeo: linear and nonlinear MPC's depend on each other\n");
 228       printf("        This is not allowed in a explicit dynamic calculation\n");
 229       FORTRAN(stop,());
 230   }
 231   
 232   /* check whether the submodel is meant for a fluid-structure
 233      interaction */
 234   
 235   strcpy(fnffrd,jobnamec);
 236   strcat(fnffrd,"f.frd");
 237   if((jobnamec[396]!='\0')&&(strcmp1(&jobnamec[396],fnffrd)==0)){
 238       
 239       /* fluid-structure interaction: wait till after the compfluid
 240          call */
 241       
 242       NNEW(integerglob,ITG,1);
 243       NNEW(doubleglob,double,1);
 244   }else{
 245       
 246       /* determining the global values to be used as boundary conditions
 247          for a submodel */
 248       
 249       getglobalresults(jobnamec,&integerglob,&doubleglob,nboun,iamboun,xboun,
 250                        nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
 251                        ithermal,nk,t1,iamt1);
 252   }
 253   
 254   /* invert nactdof */
 255   
 256   NNEW(nactdofinv,ITG,mt**nk);
 257   NNEW(nodorig,ITG,*nk);
 258   FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
 259                          ipkon,lakon,kon,ne));
 260   SFREE(nodorig);
 261   
 262   /* allocating a field for the stiffness matrix */
 263   
 264   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 265   
 266   /* allocating force fields */
 267   
 268   NNEW(f,double,neq[1]);
 269   NNEW(fext,double,neq[1]);
 270   
 271   NNEW(b,double,neq[1]);
 272   NNEW(vini,double,mt**nk);
 273   
 274   NNEW(aux,double,7*maxlenmpc);
 275   NNEW(iaux,ITG,2*maxlenmpc);
 276   
 277   /* allocating fields for the actual external loading */
 278   
 279   NNEW(xbounact,double,*nboun);
 280   NNEW(xbounini,double,*nboun);
 281   for(k=0;k<*nboun;++k){xbounact[k]=xbounold[k];}
 282   NNEW(xforcact,double,*nforc);
 283   NNEW(xloadact,double,2**nload);
 284   NNEW(xbodyact,double,7**nbody);
 285   /* copying the rotation axis and/or acceleration vector */
 286   for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
 287   
 288   /* assigning the body forces to the elements */ 
 289   
 290   if(*nbody>0){
 291       ifreebody=*ne+1;
 292       NNEW(ipobody,ITG,2*ifreebody**nbody);
 293       for(k=1;k<=*nbody;k++){
 294           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 295                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 296           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 297       }
 298       RENEW(ipobody,ITG,2*(ifreebody-1));
 299       if(inewton==1){NNEW(cgr,double,4**ne);}
 300   }
 301   
 302   /* for mechanical calculations: updating boundary conditions
 303      calculated in a previous thermal step */
 304   
 305   if(*ithermal<2) FORTRAN(gasmechbc,(vold,nload,sideload,
 306                                      nelemload,xload,mi));
 307   
 308   /* for thermal calculations: forced convection and cavity
 309      radiation*/
 310   
 311   if(*ithermal>1){
 312       NNEW(itg,ITG,*nload+3**nflow);
 313       NNEW(ieg,ITG,*nflow);
 314       /* max 6 triangles per face, 4 entries per triangle */
 315       NNEW(kontri,ITG,24**nload);
 316       NNEW(nloadtr,ITG,*nload);
 317       NNEW(nacteq,ITG,4**nk);
 318       NNEW(nactdog,ITG,4**nk);
 319       NNEW(v,double,mt**nk);
 320       FORTRAN(envtemp,(itg,ieg,&ntg,&ntr,sideload,nelemload,
 321                        ipkon,kon,lakon,ielmat,ne,nload,
 322                        kontri,&ntri,nloadtr,nflow,ndirboun,nactdog,
 323                        nodeboun,nacteq,nboun,ielprop,prop,&nteq,
 324                        v,&network,physcon,shcon,ntmat_,co,
 325                        vold,set,nshcon,rhcon,nrhcon,mi,nmpc,nodempc,
 326                        ipompc,labmpc,ikboun,&nasym,iaxial));
 327       SFREE(v);
 328       
 329       if((*mcs>0)&&(ntr>0)){
 330           NNEW(inocs,ITG,*nk);
 331           radcyc(nk,kon,ipkon,lakon,ne,cs,mcs,nkon,ialset,istartset,
 332                  iendset,&kontri,&ntri,&co,&vold,&ntrit,inocs,mi);
 333       }
 334       else{ntrit=ntri;}
 335       
 336       nzsrad=100*ntr;
 337       NNEW(mast1rad,ITG,nzsrad);
 338       NNEW(irowrad,ITG,nzsrad);
 339       NNEW(icolrad,ITG,ntr);
 340       NNEW(jqrad,ITG,ntr+1);
 341       NNEW(ipointerrad,ITG,ntr);
 342       
 343       if(ntr>0){
 344           mastructrad(&ntr,nloadtr,sideload,ipointerrad,
 345                       &mast1rad,&irowrad,&nzsrad,
 346                       jqrad,icolrad);
 347       }
 348       
 349       SFREE(ipointerrad);SFREE(mast1rad);
 350       RENEW(irowrad,ITG,nzsrad);
 351       
 352       RENEW(itg,ITG,ntg);
 353       NNEW(ineighe,ITG,ntg);
 354       RENEW(kontri,ITG,4*ntrit);
 355       RENEW(nloadtr,ITG,ntr);
 356       
 357       NNEW(adview,double,ntr);
 358       NNEW(auview,double,2*nzsrad);
 359       NNEW(tarea,double,ntr);
 360       NNEW(tenv,double,ntr);
 361       NNEW(fenv,double,ntr);
 362       NNEW(erad,double,ntr);
 363       
 364       NNEW(ac,double,nteq*nteq);
 365       NNEW(bc,double,nteq);
 366       NNEW(ipiv,ITG,nteq);
 367       NNEW(adrad,double,ntr);
 368       NNEW(aurad,double,2*nzsrad);
 369       NNEW(bcr,double,ntr);
 370       NNEW(ipivr,ITG,ntr);
 371   }
 372   
 373   /* check for fluid elements */
 374   
 375   NNEW(nactdoh,ITG,*ne);
 376   NNEW(nactdohinv,ITG,*ne);
 377   for(i=0;i<*ne;++i){
 378       if(ipkon[i]<0) continue;
 379       indexe=ipkon[i];
 380       if(strcmp1(&lakon[8*i],"F")==0){
 381           icfd=1;nactdohinv[nef]=i+1;nef++;nactdoh[i]=nef;}
 382   }
 383   if(icfd==1){
 384       RENEW(nactdohinv,ITG,nef);
 385       NNEW(ipkonf,ITG,nef);
 386       NNEW(lakonf,char,8*nef);
 387       NNEW(ielmatf,ITG,mi[2]*nef);
 388       if(*norien>0) NNEW(ielorienf,ITG,mi[2]*nef);
 389       NNEW(ipoface,ITG,*nk);
 390       NNEW(nodface,ITG,5*6*nef);
 391       NNEW(ipnei,ITG,*ne);
 392       NNEW(neifa,ITG,6**ne);
 393       NNEW(neiel,ITG,6**ne);
 394       NNEW(neij,ITG,6**ne);
 395       NNEW(ielfa,ITG,24**ne);
 396       NNEW(ifaext,ITG,6**ne);
 397       NNEW(isolidsurf,ITG,6**ne);
 398       NNEW(vel,double,6**ne);
 399       NNEW(vfa,double,6*6**ne);
 400 
 401       FORTRAN(precfd,(ne,ipkon,kon,lakon,ipnei,neifa,neiel,ipoface,
 402        nodface,ielfa,&nflnei,&nface,ifaext,&nfaext,
 403        isolidsurf,&nsolidsurf,set,nset,istartset,iendset,ialset,
 404        vel,vold,mi,neij,&nef,nactdoh,ipkonf,lakonf,ielmatf,ielmat,
 405        ielorienf,ielorien,norien));
 406 
 407       SFREE(ipoface);
 408       SFREE(nodface);
 409       RENEW(neifa,ITG,nflnei);
 410       RENEW(neiel,ITG,nflnei);
 411       RENEW(neij,ITG,nflnei);
 412       RENEW(ielfa,ITG,4*nface);
 413       RENEW(ifaext,ITG,nfaext);
 414       RENEW(isolidsurf,ITG,nsolidsurf);
 415       RENEW(vfa,double,6*nface);
 416   }else{
 417       SFREE(nactdoh);SFREE(nactdohinv);
 418   }
 419   if(*ithermal>1){NNEW(qfx,double,3*mi[0]**ne);}
 420   
 421   /* contact conditions */
 422   
 423   inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
 424           lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
 425           &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
 426           mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
 427           nmpc,&mpcfree,&memmpc_,
 428           &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 429           iperturb,ikboun,nboun,co,istep,&xnoels);
 430   
 431   if(ncont!=0){
 432       
 433       NNEW(cg,double,3*ncont);
 434       NNEW(straight,double,16*ncont);
 435           
 436       /* 11 instead of 10: last position is reserved for the
 437          local contact spring element number; needed as
 438          pointer into springarea */
 439       
 440       if(*mortar==0){
 441           RENEW(kon,ITG,*nkon+11**nslavs);
 442           NNEW(springarea,double,2**nslavs);
 443           if(*nener==1){
 444               RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
 445           }
 446           RENEW(ipkon,ITG,*ne+*nslavs);
 447           RENEW(lakon,char,8*(*ne+*nslavs));
 448           
 449           if(*norien>0){
 450               RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
 451               for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
 452           }
 453 
 454           RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
 455           for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
 456 
 457           if((maxprevcontel==0)&&(*nslavs!=0)){
 458               RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
 459               for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
 460                   xstate[k]=0.;
 461               }
 462           }
 463           maxprevcontel=*nslavs;
 464 
 465           NNEW(areaslav,double,*ifacecount);
 466       }else if(*mortar==1){
 467           NNEW(islavact,ITG,nslavnode[*ntie]);
 468           if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
 469 
 470           /* check whether at least one contact definition involves true contact
 471              and not just tied contact */
 472 
 473           FORTRAN(checktruecontact,(ntie,tieset,tietol,elcon,&itruecontact,
 474                                     ncmat_,ntmat_));
 475       }
 476 
 477       NNEW(xmastnor,double,3*nmastnode[*ntie]);
 478   }
 479   
 480   if((icascade==2)||(ncont!=0)){
 481       memmpcref_=memmpc_;mpcfreeref=mpcfree;maxlenmpcref=maxlenmpc;
 482       NNEW(nodempcref,ITG,3*memmpc_);
 483       for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
 484       NNEW(coefmpcref,double,memmpc_);
 485       for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
 486   }
 487   
 488   if((*ithermal==1)||(*ithermal>=3)){
 489       NNEW(t1ini,double,*nk);
 490       NNEW(t1act,double,*nk);
 491       for(k=0;k<*nk;++k){t1act[k]=t1old[k];}
 492   }
 493   
 494   /* allocating a field for the instantaneous amplitude */
 495   
 496   NNEW(ampli,double,*nam);
 497   
 498   /* fini is also needed in static calculations if ilin=1
 499      to get correct values of f after a divergent increment */
 500 
 501   NNEW(fini,double,neq[1]);
 502   
 503   /* allocating fields for nonlinear dynamics */
 504   
 505   if(*nmethod==4){
 506       mass[0]=1;
 507       mass[1]=1;
 508       NNEW(aux2,double,neq[1]);
 509 //      NNEW(fini,double,neq[1]);
 510       NNEW(fextini,double,neq[1]);
 511       NNEW(veini,double,mt**nk);
 512       NNEW(accini,double,mt**nk);
 513       NNEW(adb,double,neq[1]);
 514       NNEW(aub,double,nzs[1]);
 515       NNEW(cvini,double,neq[1]);
 516       NNEW(cv,double,neq[1]);
 517   }
 518 
 519   if((*nstate_!=0)&&((*mortar!=1)||(ncont==0))){
 520       NNEW(xstateini,double,*nstate_*mi[0]*(*ne+*nslavs));
 521       for(k=0;k<*nstate_*mi[0]*(*ne+*nslavs);++k){
 522           xstateini[k]=xstate[k];
 523       }
 524   }
 525   if((*nstate_!=0)&&(*mortar==1)) NNEW(xstateini,double,1);
 526   NNEW(eei,double,6*mi[0]**ne);
 527   NNEW(stiini,double,6*mi[0]**ne);
 528   NNEW(emeini,double,6*mi[0]**ne);
 529   if(*nener==1)
 530       NNEW(enerini,double,mi[0]**ne);
 531   
 532   qa[0]=qaold[0];
 533   qa[1]=qaold[1];
 534   
 535   /* normalizing the time */
 536   
 537   FORTRAN(checktime,(itpamp,namta,tinc,ttime,amta,tmin,inext,&itp,istep));
 538   dtheta=(*tinc)/(*tper);
 539 
 540   /* taking care of a small increment at the end of the step
 541      for face-to-face penalty contact */
 542 
 543   dthetaref=dtheta;
 544   if((dtheta<=1.e-6)&&(*iexpl<=1)){
 545       printf("\n *ERROR in nonlingeo\n");
 546       printf(" increment size smaller than one millionth of step size\n");
 547       printf(" increase increment size\n\n");
 548   }
 549   *tmin=*tmin/(*tper);
 550   *tmax=*tmax/(*tper);
 551   theta=0.;
 552   
 553   /* calculating an initial flux norm */
 554   
 555   if(*ithermal!=2){
 556       if(qau>1.e-10){qam[0]=qau;}
 557       else if(qa0>1.e-10){qam[0]=qa0;}
 558       else if(qa[0]>1.e-10){qam[0]=qa[0];}
 559       else {qam[0]=1.e-2;}
 560   }
 561   if(*ithermal>1){
 562       if(qau>1.e-10){qam[1]=qau;}
 563       else if(qa0>1.e-10){qam[1]=qa0;}
 564       else if(qa[1]>1.e-10){qam[1]=qa[1];}
 565       else {qam[1]=1.e-2;}
 566   }
 567   
 568   /* storing the element and topology information before introducing 
 569      contact elements */
 570   
 571   ne0=*ne;nkon0=*nkon;neold=*ne;
 572   
 573   /*********************************************************************/
 574   
 575   /* calculating of the acceleration due to force discontinuities
 576      (external - internal force) at the start of a step */
 577   
 578   /*********************************************************************/
 579   
 580   if((*nmethod==4)&&(*ithermal!=2)){
 581       bet=(1.-*alpha)*(1.-*alpha)/4.;
 582       gam=0.5-*alpha;
 583       
 584       /* calculating the stiffness and mass matrix */
 585       
 586       reltime=0.;
 587       time=0.;
 588       dtime=0.;
 589       
 590       FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
 591               xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
 592               t1old,t1,t1act,iamt1,nk,amta,
 593               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
 594               xbounold,xboun,xbounact,iamboun,nboun,
 595               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 596               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 597               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 598               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 599       
 600       time=0.;
 601       dtime=1.;
 602     
 603       /*  updating the nonlinear mpc's (also affects the boundary
 604           conditions through the nonhomogeneous part of the mpc's)
 605           if contact arises the number of MPC's can also change */
 606       
 607       cam[0]=0.;cam[1]=0.;cam[2]=0.;
 608       
 609       if(icascade==2){
 610           memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
 611           RENEW(nodempc,ITG,3*memmpcref_);
 612           for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
 613           RENEW(coefmpc,double,memmpcref_);
 614           for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
 615       }
 616 
 617       newstep=0;
 618       FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
 619                          nmpc,ikboun,ilboun,nboun,xbounold,aux,iaux,
 620                          &maxlenmpc,ikmpc,ilmpc,&icascade,
 621                          kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
 622                          &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
 623       if(icascade==2){
 624 //        memmpcref_=memmpc_;mpcfreeref=mpcfree;
 625 //        RENEW(nodempcref,ITG,3*memmpc_);
 626           for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
 627 //        RENEW(coefmpcref,double,memmpc_);
 628           for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
 629       }
 630       
 631       if(icascade>0) remastruct(ipompc,&coefmpc,&nodempc,nmpc,
 632               &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
 633               labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
 634               kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
 635               neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
 636               &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
 637               typeboun,&cv,&cvini,&iit);
 638       
 639       /* invert nactdof */
 640 
 641       SFREE(nactdofinv);
 642       NNEW(nactdofinv,ITG,1);
 643       
 644       iout=-1;
 645       ielas=1;
 646       
 647       NNEW(fn,double,mt**nk);
 648       NNEW(stx,double,6*mi[0]**ne);
 649       
 650       NNEW(inum,ITG,*nk);
 651       results(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,stx,
 652               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 653               ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
 654               prestr,iprestr,filab,eme,emn,een,iperturb,
 655               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 656               ndirboun,xbounold,nboun,ipompc,
 657               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,&bet,
 658               &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 659               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 660               ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,emeini,xstaten,
 661               eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 662               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 663               nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
 664               &ne0,xforc,nforc,thicke,shcon,nshcon,
 665               sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 666               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 667               islavsurf,ielprop,prop);
 668       
 669       SFREE(fn);SFREE(stx);SFREE(inum);
 670       
 671       iout=0;
 672       ielas=0;
 673       
 674       reltime=0.;
 675       time=0.;
 676       dtime=0.;
 677       
 678       if(*iexpl<=1){intscheme=1;}
 679       
 680       /* in mafillsm the stiffness and mass matrix are computed;
 681          The primary aim is to calculate the mass matrix (not 
 682          lumped for an implicit dynamic calculation, lumped for an
 683          explicit dynamic calculation). However:
 684          - for an implicit calculation the mass matrix is "doped" with
 685          a small amount of stiffness matrix, therefore the calculation
 686          of the stiffness matrix is needed.
 687          - for an explicit calculation the stiffness matrix is not 
 688          needed at all. Since the calculation of the mass matrix alone
 689          is not possible in mafillsm, the determination of the stiffness
 690          matrix is taken as unavoidable "ballast". */
 691       
 692       NNEW(ad,double,neq[1]);
 693       NNEW(au,double,nzs[1]);
 694       
 695       FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounold,nboun,
 696               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
 697               nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
 698               nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
 699               nmethod,ikmpc,ilmpc,ikboun,ilboun,
 700               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
 701               ielmat,ielorien,norien,orab,ntmat_,
 702               t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
 703               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 704               xstiff,npmat_,&dtime,matname,mi,
 705               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
 706               physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
 707               &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
 708               xstateini,xstate,thicke,integerglob,doubleglob,
 709               tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 710               pmastsurf,mortar,clearini,ielprop,prop));
 711       
 712       if(*nmethod==0){
 713           
 714           /* error occurred in mafill: storing the geometry in frd format */
 715           
 716           ++*kode;
 717           if(strcmp1(&filab[1044],"ZZS")==0){
 718               NNEW(neigh,ITG,40**ne);
 719               NNEW(ipneigh,ITG,*nk);
 720           }
 721           
 722           ptime=*ttime+time;
 723           frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
 724               kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
 725               nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
 726               ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 727               mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 728               cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 729               thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
 730           
 731           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}      
 732 #ifdef COMPANY
 733           FORTRAN(uout,(v,mi,ithermal));
 734 #endif    
 735           FORTRAN(stop,());
 736           
 737       }
 738       
 739       /* mass x acceleration = f(external)-f(internal) 
 740          only for the mechanical loading*/
 741       
 742       for(k=0;k<neq[0];++k){
 743           b[k]=fext[k]-f[k];
 744       }
 745       
 746       if(*iexpl<=1){
 747           
 748           /* a small amount of stiffness is added to the mass matrix
 749              otherwise the system leads to huge accelerations in 
 750              case of discontinuous load changes at the start of the step */
 751           
 752           dtime=*tinc/10.;
 753           scal1=bet*dtime*dtime*(1.+*alpha);
 754           for(k=0;k<neq[0];++k){
 755               ad[k]=adb[k]+scal1*ad[k];
 756           }
 757           for(k=0;k<nzs[0];++k){
 758               au[k]=aub[k]+scal1*au[k];
 759           }
 760           if(*isolver==0){
 761 #ifdef SPOOLES
 762               spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
 763                       &symmetryflag,&inputformat,&nzs[2]);
 764 #else
 765               printf(" *ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
 766               FORTRAN(stop,());
 767 #endif
 768           }
 769           else if((*isolver==2)||(*isolver==3)){
 770               preiter(ad,&au,b,&icol,&irow,&neq[0],&nzs[0],isolver,iperturb);
 771           }
 772           else if(*isolver==4){
 773 #ifdef SGI
 774               token=1;
 775               sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
 776 #else
 777               printf(" *ERROR in nonlingeo: the SGI library is not linked\n\n");
 778               FORTRAN(stop,());
 779 #endif
 780           }
 781           else if(*isolver==5){
 782 #ifdef TAUCS
 783               tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[0],&nzs[0]);
 784 #else
 785               printf(" *ERROR in nonlingeo: the TAUCS library is not linked\n\n");
 786               FORTRAN(stop,());
 787 #endif
 788           }
 789           else if(*isolver==7){
 790 #ifdef PARDISO
 791               pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
 792                            &symmetryflag,&inputformat,jq,&nzs[2]);
 793 #else
 794               printf(" *ERROR in nonlingeo: the PARDISO library is not linked\n\n");
 795               FORTRAN(stop,());
 796 #endif
 797           }
 798       }
 799       
 800       else{
 801           for(k=0;k<neq[0];++k){
 802               b[k]=(fext[k]-f[k])/adb[k];
 803           }
 804       }
 805       
 806       /* for thermal loading the acceleration is set to zero */
 807       
 808       for(k=neq[0];k<neq[1];++k){
 809           b[k]=0.;
 810       }
 811       
 812       /* calculating the displacements, stresses and forces */
 813       
 814       NNEW(v,double,mt**nk);
 815       memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
 816       
 817       NNEW(stx,double,6*mi[0]**ne);
 818       NNEW(fn,double,mt**nk);
 819       
 820       NNEW(inum,ITG,*nk);
 821       dtime=1.e-20;
 822       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 823               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 824               ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
 825               prestr,iprestr,filab,eme,emn,een,iperturb,
 826               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 827               ndirboun,xbounact,nboun,ipompc,
 828               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 829               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 830               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 831               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 832               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 833               iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 834               fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 835               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 836               sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 837               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 838               islavsurf,ielprop,prop);
 839       SFREE(inum);
 840       dtime=0.;
 841 
 842       memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
 843       if(*ithermal!=2){
 844           for(k=0;k<6*mi[0]*ne0;++k){
 845               sti[k]=stx[k];
 846           }
 847       }
 848 
 849       SFREE(v);SFREE(stx);SFREE(fn);
 850       SFREE(ad);SFREE(au);
 851       
 852       /* the mass matrix is kept for subsequent calculations, therefore,
 853          no new mass calculation is necessary for the remaining iterations
 854          in the present step */
 855       
 856       mass[0]=0;intscheme=0;
 857 /*      DMEMSET(f,0,neq[1],0.);
 858         DMEMSET(fext,0,neq[1],0.);*/
 859       
 860   }
 861   
 862   if(*iexpl>1) icmd=3;
 863   
 864   /**************************************************************/
 865   /* starting the loop over the increments                      */
 866   /**************************************************************/
 867   
 868   newstep=1;
 869   
 870   while((1.-theta>1.e-6)||(negpres==1)){
 871       
 872       if(icutb==0){
 873           
 874           /* previous increment converged: update the initial values */
 875           
 876           iinc++;
 877           jprint++;
 878           
 879           /* vold is copied into vini */
 880           
 881           memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
 882           
 883           for(k=0;k<*nboun;++k){xbounini[k]=xbounact[k];}
 884           if((*ithermal==1)||(*ithermal>=3)){
 885               for(k=0;k<*nk;++k){t1ini[k]=t1act[k];}
 886           }
 887           for(k=0;k<neq[1];++k){
 888               fini[k]=f[k];
 889           }
 890           if(*nmethod==4){
 891               for(k=0;k<mt**nk;++k){
 892                   veini[k]=veold[k];
 893                   accini[k]=accold[k];
 894               }
 895               for(k=0;k<neq[1];++k){
 896 //                fini[k]=f[k];
 897                   fextini[k]=fext[k];
 898                   cvini[k]=cv[k];
 899               }
 900           }
 901           if(*ithermal!=2){
 902               for(k=0;k<6*mi[0]*ne0;++k){
 903                   stiini[k]=sti[k];
 904                   emeini[k]=eme[k];
 905               }
 906           }
 907           if(*nener==1)
 908               for(k=0;k<mi[0]*ne0;++k){enerini[k]=ener[k];}
 909 
 910           if(*mortar!=1){
 911               if(*nstate_!=0){
 912                   for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
 913                       xstateini[k]=xstate[k];
 914                   }
 915               }
 916           }     
 917       }
 918       
 919       /* check for max. # of increments */
 920       
 921       if(iinc>*jmax){
 922           printf(" *ERROR: max. # of increments reached\n\n");
 923           FORTRAN(stop,());
 924       }
 925       printf(" increment %" ITGFORMAT " attempt %" ITGFORMAT " \n",iinc,icutb+1);
 926       printf(" increment size= %e\n",dtheta**tper);
 927       printf(" sum of previous increments=%e\n",theta**tper);
 928       printf(" actual step time=%e\n",(theta+dtheta)**tper);
 929       printf(" actual total time=%e\n\n",*ttime+(theta+dtheta)**tper);
 930       
 931       printf(" iteration 1\n\n");
 932       
 933       qamold[0]=qam[0];
 934       qamold[1]=qam[1];
 935       
 936       /* determining the actual loads at the end of the new increment*/
 937       
 938       reltime=theta+dtheta;
 939       time=reltime**tper;
 940       dtime=dtheta**tper;
 941       
 942       FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
 943               xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
 944               t1old,t1,t1act,iamt1,nk,amta,
 945               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
 946               xbounold,xboun,xbounact,iamboun,nboun,
 947               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 948               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 949               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 950               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 951       
 952       for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
 953       if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,bcr,ipivr,
 954        ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,ntmat_,vold,
 955        shcon,nshcon,ipkon,kon,co,
 956        kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
 957        nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
 958        cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
 959        ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
 960        nactdog,nacteq,nodeboun,ndirboun,&network,
 961        rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
 962        ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,nset,
 963        ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,iamload,
 964        jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
 965       }
 966       
 967       if(icfd==1){
 968           compfluid(&co,nk,&ipkonf,&kon,&lakonf,&sideface,
 969             ifreestream,&nfreestream,isolidsurf,neighsolidsurf,&nsolidsurf,
 970             &iponoel,&inoel,nshcon,shcon,nrhcon,rhcon,&vold,ntmat_,nodeboun,
 971             ndirboun,nboun,&ipompc,&nodempc,nmpc,&ikmpc,&ilmpc,ithermal,
 972             ikboun,ilboun,&iturbulent,isolver,iexpl,vcontu,ttime,
 973             &time,&dtime,nodeforc,ndirforc,xforc,nforc,nelemload,sideload,
 974             xload,nload,xbody,ipobody,nbody,ielmatf,matname,mi,ncmat_,
 975             physcon,istep,&iinc,ibody,xloadold,xboun,&coefmpc,
 976             nmethod,xforcold,xforcact,iamforc,iamload,xbodyold,xbodyact,
 977             t1old,t1,t1act,iamt1,amta,namta,nam,ampli,xbounold,xbounact,
 978             iamboun,itg,&ntg,amname,t0,&nelemface,&nface,cocon,ncocon,xloadact,
 979             tper,jmax,jout,set,nset,istartset,iendset,ialset,prset,prlab,
 980             nprint,trab,inotr,ntrans,filab,&labmpc,sti,norien,orab,jobnamef,
 981             tieset,ntie,mcs,ics,cs,nkon,&mpcfree,&memmpc_,&fmpc,&nef,&inomat,
 982             qfx,neifa,neiel,ielfa,ifaext,vfa,vel,ipnei,&nflnei,&nfaext,
 983             typeboun,neij,tincf,nactdoh,nactdohinv,ielorienf);
 984 
 985           /* determining the global values to be used as boundary conditions
 986              for a submodel */
 987           
 988           SFREE(integerglob);SFREE(doubleglob);
 989           getglobalresults(jobnamec,&integerglob,&doubleglob,nboun,iamboun,
 990                            xboun,nload,sideload,iamload,&iglob,nforc,iamforc,
 991                            xforc,ithermal,nk,t1,iamt1);
 992       }
 993       
 994       if((icascade==2)||(ncont!=0)){
 995           memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
 996           RENEW(nodempc,ITG,3*memmpcref_);
 997           for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
 998           RENEW(coefmpc,double,memmpcref_);
 999           for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
1000       }
1001 
1002       /* generating contact elements */
1003       
1004       if((ncont!=0)&&(*mortar<=1)){
1005           *ne=ne0;*nkon=nkon0;
1006 
1007           /* at start of new increment: 
1008              - copy state variables (node-to-face)
1009              - determine slave integration points (face-to-face)
1010              - interpolate state variables (face-to-face) */
1011 
1012           if(icutb==0){
1013               if(*mortar==1){
1014 
1015                   if(*nstate_!=0){
1016                       if(maxprevcontel!=0){
1017                           if(iit!=-1){
1018                               NNEW(islavsurfold,ITG,2**ifacecount+2);
1019                               NNEW(pslavsurfold,double,3**nintpoint);
1020                               memcpy(&islavsurfold[0],&islavsurf[0],
1021                                      sizeof(ITG)*(2**ifacecount+2));
1022                               memcpy(&pslavsurfold[0],&pslavsurf[0],
1023                                      sizeof(double)*(3**nintpoint));
1024                           }
1025                       }
1026                   }
1027 
1028                   *nintpoint=0;
1029 
1030                   precontact(&ncont,ntie,tieset,nset,set,istartset,
1031                      iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
1032                      cg,straight,co,vold,istep,&iinc,&iit,itiefac,
1033                      islavsurf,islavnode,imastnode,nslavnode,nmastnode,
1034                      imastop,mi,ipe,ime,tietol,&iflagact,
1035                      nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
1036                      nslavs);
1037                   
1038                   /* changing the dimension of element-related fields */
1039                   
1040                   RENEW(kon,ITG,*nkon+22**nintpoint);
1041                   RENEW(springarea,double,2**nintpoint);
1042                   RENEW(pmastsurf,double,6**nintpoint);
1043                   
1044                   if(*nener==1){
1045                       RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
1046                   }
1047                   RENEW(ipkon,ITG,*ne+*nintpoint);
1048                   RENEW(lakon,char,8*(*ne+*nintpoint));
1049                   
1050                   if(*norien>0){
1051                       RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
1052                       for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
1053                   }
1054                   RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
1055                   for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
1056 
1057               }
1058           }
1059 
1060           contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1061                   ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
1062                   co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
1063                   &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
1064                   &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
1065                   iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,
1066                   islavsurf,
1067                   itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
1068                   imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
1069                   xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta);
1070 
1071           /* check whether, for a dynamic calculation, damping is involved */
1072 
1073           if(*nmethod==4){
1074               if(*iexpl<=1){
1075                   if(idampingwithoutcontact==0){
1076                       for(i=0;i<*ne;i++){
1077                           if(ipkon[i]<0) continue;
1078                           if(*ncmat_>=5){
1079                               if(strcmp1(&lakon[i*8],"ES")==0){
1080                                   if(strcmp1(&lakon[i*8+6],"C")==0){
1081                                       imat=ielmat[i*mi[2]];
1082                                       if(elcon[(*ncmat_+1)**ntmat_*(imat-1)+4]>0.){
1083                                           idamping=1;break;
1084                                       }
1085                                   }
1086                               }
1087                           }
1088                       }
1089                   }
1090               }
1091           }
1092           
1093           printf(" Number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
1094                   
1095           /* interpolating the state variables */
1096 
1097           if(icutb==0){
1098               if(*mortar==1){
1099                   if(*nstate_!=0){
1100                       if(maxprevcontel!=0){
1101                           RENEW(xstateini,double,
1102                                 *nstate_*mi[0]*(ne0+maxprevcontel));
1103                           for(k=*nstate_*mi[0]*ne0;
1104                                 k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
1105                               xstateini[k]=xstate[k];
1106                           }
1107                       }
1108                       
1109                       RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
1110                       for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
1111                           xstate[k]=0.;
1112                       }
1113                       
1114                       if((*nintpoint>0)&&(maxprevcontel>0)){
1115                           iex=2;
1116                           
1117                           /* interpolation of xstate */
1118                           
1119                           FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
1120                                &ne0,mi,xstate,pslavsurf,nstate_,
1121                                xstateini,islavsurf,islavsurfold,
1122                                pslavsurfold));
1123                           
1124                       }
1125 
1126                       if(maxprevcontel!=0){
1127                           SFREE(islavsurfold);SFREE(pslavsurfold);
1128                       }
1129 
1130                       maxprevcontel=*nintpoint;
1131 
1132                       RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
1133                       for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
1134                           xstateini[k]=xstate[k];
1135                       }
1136                   }
1137               }
1138           }
1139           
1140           if(icascade<1)icascade=1;
1141       }
1142       
1143       /*  updating the nonlinear mpc's (also affects the boundary
1144           conditions through the nonhomogeneous part of the mpc's) */
1145       
1146       FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1147                          nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1148                          &maxlenmpc,ikmpc,ilmpc,&icascade,
1149                          kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1150                          &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1151       
1152       if((icascade==2)||(ncont!=0)){
1153 //        memmpcref_=memmpc_;mpcfreeref=mpcfree;
1154 //        RENEW(nodempcref,ITG,3*memmpc_);
1155           for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
1156 //        RENEW(coefmpcref,double,memmpc_);
1157           for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
1158       }
1159       
1160       if(icascade>0) remastruct(ipompc,&coefmpc,&nodempc,nmpc,
1161           &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
1162           labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
1163           kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
1164           neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
1165           &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
1166           typeboun,&cv,&cvini,&iit);
1167       
1168       /* invert nactdof */
1169       
1170       SFREE(nactdofinv);
1171       NNEW(nactdofinv,ITG,mt**nk);
1172       NNEW(nodorig,ITG,*nk);
1173       FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
1174                              ipkon,lakon,kon,ne));
1175       SFREE(nodorig);
1176       
1177       /* check whether the forced displacements changed; if so, and
1178          if the procedure is static, the first iteration has to be
1179          purely linear elastic, in order to get an equilibrium
1180          displacement field; otherwise huge (maybe nonelastic)
1181          stresses may occur, jeopardizing convergence */
1182       
1183       ilin=0;
1184       
1185       /* only for iinc=1 a linearized calculation is performed, since
1186          for iinc>1 a reasonable displacement field is predicted by using the
1187          initial velocity field at the end of the last increment */
1188       
1189       if((iinc==1)&&(*ithermal<2)){
1190           dev=0.;
1191           for(k=0;k<*nboun;++k){
1192               err=fabs(xbounact[k]-xbounini[k]);
1193               if(err>dev){dev=err;}
1194           }
1195           if(dev>1.e-5) ilin=1;
1196       }
1197       
1198       /* prediction of the kinematic vectors  */
1199       
1200       NNEW(v,double,mt**nk);
1201       
1202       if(ncont>idiscon)idiscon=ncont;
1203       prediction(uam,nmethod,&bet,&gam,&dtime,ithermal,nk,veold,accold,v,
1204                  &iinc,&idiscon,vold,nactdof,mi);
1205       
1206       NNEW(fn,double,mt**nk);
1207       NNEW(stx,double,6*mi[0]**ne);
1208       
1209       /* determining the internal forces at the start of the increment
1210          
1211          for a static calculation with increased forced displacements
1212          the linear strains are calculated corresponding to
1213          
1214          the displacements at the end of the previous increment, extrapolated
1215          if appropriate (for nondispersive media) +
1216          the forced displacements at the end of the present increment +
1217          the temperatures at the end of the present increment (this sum is
1218          v) -
1219          the displacements at the end of the previous increment (this is vold)
1220          
1221          these linear strains are converted in stresses by multiplication
1222          with the tangent element stiffness matrix and converted into nodal
1223          forces. 
1224          
1225          this boils down to the fact that the effect of forced displacements
1226          should be handled in a purely linear way at the
1227          start of a new increment, in order to speed up the convergence and
1228          (for dissipative media) guarantee smooth loading within the increment.
1229          
1230          for all other cases the nodal force calculation is based on
1231          the true stresses derived from the appropriate strain tensor taking
1232          into account the extrapolated displacements at the end of the 
1233          previous increment + the forced displacements and the temperatures
1234          at the end of the present increment */
1235       
1236       iout=-1;
1237       iperturb_sav[0]=iperturb[0];
1238       iperturb_sav[1]=iperturb[1];
1239       
1240       /* first iteration in first increment: elastic tangent */
1241       
1242       if((*nmethod!=4)&&(ilin==1)){
1243           
1244           ielas=1;
1245           
1246           iperturb[0]=-1;
1247           iperturb[1]=0;
1248           
1249           for(k=0;k<neq[1];++k){b[k]=f[k];}
1250           NNEW(inum,ITG,*nk);
1251           results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1252                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1253                   ielorien,norien,orab,ntmat_,t1ini,t1act,ithermal,
1254                   prestr,iprestr,filab,eme,emn,een,iperturb,
1255                   f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1256                   ndirboun,xbounact,nboun,ipompc,
1257                   nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1258                   &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1259                   xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1260                   &icmd, ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1261                   emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1262                   iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1263                   fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1264                   &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1265                   sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1266                   mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1267                   islavsurf,ielprop,prop);
1268           iperturb[0]=0;SFREE(inum);
1269           
1270           /* check whether any displacements or temperatures are changed
1271              in the new increment */
1272           
1273           for(k=0;k<neq[1];++k){f[k]=f[k]+b[k];}
1274           
1275       }
1276       else{
1277           
1278           NNEW(inum,ITG,*nk);
1279           results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1280                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1281                   ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1282                   prestr,iprestr,filab,eme,emn,een,iperturb,
1283                   f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1284                   ndirboun,xbounact,nboun,ipompc,
1285                   nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1286                   &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1287                   xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1288                   &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1289                   emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1290                   iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1291                   fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1292                   &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1293                   sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1294                   mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1295                   islavsurf,ielprop,prop);
1296           SFREE(inum);
1297           
1298           memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1299           
1300           if(*ithermal!=2){
1301               for(k=0;k<6*mi[0]*ne0;++k){
1302                   sti[k]=stx[k];
1303               }
1304           }
1305           
1306       }
1307       
1308       ielas=0;
1309       iout=0;
1310       
1311       SFREE(fn);SFREE(stx);SFREE(v);
1312       
1313       /***************************************************************/
1314       /* iteration counter and start of the loop over the iterations */
1315       /***************************************************************/
1316 
1317     iit=1;
1318     icntrl=0;
1319     ctrl[0]=i0ref;ctrl[1]=irref;ctrl[3]=icref;
1320     if(uncoupled){
1321         *ithermal=2;
1322         NNEW(iruc,ITG,nzs[1]-nzs[0]);
1323         for(k=0;k<nzs[1]-nzs[0];k++) iruc[k]=irow[k+nzs[0]]-neq[0];
1324     }
1325 
1326     while(icntrl==0){
1327 
1328     /*  updating the nonlinear mpc's (also affects the boundary
1329         conditions through the nonhomogeneous part of the mpc's) */
1330 
1331       if((iit!=1)||((uncoupled)&&(*ithermal==1))){
1332 
1333           printf(" iteration %" ITGFORMAT "\n\n",iit);
1334 
1335           FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
1336               xloadold,xload,
1337               xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
1338               t1old,t1,t1act,iamt1,nk,amta,
1339               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
1340               xbounold,xboun,xbounact,iamboun,nboun,
1341               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
1342               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
1343               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
1344               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
1345 
1346           for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
1347           if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,bcr,ipivr,
1348              ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,
1349              ntmat_,vold,shcon,nshcon,ipkon,kon,co,
1350              kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
1351              nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
1352              cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
1353              ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
1354              nactdog,nacteq,nodeboun,ndirboun,&network,
1355              rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
1356              ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,
1357              nset,ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,
1358              iamload,jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
1359           }
1360 
1361           if((icascade==2)||
1362              ((ncont!=0)&&(ismallsliding==0))){
1363               memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
1364               RENEW(nodempc,ITG,3*memmpcref_);
1365               for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
1366               RENEW(coefmpc,double,memmpcref_);
1367               for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
1368           }
1369 
1370           if((ncont!=0)&&(*mortar<=1)&&(ismallsliding==0)&&((iit<=8)||(*mortar==1))){
1371               neold=*ne;
1372               *ne=ne0;*nkon=nkon0;
1373               contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1374                       ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,
1375                       straight,nkon,co,vold,ielmat,cs,elcon,istep,
1376                       &iinc,&iit,ncmat_,ntmat_,&ne0,
1377                       vini,nmethod,nmpc,&mpcfree,&memmpc_,&ipompc,&labmpc,
1378                       &ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,iperturb,
1379                       ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
1380                       itiefac,areaslav,iponoels,inoels,springarea,tietol,
1381                       &reltime,imastnode,nmastnode,xmastnor,
1382                       filab,mcs,ics,&nasym,xnoels,mortar,pslavsurf,pmastsurf,
1383                       clearini,&theta);
1384 
1385               /* check whether, for a dynamic calculation, damping is involved */
1386               
1387               if(*nmethod==4){
1388                   if(*iexpl<=1){
1389                       if(idampingwithoutcontact==0){
1390                           for(i=0;i<*ne;i++){
1391                               if(ipkon[i]<0) continue;
1392                               if(*ncmat_>=5){
1393                                   if(strcmp1(&lakon[i*8],"ES")==0){
1394                                       if(strcmp1(&lakon[i*8+6],"C")==0){
1395                                           imat=ielmat[i*mi[2]];
1396                                           if(elcon[(*ncmat_+1)**ntmat_*(imat-1)+4]>0.){
1397                                               idamping=1;break;
1398                                           }
1399                                       }
1400                                   }
1401                               }
1402                           }
1403                       }
1404                   }
1405               }
1406               
1407               if(*mortar==0){
1408                  if(*ne!=neold){iflagact=1;}
1409               }else if(*mortar==1){
1410                  if(((*ne-ne0)<(neold-ne0)*0.999)||
1411                     ((*ne-ne0)>(neold-ne0)*1.001)){iflagact=1;}
1412               }
1413 
1414               printf(" Number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
1415 
1416               if(icascade<1)icascade=1;
1417           }
1418           
1419           if(*ithermal==3){
1420               for(k=0;k<*nk;++k){t1act[k]=vold[mt*k];}
1421           }
1422 
1423           FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1424                 nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1425                 &maxlenmpc,ikmpc,ilmpc,&icascade,
1426                 kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,&iit,
1427                 &idiscon,&ncont,trab,ntrans,ithermal,mi));
1428 
1429           if((icascade==2)||
1430               ((ncont!=0)&&(ismallsliding==0))){
1431 //            memmpcref_=memmpc_;mpcfreeref=mpcfree;
1432 //            RENEW(nodempcref,ITG,3*memmpc_);
1433               for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
1434 //            RENEW(coefmpcref,double,memmpc_);
1435               for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
1436           }
1437 
1438           if(icascade>0){
1439               remastruct(ipompc,&coefmpc,&nodempc,nmpc,
1440                 &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
1441                 labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
1442                 kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
1443                 neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
1444                 &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
1445                 typeboun,&cv,&cvini,&iit);
1446 
1447               /* invert nactdof */
1448               
1449               SFREE(nactdofinv);
1450               NNEW(nactdofinv,ITG,mt**nk);
1451               NNEW(nodorig,ITG,*nk);
1452               FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
1453                                      ipkon,lakon,kon,ne));
1454               SFREE(nodorig);
1455               
1456               NNEW(v,double,mt**nk);
1457               NNEW(stx,double,6*mi[0]**ne);
1458               NNEW(fn,double,mt**nk);
1459       
1460               memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1461               iout=-1;
1462               
1463               NNEW(inum,ITG,*nk);
1464               results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1465                 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1466                 ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1467                 prestr,iprestr,filab,eme,emn,een,iperturb,
1468                 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1469                 ndirboun,xbounact,nboun,ipompc,
1470                 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1471                 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1472                 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1473                 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1474                 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1475                 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1476                 nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1477                 &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1478                 sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1479                 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1480                 islavsurf,ielprop,prop);
1481 
1482               /*for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}*/
1483               
1484               SFREE(v);SFREE(stx);SFREE(fn);SFREE(inum);
1485               iout=0;
1486               
1487           }else{
1488 
1489               /*for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}*/
1490           }
1491       }
1492       
1493       if(*iexpl<=1){
1494 
1495         /* calculating the local stiffness matrix and external loading */
1496 
1497         NNEW(ad,double,neq[1]);
1498         NNEW(au,double,nzs[1]);
1499 
1500         FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounact,nboun,
1501                   ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1502                   nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1503                   nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
1504                   nmethod,ikmpc,ilmpc,ikboun,ilboun,
1505                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1506                   ielmat,ielorien,norien,orab,ntmat_,
1507                   t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1508                   nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1509                   xstiff,npmat_,&dtime,matname,mi,
1510                   ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
1511                   physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
1512                   &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
1513                   xstateini,xstate,thicke,integerglob,doubleglob,
1514                   tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
1515                   pmastsurf,mortar,clearini,ielprop,prop));
1516 
1517         if(nasym==1){
1518             RENEW(au,double,2*nzs[1]);
1519             if(*nmethod==4) RENEW(aub,double,2*nzs[1]);
1520             symmetryflag=2;
1521             inputformat=1;
1522 
1523             FORTRAN(mafillsmas,(co,nk,kon,ipkon,lakon,ne,nodeboun,
1524                   ndirboun,xbounact,nboun,
1525                   ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1526                   nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1527                   nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
1528                   nmethod,ikmpc,ilmpc,ikboun,ilboun,
1529                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1530                   ielmat,ielorien,norien,orab,ntmat_,
1531                   t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1532                   nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1533                   xstiff,npmat_,&dtime,matname,mi,
1534                   ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
1535                   physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
1536                   &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
1537                   xstateini,xstate,thicke,
1538                   integerglob,doubleglob,tieset,istartset,iendset,
1539                   ialset,ntie,&nasym,pslavsurf,pmastsurf,mortar,clearini,
1540                   ielprop,prop));
1541         }
1542             
1543 
1544         iperturb[0]=iperturb_sav[0];
1545         iperturb[1]=iperturb_sav[1];
1546 
1547       }else{
1548 
1549         /* calculating the external loading 
1550 
1551            This is only done once per increment. In reality, the
1552            external loading is a function of vold (specifically,
1553            the body forces and surface loading). This effect is
1554            neglected, since the increment size in dynamic explicit
1555            calculations is usually small */
1556 
1557           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1558                   ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1559                   nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1560                   nbody,cgr,fext,nactdof,&neq[1],
1561                   nmethod,ikmpc,ilmpc,
1562                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1563                   ielmat,ielorien,norien,orab,ntmat_,
1564                   t0,t1act,ithermal,iprestr,vold,iperturb,
1565                   iexpl,plicon,nplicon,plkcon,nplkcon,
1566                   npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1567                   xbodyold,&reltime,veold,matname,mi,ikactmech,
1568                   &nactmech,ielprop,prop));
1569       }
1570       
1571 /*      for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}
1572       for(k=0;k<neq[1];++k){printf("fext=%" ITGFORMAT ",%f\n",k,fext[k]);}
1573       for(k=0;k<neq[1];++k){printf("ad=%" ITGFORMAT ",%f\n",k,ad[k]);}
1574       for(k=0;k<nzs[1];++k){printf("au=%" ITGFORMAT ",%f\n",k,au[k]);}*/
1575 
1576       /* calculating the damping matrix for implicit dynamic
1577          calculations */
1578 
1579       if(idamping==1){
1580 
1581           /* Rayleigh damping */
1582 
1583           NNEW(adc,double,neq[1]);
1584           for(k=0;k<neq[0];k++){adc[k]=alpham*adb[k]+betam*ad[k];}
1585           if(nasym==0){
1586               NNEW(auc,double,nzs[1]);
1587               for(k=0;k<nzs[0];k++){auc[k]=alpham*aub[k]+betam*au[k];}
1588           }else{
1589               NNEW(auc,double,2*nzs[1]);
1590               for(k=0;k<2*nzs[0];k++){auc[k]=alpham*aub[k]+betam*au[k];}
1591           }
1592  
1593           /* dashpots and contact damping */
1594 
1595           FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
1596               ndirboun,xbounact,nboun,
1597               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1598               nforc,nelemload,sideload,xloadact,nload,xbodyact,
1599               ipobody,nbody,cgr,
1600               adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
1601               ikmpc,ilmpc,ikboun,ilboun,
1602               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1603               ielorien,norien,orab,ntmat_,
1604               t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1605               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1606               xstiff,npmat_,&dtime,matname,mi,ncmat_,
1607               ttime,&time,istep,&iinc,ibody,clearini,mortar,springarea,
1608               pslavsurf,pmastsurf,&reltime,&nasym));
1609       }
1610 
1611       /* calculating the residual */
1612 
1613       calcresidual(nmethod,neq,b,fext,f,iexpl,nactdof,aux2,vold,
1614          vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1615          islavnode,nslavnode,mortar,ntie,f_cm,f_cs,mi,
1616          nzs,&nasym,&idamping,veold,adc,auc,cvini,cv);
1617           
1618       newstep=0;
1619       
1620       if(*nmethod==0){
1621           
1622           /* error occurred in mafill: storing the geometry in frd format */
1623           
1624           *nmethod=0;
1625           ++*kode;
1626           NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
1627           if(strcmp1(&filab[1044],"ZZS")==0){
1628               NNEW(neigh,ITG,40**ne);
1629               NNEW(ipneigh,ITG,*nk);
1630           }
1631           
1632           ptime=*ttime+time;
1633           frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
1634                   kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1635                   nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1636                   ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1637                   mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1638                   cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1639                   thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
1640 
1641           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);} 
1642      #ifdef COMPANY
1643           FORTRAN(uout,(v,mi,ithermal));
1644 #endif
1645           SFREE(inum);FORTRAN(stop,());
1646           
1647       }
1648       
1649       /* implicit step (static or dynamic */
1650       
1651       if(*iexpl<=1){
1652           if((*nmethod==4)&&(*mortar<2)){
1653               
1654               /* mechanical part */
1655               
1656               if(*ithermal!=2){
1657                   scal1=bet*dtime*dtime*(1.+*alpha);
1658                   for(k=0;k<neq[0];++k){
1659                       ad[k]=adb[k]+scal1*ad[k];
1660                   }
1661                   for(k=0;k<nzs[0];++k){
1662                       au[k]=aub[k]+scal1*au[k];
1663                   }
1664                   
1665                   /* upper triangle of asymmetric matrix */
1666                   
1667                   if(nasym>0){
1668                       for(k=nzs[2];k<nzs[2]+nzs[0];++k){
1669                           au[k]=aub[k]+scal1*au[k];
1670                       }
1671                   }
1672 
1673                   /* damping */
1674                   
1675                   if(idamping==1){
1676                       scal1=gam*dtime*(1.+*alpha);
1677                       for(k=0;k<neq[0];++k){
1678                           ad[k]+=scal1*adc[k];
1679                       }
1680                       for(k=0;k<nzs[0];++k){
1681                           au[k]+=scal1*auc[k];
1682                       }
1683                       
1684                       /* upper triangle of asymmetric matrix */
1685                       
1686                       if(nasym>0){
1687                           for(k=nzs[2];k<nzs[2]+nzs[0];++k){
1688                               au[k]+=scal1*auc[k];
1689                           }
1690                       }
1691                   }
1692 
1693               }
1694               
1695               /* thermal part */
1696               
1697               if(*ithermal>1){
1698                   for(k=neq[0];k<neq[1];++k){
1699                       ad[k]=adb[k]/dtime+ad[k];
1700                   }
1701                   for(k=nzs[0];k<nzs[1];++k){
1702                       au[k]=aub[k]/dtime+au[k];
1703                   }
1704                   
1705                   /* upper triangle of asymmetric matrix */
1706                   
1707                   if(nasym>0){
1708                       for(k=nzs[2]+nzs[0];k<nzs[2]+nzs[1];++k){
1709                           au[k]=aub[k]/dtime+au[k];
1710                       }
1711                   }
1712               }
1713           }
1714           
1715           if(*isolver==0){
1716 #ifdef SPOOLES
1717               if(*ithermal<2){
1718                   spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1719                           &symmetryflag,&inputformat,&nzs[2]);
1720               }else if((*ithermal==2)&&(uncoupled)){
1721                   n1=neq[1]-neq[0];
1722                   n2=nzs[1]-nzs[0];
1723                   spooles(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1724                           &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1725                           &n1,&n2,&symmetryflag,&inputformat,&nzs[2]);
1726               }else{
1727                   spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1728                           &symmetryflag,&inputformat,&nzs[2]);
1729               }
1730 #else
1731               printf(" *ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
1732               FORTRAN(stop,());
1733 #endif
1734           }
1735           else if((*isolver==2)||(*isolver==3)){
1736               if(nasym>0){
1737                   printf(" *ERROR in nonlingeo: the iterative solver cannot be used for asymmetric matrices\n\n");
1738                   FORTRAN(stop,());
1739               }
1740               preiter(ad,&au,b,&icol,&irow,&neq[1],&nzs[1],isolver,iperturb);
1741           }
1742           else if(*isolver==4){
1743 #ifdef SGI
1744               if(nasym>0){
1745                   printf(" *ERROR in nonlingeo: the SGI solver cannot be used for asymmetric matrices\n\n");
1746                   FORTRAN(stop,());
1747               }
1748               token=1;
1749               if(*ithermal<2){
1750                   sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
1751               }else if((*ithermal==2)&&(uncoupled)){
1752                   n1=neq[1]-neq[0];
1753                   n2=nzs[1]-nzs[0];
1754                   sgi_main(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1755                            &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1756                            &n1,&n2,token);
1757               }else{
1758                   sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],token);
1759               }
1760 #else
1761               printf(" *ERROR in nonlingeo: the SGI library is not linked\n\n");
1762               FORTRAN(stop,());
1763 #endif
1764           }
1765           else if(*isolver==5){
1766               if(nasym>0){
1767                   printf(" *ERROR in nonlingeo: the TAUCS solver cannot be used for asymmetric matrices\n\n");
1768                   FORTRAN(stop,());
1769               }
1770 #ifdef TAUCS
1771               tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[1],&nzs[1]);
1772 #else
1773               printf(" *ERROR in nonlingeo: the TAUCS library is not linked\n\n");
1774               FORTRAN(stop,());
1775 #endif
1776           }
1777           else if(*isolver==7){
1778 #ifdef PARDISO
1779               if(*ithermal<2){
1780                   pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1781                                &symmetryflag,&inputformat,jq,&nzs[2]);
1782               }else if((*ithermal==2)&&(uncoupled)){
1783                   n1=neq[1]-neq[0];
1784                   n2=nzs[1]-nzs[0];
1785                   pardiso_main(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1786                                &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1787                                &n1,&n2,&symmetryflag,&inputformat,jq,&nzs[2]);
1788               }else{
1789                   pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1790                                &symmetryflag,&inputformat,jq,&nzs[2]);
1791               }
1792 #else
1793               printf(" *ERROR in nonlingeo: the PARDISO library is not linked\n\n");
1794               FORTRAN(stop,());
1795 #endif
1796           }
1797           
1798           if(*mortar<=1){SFREE(ad);SFREE(au);} 
1799       }
1800       
1801       /* explicit dynamic step */
1802       
1803       else{
1804           if(*ithermal!=2){
1805               for(k=0;k<neq[0];++k){
1806                   b[k]=b[k]/adb[k];
1807               }
1808           }
1809           if(*ithermal>1){
1810               for(k=neq[0];k<neq[1];++k){
1811                   b[k]=b[k]*dtime/adb[k];
1812               }
1813           }
1814       }
1815 
1816       /* calculating the displacements, stresses and forces */
1817       
1818       NNEW(v,double,mt**nk);
1819       memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1820       
1821       NNEW(stx,double,6*mi[0]**ne);
1822       NNEW(fn,double,mt**nk);
1823       
1824       NNEW(inum,ITG,*nk);
1825       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1826               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1827               ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1828               prestr,iprestr,filab,eme,emn,een,iperturb,
1829               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1830               ndirboun,xbounact,nboun,ipompc,
1831               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1832               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1833               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1834               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1835               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1836               iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1837               fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1838               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1839               sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1840               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1841               islavsurf,ielprop,prop);
1842       SFREE(inum);
1843 
1844       if(*ithermal!=2){
1845           if(cam[0]>uam[0]){uam[0]=cam[0];}      
1846           if(qau<1.e-10){
1847               if(qa[0]>ea*qam[0]){qam[0]=(qamold[0]*jnz+qa[0])/(jnz+1);}
1848               else {qam[0]=qamold[0];}
1849           }
1850       }
1851       if(*ithermal>1){
1852           if(cam[1]>uam[1]){uam[1]=cam[1];}      
1853           if(qau<1.e-10){
1854               if(qa[1]>ea*qam[1]){qam[1]=(qamold[1]*jnz+qa[1])/(jnz+1);}
1855               else {qam[1]=qamold[1];}
1856           }
1857       }
1858 
1859       memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1860       if(*ithermal!=2){
1861           for(k=0;k<6*mi[0]*ne0;++k){
1862               sti[k]=stx[k];
1863           }
1864       }
1865       
1866       /* calculating the ratio of the smallest to largest pressure
1867          for face-to-face contact
1868          only done at the end of a step */
1869 
1870       if((*mortar==1)&&(1.-theta-dtheta<=1.e-6)){
1871           FORTRAN(negativepressure,(&ne0,ne,mi,stx,&pressureratio));
1872       }else{pressureratio=0.;}
1873 
1874       SFREE(v);SFREE(stx);SFREE(fn);
1875       
1876       /* calculating the residual */
1877       
1878       calcresidual(nmethod,neq,b,fext,f,iexpl,nactdof,aux2,vold,
1879          vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1880          islavnode,nslavnode,mortar,ntie,f_cm,f_cs,mi,
1881          nzs,&nasym,&idamping,veold,adc,auc,cvini,cv);
1882 
1883       if(idamping==1){SFREE(adc);SFREE(auc);}
1884       
1885       /* calculating the maximum residual */
1886 
1887       for(k=0;k<2;++k){
1888           ram2[k]=ram1[k];
1889           ram1[k]=ram[k];
1890           ram[k]=0.;
1891       }
1892       if(*ithermal!=2){
1893           for(k=0;k<neq[0];++k){
1894               err=fabs(b[k]);
1895               if(err>ram[0]){ram[0]=err;ram[2]=k+0.5;}
1896           }
1897       }
1898       if(*ithermal>1){
1899           for(k=neq[0];k<neq[1];++k){
1900               err=fabs(b[k]);
1901               if(err>ram[1]){ram[1]=err;ram[3]=k+0.5;}
1902           }
1903       }
1904 
1905   /*   Divergence criteria for face-to-face penalty is different */
1906 
1907       if(*mortar==1){
1908           for(k=4;k<8;++k){
1909               ram2[k]=ram1[k];
1910               ram1[k]=ram[k];
1911           } 
1912           ram[4]=ram[0]+ram1[0];
1913           if((iflagact==0)&&(iit>1)){
1914               ram[5]=1.5;
1915           }else{ram[5]=0.5;}
1916           ram[6]=abs((neold-ne0)-(*ne-ne0))+0.5;
1917           if(iit>3){
1918               if((ram[6]>=ram1[6])&&(ram[6]>=ram2[6])){
1919                   ram[7]=1.5;
1920               }else{ram[7]=0.5;}
1921           }
1922           
1923 /*        for(k=2;k<4;++k){
1924               ram2[k]=ram1[k];
1925               ram1[k]=ram[k+2];
1926           } 
1927           ram[4]=ram[0]+ram1[0];
1928           
1929           ram[5]=abs((neold-ne0)-(*ne-ne0));
1930           
1931           if(iit==2){ram[6]=ram[4];}
1932           if((iit>=3)&&(ram[6]>ram[4])){ram[6]=ram[4];}*/
1933       }
1934       
1935       /* next line is inserted to cope with stress-less
1936          temperature calculations */
1937       
1938       if(*ithermal!=2){
1939           if(ram[0]<1.e-6) ram[0]=0.;      
1940           printf(" average force= %f\n",qa[0]);
1941           printf(" time avg. forc= %f\n",qam[0]);
1942           if((ITG)((double)nactdofinv[(ITG)ram[2]]/mt)+1==0){
1943               printf(" largest residual force= %f\n",
1944                  ram[0]);
1945           }else{
1946               inode=(ITG)((double)nactdofinv[(ITG)ram[2]]/mt)+1;
1947               idir=nactdofinv[(ITG)ram[2]]-mt*(inode-1);
1948               printf(" largest residual force= %f in node %" ITGFORMAT " and dof %" ITGFORMAT "\n",
1949                      ram[0],inode,idir);
1950           }
1951           printf(" largest increment of disp= %e\n",uam[0]);
1952           if((ITG)cam[3]==0){
1953               printf(" largest correction to disp= %e\n\n",
1954                  cam[0]);
1955           }else{
1956               inode=(ITG)((double)nactdofinv[(ITG)cam[3]]/mt)+1;
1957               idir=nactdofinv[(ITG)cam[3]]-mt*(inode-1);
1958               printf(" largest correction to disp= %e in node %" ITGFORMAT " and dof %" ITGFORMAT "\n\n",cam[0],inode,idir);
1959           }
1960       }
1961       if(*ithermal>1){
1962           if(ram[1]<1.e-6) ram[1]=0.;      
1963           printf(" average flux= %f\n",qa[1]);
1964           printf(" time avg. flux= %f\n",qam[1]);
1965           if((ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1==0){
1966               printf(" largest residual flux= %f\n",
1967                  ram[1]);
1968           }else{
1969               inode=(ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1;
1970               idir=nactdofinv[(ITG)ram[3]]-mt*(inode-1);
1971               printf(" largest residual flux= %f in node %" ITGFORMAT " and dof %" ITGFORMAT "\n",ram[1],inode,idir);
1972           }
1973           printf(" largest increment of temp= %e\n",uam[1]);
1974           if((ITG)cam[4]==0){
1975               printf(" largest correction to temp= %e\n\n",
1976                  cam[1]);
1977           }else{
1978               inode=(ITG)((double)nactdofinv[(ITG)cam[4]]/mt)+1;
1979               idir=nactdofinv[(ITG)cam[4]]-mt*(inode-1);
1980               printf(" largest correction to temp= %e in node %" ITGFORMAT " and dof %" ITGFORMAT "\n\n",cam[1],inode,idir);
1981           }
1982       }
1983       fflush(stdout);
1984 
1985       FORTRAN(writecvg,(istep,&iinc,&iit,ne,&ne0,ram,qam,cam,uam,
1986                         ithermal));
1987 
1988       checkconvergence(co,nk,kon,ipkon,lakon,ne,stn,nmethod, 
1989           kode,filab,een,t1act,&time,epn,ielmat,matname,enern, 
1990           xstaten,nstate_,istep,&iinc,iperturb,ener,mi,output,
1991           ithermal,qfn,&mode,&noddiam,trab,inotr,ntrans,orab,
1992           ielorien,norien,description,sti,&icutb,&iit,&dtime,qa,
1993           vold,qam,ram1,ram2,ram,cam,uam,&ntg,ttime,&icntrl,
1994           &theta,&dtheta,veold,vini,idrct,tper,&istab,tmax, 
1995           nactdof,b,tmin,ctrl,amta,namta,itpamp,inext,&dthetaref,
1996           &itp,&jprint,jout,&uncoupled,t1,&iitterm,nelemload,
1997           nload,nodeboun,nboun,itg,ndirboun,&deltmx,&iflagact,
1998           set,nset,istartset,iendset,ialset,emn,thicke,jobnamec,
1999           mortar,nmat,ielprop,prop);
2000 
2001     }
2002 
2003     /*********************************************************/
2004     /*   end of the iteration loop                          */
2005     /*********************************************************/
2006 
2007     /* icutb=0 means that the iterations in the increment converged,
2008        icutb!=0 indicates that the increment has to be reiterated with
2009                 another increment size (dtheta) */
2010    
2011     if(uncoupled){
2012         SFREE(iruc);
2013     }
2014 
2015     if(((qa[0]>ea*qam[0])||(qa[1]>ea*qam[1]))&&(icutb==0)){jnz++;}
2016     iit=0;
2017 
2018     if(icutb!=0){
2019       memcpy(&vold[0],&vini[0],sizeof(double)*mt**nk);
2020 
2021       for(k=0;k<*nboun;++k){xbounact[k]=xbounini[k];}
2022       if((*ithermal==1)||(*ithermal>=3)){
2023         for(k=0;k<*nk;++k){t1act[k]=t1ini[k];}
2024       }
2025       for(k=0;k<neq[1];++k){
2026           f[k]=fini[k];
2027       }
2028       if(*nmethod==4){
2029         for(k=0;k<mt**nk;++k){
2030           veold[k]=veini[k];
2031           accold[k]=accini[k];
2032         }
2033         for(k=0;k<neq[1];++k){
2034 //        f[k]=fini[k];
2035           fext[k]=fextini[k];
2036           cv[k]=cvini[k];
2037         }
2038       }
2039       if(*ithermal!=2){
2040           for(k=0;k<6*mi[0]*ne0;++k){
2041               sti[k]=stiini[k];
2042               eme[k]=emeini[k];
2043           }
2044       }
2045       if(*nener==1)
2046           for(k=0;k<mi[0]*ne0;++k){ener[k]=enerini[k];}
2047 
2048       for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
2049           xstate[k]=xstateini[k];
2050       }   
2051 
2052       qam[0]=qamold[0];
2053       qam[1]=qamold[1];
2054     }
2055     
2056     /* face-to-face penalty */
2057 
2058     if((*mortar==1)&&(icutb==0)){
2059         
2060         ntrimax=0;
2061         for(i=0;i<*ntie;i++){       
2062             if(itietri[2*i+1]-itietri[2*i]+1>ntrimax)           
2063                 ntrimax=itietri[2*i+1]-itietri[2*i]+1;          
2064         }
2065         NNEW(xo,double,ntrimax);            
2066         NNEW(yo,double,ntrimax);            
2067         NNEW(zo,double,ntrimax);            
2068         NNEW(x,double,ntrimax);     
2069         NNEW(y,double,ntrimax);     
2070         NNEW(z,double,ntrimax);    
2071         NNEW(nx,ITG,ntrimax);      
2072         NNEW(ny,ITG,ntrimax);       
2073         NNEW(nz,ITG,ntrimax);
2074       
2075         /*  Determination of active nodes (islavact) */
2076       
2077         FORTRAN(islavactive,(tieset,ntie,itietri,cg,straight,
2078                    co,vold,xo,yo,zo,x,y,z,nx,ny,nz,mi,
2079                    imastop,nslavnode,islavnode,islavact));
2080 
2081         SFREE(xo);SFREE(yo);SFREE(zo);SFREE(x);SFREE(y);SFREE(z);SFREE(nx);         
2082         SFREE(ny);SFREE(nz);
2083 
2084         if(negpres==0){
2085             if((*mortar==1)&&(1.-theta-dtheta<=1.e-6)&&(itruecontact==1)){
2086                 printf(" pressure ratio (smallest/largest pressure over all contact areas) =%e\n\n",pressureratio);
2087                 if(pressureratio<-0.05){
2088                     printf(" zero-size increment is appended\n\n");
2089                     negpres=1;dtheta=0.;
2090                 }
2091             }
2092         }else{negpres=0;}
2093 
2094     }
2095 
2096     /* output */
2097 
2098     if((jout[0]==jprint)&&(icutb==0)){
2099 
2100       jprint=0;
2101 
2102       /* calculating the displacements and the stresses and storing */
2103       /* the results in frd format  */
2104         
2105       NNEW(v,double,mt**nk);
2106       NNEW(fn,double,mt**nk);
2107       NNEW(stn,double,6**nk);
2108       if(*ithermal>1) NNEW(qfn,double,3**nk);
2109       NNEW(inum,ITG,*nk);
2110       NNEW(stx,double,6*mi[0]**ne);
2111       
2112       if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
2113       if(strcmp1(&filab[435],"PEEQ")==0) NNEW(epn,double,*nk);
2114       if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2115       if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
2116       if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
2117       if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
2118 
2119       memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
2120 
2121       iout=2;
2122       icmd=3;
2123       
2124 #ifdef COMPANY
2125       FORTRAN(uinit,());
2126 #endif
2127       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
2128               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2129               ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
2130               prestr,iprestr,filab,eme,emn,een,iperturb,
2131               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
2132               ndirboun,xbounact,nboun,ipompc,
2133               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
2134               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
2135               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
2136               ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
2137               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
2138               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
2139               nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
2140               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
2141               sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2142               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
2143               islavsurf,ielprop,prop);
2144       
2145       memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
2146 
2147       iout=0;
2148       if(*iexpl<=1) icmd=0;
2149 //      FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
2150 //      for(k=0;k<ntg;k++)if(inum[itg[k]-1]>0){inum[itg[k]-1]*=-1;}
2151       
2152       ++*kode;
2153       if(*mcs!=0){
2154         ptime=*ttime+time;
2155         frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
2156                t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
2157                nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
2158                ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
2159                norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,&ne0,
2160                cdn,mortar,nmat);
2161 #ifdef COMPANY
2162         FORTRAN(uout,(v,mi,ithermal));
2163 #endif
2164       }
2165       else{
2166           if(strcmp1(&filab[1044],"ZZS")==0){
2167               NNEW(neigh,ITG,40**ne);
2168               NNEW(ipneigh,ITG,*nk);
2169           }
2170 
2171           ptime=*ttime+time;
2172           frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
2173             kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2174             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2175             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2176             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
2177             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2178             thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
2179 
2180           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2181 #ifdef COMPANY
2182           FORTRAN(uout,(v,mi,ithermal));
2183 #endif
2184       }
2185       
2186       SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
2187       if(*ithermal>1){SFREE(qfn);}
2188       
2189       if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
2190       if(strcmp1(&filab[435],"PEEQ")==0) SFREE(epn);
2191       if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2192       if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
2193       if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
2194       if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
2195    }
2196     
2197   }
2198 
2199   /*********************************************************/
2200   /*   end of the increment loop                          */
2201   /*********************************************************/
2202 
2203   if(jprint!=0){
2204 
2205   /* calculating the displacements and the stresses and storing  
2206      the results in frd format */
2207   
2208     NNEW(v,double,mt**nk);
2209     NNEW(fn,double,mt**nk);
2210     NNEW(stn,double,6**nk);
2211     if(*ithermal>1) NNEW(qfn,double,3**nk);
2212     NNEW(inum,ITG,*nk);
2213     NNEW(stx,double,6*mi[0]**ne);
2214   
2215     if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
2216     if(strcmp1(&filab[435],"PEEQ")==0) NNEW(epn,double,*nk);
2217     if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2218     if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
2219     if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
2220     if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
2221     
2222     memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
2223     iout=2;
2224     icmd=3;
2225 
2226 #ifdef COMPANY
2227     FORTRAN(uinit,());
2228 #endif
2229     results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
2230             elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2231             ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
2232             prestr,iprestr,filab,eme,emn,een,iperturb,
2233             f,fn,nactdof,&iout,qa,vold,b,nodeboun,
2234             ndirboun,xbounact,nboun,ipompc,
2235             nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
2236             &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
2237             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
2238             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
2239             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
2240             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
2241             nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
2242             &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
2243             sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2244             mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
2245             islavsurf,ielprop,prop);
2246     
2247     memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
2248 
2249     iout=0;
2250     if(*iexpl<=1) icmd=0;
2251 //    FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
2252 //    for(k=0;k<ntg;k++)if(inum[itg[k]-1]>0){inum[itg[k]-1]*=-1;}
2253     
2254     ++*kode;
2255     if(*mcs>0){
2256         ptime=*ttime+time;
2257       frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
2258              t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
2259              nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
2260              ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
2261              norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,&ne0,
2262              cdn,mortar,nmat);
2263 #ifdef COMPANY
2264       FORTRAN(uout,(v,mi,ithermal));
2265 #endif
2266 
2267     }
2268     else{
2269         if(strcmp1(&filab[1044],"ZZS")==0){
2270             NNEW(neigh,ITG,40**ne);
2271             NNEW(ipneigh,ITG,*nk);
2272         }
2273 
2274         ptime=*ttime+time;
2275         frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
2276             kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2277             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2278             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2279             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
2280             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2281             thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
2282 
2283         if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2284 #ifdef COMPANY
2285         FORTRAN(uout,(v,mi,ithermal));
2286 #endif
2287     }
2288     
2289     SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
2290     if(*ithermal>1){SFREE(qfn);}
2291     
2292     if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
2293     if(strcmp1(&filab[435],"PEEQ")==0) SFREE(epn);
2294     if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2295     if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
2296     if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
2297     if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
2298 
2299   }
2300 
2301   /* setting the velocity to zero at the end of a quasistatic or stationary
2302      step */
2303 
2304   if(abs(*nmethod)==1){
2305     for(k=0;k<mt**nk;++k){veold[k]=0.;}
2306   }
2307 
2308   /* updating the loading at the end of the step; 
2309      important in case the amplitude at the end of the step
2310      is not equal to one */
2311 
2312   for(k=0;k<*nboun;++k){
2313 
2314       /* thermal boundary conditions are updated only if the
2315          step was thermal or thermomechanical */
2316 
2317       if(ndirboun[k]==0){
2318           if(*ithermal<2) continue;
2319 
2320           /* mechanical boundary conditions are updated only
2321              if the step was not thermal or the node is a
2322              network node */
2323 
2324       }else if((ndirboun[k]>0)&&(ndirboun[k]<4)){
2325           node=nodeboun[k];
2326           FORTRAN(nident,(itg,&node,&ntg,&id));
2327           networknode=0;
2328           if(id>0){
2329               if(itg[id-1]==node) networknode=1;
2330           }
2331           if((*ithermal==2)&&(networknode==0)) continue;
2332       }
2333       xbounold[k]=xbounact[k];
2334   }
2335   for(k=0;k<*nforc;++k){xforcold[k]=xforcact[k];}
2336   for(k=0;k<2**nload;++k){xloadold[k]=xloadact[k];}
2337   for(k=0;k<7**nbody;k=k+7){xbodyold[k]=xbodyact[k];}
2338   if(*ithermal==1){
2339     for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
2340     for(k=0;k<*nk;++k){vold[mt*k]=t1act[k];}
2341   }
2342   else if(*ithermal>1){
2343     for(k=0;k<*nk;++k){t1[k]=vold[mt*k];}
2344     if(*ithermal>=3){
2345         for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
2346     }
2347   }
2348 
2349   qaold[0]=qa[0];
2350   qaold[1]=qa[1];
2351 
2352   SFREE(f);SFREE(b);
2353   SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(xbodyact);
2354   if(*nbody>0) SFREE(ipobody);if(inewton==1){SFREE(cgr);}
2355   SFREE(fext);SFREE(ampli);SFREE(xbounini);SFREE(xstiff);
2356   if((*ithermal==1)||(*ithermal>=3)){SFREE(t1act);SFREE(t1ini);}
2357 
2358   if(*ithermal>1){
2359       SFREE(itg);SFREE(ieg);SFREE(kontri);SFREE(nloadtr);
2360       SFREE(nactdog);SFREE(nacteq);SFREE(ineighe);
2361       SFREE(tarea);SFREE(tenv);SFREE(fenv);SFREE(qfx);
2362       SFREE(erad);SFREE(ac);SFREE(bc);SFREE(ipiv);
2363       SFREE(bcr);SFREE(ipivr);SFREE(adview);SFREE(auview);SFREE(adrad);
2364       SFREE(aurad);SFREE(irowrad);SFREE(jqrad);SFREE(icolrad);
2365       if((*mcs>0)&&(ntr>0)){SFREE(inocs);}
2366   }
2367 
2368   if(icfd==1){
2369       SFREE(neifa);SFREE(neiel);SFREE(neij);SFREE(ielfa);SFREE(ifaext);
2370       SFREE(vel);SFREE(vfa);SFREE(nactdoh);SFREE(nactdohinv);
2371       SFREE(ipkonf);SFREE(lakonf);SFREE(ielmatf);
2372       if(*norien>0) SFREE(ielorienf);
2373   }
2374 
2375   SFREE(fini);
2376   if(*nmethod==4){
2377     SFREE(aux2);SFREE(fextini);SFREE(veini);SFREE(accini);
2378     SFREE(adb);SFREE(aub);SFREE(cvini);SFREE(cv);
2379   }
2380   SFREE(eei);SFREE(stiini);SFREE(emeini);
2381   if(*nener==1)SFREE(enerini);
2382   if(*nstate_!=0){SFREE(xstateini);}
2383 
2384   SFREE(aux);SFREE(iaux);SFREE(vini);
2385 
2386   if((icascade==2)||(ncont!=0)){
2387       memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
2388       RENEW(nodempc,ITG,3*memmpcref_);
2389       for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
2390       RENEW(coefmpc,double,memmpcref_);
2391       for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
2392       SFREE(nodempcref);SFREE(coefmpcref);
2393   }
2394 
2395   if(ncont!=0){
2396       *ne=ne0;*nkon=nkon0;
2397       if(*nener==1){
2398         RENEW(ener,double,mi[0]**ne*2);
2399       }
2400       RENEW(ipkon,ITG,*ne);
2401       RENEW(lakon,char,8**ne);
2402       RENEW(kon,ITG,*nkon);
2403       if(*norien>0){
2404           RENEW(ielorien,ITG,mi[2]**ne);
2405       }
2406       RENEW(ielmat,ITG,mi[2]**ne);
2407 
2408       SFREE(cg);SFREE(straight);
2409       SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
2410       SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
2411       SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
2412       SFREE(springarea);SFREE(xmastnor);
2413 
2414       if(*mortar==0){
2415           SFREE(areaslav);
2416       }else if(*mortar==1){
2417           SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
2418           SFREE(islavact);
2419       }
2420   }
2421 
2422   /* reset icascade */
2423 
2424   if(icascade==1){icascade=0;}
2425 
2426   mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
2427   mpcinfo[3]=maxlenmpc;
2428 
2429   if(iglob==1){SFREE(integerglob);SFREE(doubleglob);}
2430 
2431   *icolp=icol;*irowp=irow;*cop=co;*voldp=vold;
2432 
2433   *ipompcp=ipompc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2434   *fmpcp=fmpc;*nodempcp=nodempc;*coefmpcp=coefmpc;
2435 
2436   *ipkonp=ipkon;*lakonp=lakon;*konp=kon;*ielorienp=ielorien;
2437   *ielmatp=ielmat;*enerp=ener;*xstatep=xstate;
2438 
2439   *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
2440 
2441   (*tmin)*=(*tper);
2442   (*tmax)*=(*tper);
2443 
2444   SFREE(nactdofinv);
2445 
2446   (*ttime)+=(*tper);
2447   
2448   return;
2449 }

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