root/src/checkconvergence.c

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

DEFINITIONS

This source file includes following definitions.
  1. checkconvergence

   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 
  32 
  33 void checkconvergence(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon,
  34           ITG *ne, double *stn, ITG *nmethod, 
  35           ITG *kode, char *filab, double *een, double *t1act,
  36           double *time, double *epn,ITG *ielmat,char *matname,
  37           double *enern, double *xstaten, ITG *nstate_, ITG *istep,
  38           ITG *iinc, ITG *iperturb, double *ener, ITG *mi, char *output,
  39           ITG *ithermal, double *qfn, ITG *mode, ITG *noddiam, double *trab,
  40           ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien,
  41           char *description,double *sti,
  42           ITG *icutb, ITG *iit, double *dtime, double *qa, double *vold,
  43           double *qam, double *ram1, double *ram2, double *ram,
  44           double *cam, double *uam, ITG *ntg, double *ttime,
  45           ITG *icntrl, double *theta, double *dtheta, double *veold,
  46           double *vini, ITG *idrct, double *tper,ITG *istab, double *tmax, 
  47           ITG *nactdof, double *b, double *tmin, double *ctrl, double *amta,
  48           ITG *namta, ITG *itpamp, ITG *inext, double *dthetaref, ITG *itp,
  49           ITG *jprint, ITG *jout, ITG *uncoupled, double *t1, ITG *iitterm,
  50           ITG *nelemload, ITG *nload, ITG *nodeboun, ITG *nboun, ITG *itg,
  51           ITG *ndirboun, double *deltmx, ITG *iflagact,char *set,ITG *nset,
  52           ITG *istartset,ITG *iendset,ITG *ialset, double *emn, double *thicke,
  53           char *jobnamec,ITG *mortar,ITG *nmat,ITG *ielprop,double *prop){
  54 
  55     ITG i0,ir,ip,ic,il,ig,ia,iest,iest1=0,iest2=0,iconvergence,idivergence,
  56         ngraph=1,k,*ipneigh=NULL,*neigh=NULL,*inum=NULL,id,istart,iend,inew,
  57         i,j,mt=mi[1]+1,iexceed;
  58 
  59     double df,dc,db,dd,ran,can,rap,ea,cae,ral,da,*vr=NULL,*vi=NULL,*stnr=NULL,
  60         *stni=NULL,*vmax=NULL,*stnmax=NULL,*cs=NULL,c1[2],c2[2],reftime,
  61         *fn=NULL,*eenmax=NULL,*fnr=NULL,*fni=NULL,*qfx=NULL,*cdn=NULL,
  62         *cdnr=NULL,*cdni=NULL;
  63 
  64     /* next lines are active if the number of contact elements was
  65        changed in the present increment */
  66 
  67     if ((*iflagact==1)&&(*mortar==0)){
  68         if(ctrl[0]<*iit+4)ctrl[0]=*iit+4;
  69         if(ctrl[1]<*iit+8)ctrl[1]=*iit+8;
  70         ctrl[3]+=1;
  71     }
  72         
  73     i0=ctrl[0];ir=ctrl[1];ip=ctrl[2];ic=ctrl[3];il=ctrl[4];ig=ctrl[5];
  74     ia=ctrl[7];df=ctrl[10];dc=ctrl[11];db=ctrl[12];da=ctrl[13];dd=ctrl[16];
  75     ran=ctrl[18];can=ctrl[19];rap=ctrl[22];
  76     ea=ctrl[23];cae=ctrl[24];ral=ctrl[25];
  77 
  78     /* check for forced divergence (due to divergence of a user material
  79        routine */
  80 
  81     if(qa[2]>0.){idivergence=1;}else{idivergence=0;}
  82 
  83     if(*ithermal!=2){
  84         if(qa[0]>ea*qam[0]){
  85             if(*iit<=ip){c1[0]=ran;}
  86             else{c1[0]=rap;}
  87             c2[0]=can;
  88         }
  89         else{
  90             c1[0]=ea;
  91             c2[0]=cae;
  92         }
  93         if(ram1[0]<ram2[0]){ram2[0]=ram1[0];}
  94     }
  95     if(*ithermal>1){
  96         if(qa[1]>ea*qam[1]){
  97             if(*iit<=ip){c1[1]=ran;}
  98             else{c1[1]=rap;}
  99             c2[1]=can;
 100         }
 101         else{
 102             c1[1]=ea;
 103             c2[1]=cae;
 104         }
 105         if(ram1[1]<ram2[1]){ram2[1]=ram1[1];}
 106     }
 107 
 108     iconvergence=0; 
 109  
 110     /* mechanical */
 111 
 112     if(*ithermal<2){
 113         if((*iit>1)&&(ram[0]<=c1[0]*qam[0])&&(*iflagact==0)&&
 114 //      if((*iit>1)&&(ram[0]<=c1[0]*qam[0])&&
 115            ((cam[0]<=c2[0]*uam[0])||
 116             (((ram[0]*cam[0]<c2[0]*uam[0]*ram2[0])||(ram[0]<=ral*qam[0])||
 117               (qa[0]<=ea*qam[0]))&&(*ntg==0))||
 118             (cam[0]<1.e-8))) iconvergence=1;
 119     }
 120 
 121     /* thermal */
 122 
 123     if(*ithermal==2){
 124         if((ram[1]<=c1[1]*qam[1])&&
 125            (cam[2]<*deltmx)&&
 126            ((cam[1]<=c2[1]*uam[1])||
 127             (((ram[1]*cam[1]<c2[1]*uam[1]*ram2[1])||(ram[1]<=ral*qam[1])||
 128               (qa[1]<=ea*qam[1]))&&(*ntg==0))||
 129             (cam[1]<1.e-8)))iconvergence=1;
 130     }
 131 
 132     /* thermomechanical */
 133 
 134     if(*ithermal==3){
 135 //      if(((ram[0]<=c1[0]*qam[0])&&
 136         if(((*iit>1)&&(ram[0]<=c1[0]*qam[0])&&
 137             ((cam[0]<=c2[0]*uam[0])||
 138              (((ram[0]*cam[0]<c2[0]*uam[0]*ram2[0])||(ram[0]<=ral*qam[0])||
 139                (qa[0]<=ea*qam[0]))&&(*ntg==0))||
 140              (cam[0]<1.e-8)))&&
 141            ((ram[1]<=c1[1]*qam[1])&&
 142             (cam[2]<*deltmx)&&
 143             ((cam[1]<=c2[1]*uam[1])||
 144              (((ram[1]*cam[1]<c2[1]*uam[1]*ram2[1])||(ram[1]<=ral*qam[1])||
 145                (qa[1]<=ea*qam[1]))&&(*ntg==0))||
 146              (cam[1]<1.e-8))))iconvergence=1;
 147     }
 148 
 149     /* reset iflagact */
 150 
 151     *iflagact=0;
 152         
 153     /* increment convergence reached */
 154         
 155     if((iconvergence==1)&&(idivergence==0)){
 156 //      *ttime=*ttime+*dtime;
 157 
 158         /* cutting the insignificant digits from ttime */
 159 
 160 //      *ttime=*ttime+1.;
 161 //      *ttime=*ttime-1.;
 162         FORTRAN(writesummary,(istep,iinc,icutb,iit,ttime,time,dtime));
 163         if(*uncoupled){
 164             if(*ithermal==2){
 165                 *iitterm=*iit;
 166                 *ithermal=1;
 167                 for(k=0;k<*nk;++k){t1[k]=vold[mt*k];}
 168 //              *ttime=*ttime-*dtime;
 169                 *iit=1;
 170                 (ctrl[0])*=4;
 171                 printf(" thermal convergence\n\n");
 172                 return;
 173             }else{
 174                 *ithermal=3;
 175                 *iit=*iitterm;
 176                 (ctrl[0])/=4;
 177             }
 178         }
 179         
 180         *icntrl=1;
 181         *icutb=0;
 182         *theta=*theta+*dtheta;
 183         
 184         /* defining a mean "velocity" for static calculations: is used to
 185            extrapolate the present results for next increment */
 186         
 187         if(*nmethod != 4){
 188             for(i=0;i<*nk;i++){
 189                 for(j=1;j<mt;j++){
 190                     veold[mt*i+j]=(vold[mt*i+j]-vini[mt*i+j])/(*dtime);
 191                 }
 192             }
 193         }
 194         
 195         /* check whether next increment size must be decreased */
 196         
 197         if((*iit>il)&&(*idrct==0)){
 198             if(*mortar==0){
 199                 *dtheta=*dthetaref*db;
 200                 *dthetaref=*dtheta;
 201                 printf(" convergence; the increment size is decreased to %e\n\n",*dtheta**tper);
 202                 if(*dtheta<*tmin){
 203                     printf("\n *ERROR: increment size smaller than minimum\n");
 204                     printf(" best solution and residuals are in the frd file\n\n");
 205                     NNEW(fn,double,mt**nk);
 206                     NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 207                     FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
 208                       nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 209                       ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 210                     ++*kode;
 211 
 212                     (*ttime)+=(*time);
 213                     frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 214                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 215                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 216                         trab,inotr,ntrans,orab,ielorien,norien,description,
 217                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 218                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 219                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,
 220                         cdn,mortar,cdnr,cdni,nmat);
 221 
 222                     FORTRAN(stop,());
 223                 }
 224             }
 225             else{
 226                 printf("convergence\n\n");}
 227         }
 228         
 229         /* check whether next increment size can be increased */
 230         
 231         else if(*iit<=ig){
 232             if((*istab==1)&&(*idrct==0)){
 233                 *dtheta=*dthetaref*dd;
 234                 *dthetaref=*dtheta;
 235                 printf(" convergence; the increment size is increased to %e\n\n",*dtheta**tper);
 236             }
 237             else{
 238                 *istab=1;
 239                 printf(" convergence\n\n");
 240                 *dtheta=*dthetaref;
 241             }
 242         }
 243         else{
 244             *istab=0;
 245             printf(" convergence\n\n");
 246             *dtheta=*dthetaref;
 247         }
 248         
 249         /* check whether new increment size exceeds maximum increment
 250            size allowed (by the user) */
 251 
 252         if((*dtheta>*tmax)&&(*idrct==0)){
 253             *dtheta=*tmax;
 254             *dthetaref=*dtheta;
 255             printf(" the increment size exceeds thetamax and is decreased to %e\n\n",*dtheta**tper);
 256         }
 257 
 258         /* if itp=1 the increment just finished ends at a time point */
 259 
 260         if((*itpamp>0)&&(*idrct==0)){
 261             if(*itp==1){
 262                 *jprint=*jout;
 263             }else{
 264                 *jprint=*jout+1;
 265             }
 266             if(namta[3**itpamp-1]<0){
 267                 reftime=*ttime+*time+*dtheta**tper;
 268             }else{
 269                 reftime=*time+*dtheta**tper;
 270             }
 271             istart=namta[3**itpamp-3];
 272             iend=namta[3**itpamp-2];
 273             FORTRAN(identamta,(amta,&reftime,&istart,&iend,&id));
 274             if(id<istart){
 275                 inew=istart;
 276             }else{
 277                 inew=id+1;
 278             }
 279 
 280             /* inew: smallest time point exceeding time+dtheta*tper
 281                inext: smallest time point exceeding time */
 282 
 283             if(*inext<inew){
 284                 if(namta[3**itpamp-1]<0){
 285                     *dtheta=(amta[2**inext-2]-*ttime-*time)/(*tper);
 286                 }else{
 287                     *dtheta=(amta[2**inext-2]-*time)/(*tper);
 288                 }
 289                 (*inext)++;
 290                 *itp=1;
 291                 printf(" the increment size exceeds a time point and is decreased to %e\n\n",*dtheta**tper);
 292             }else{*itp=0;}
 293         }
 294 
 295         /* check whether new time point exceeds end of step */
 296 
 297         if(*dtheta>=1.-*theta){
 298             if(*dtheta>1.-*theta){iexceed=1;}else{iexceed=0;}
 299 //          if((*mortar==1)&&(1.-*theta>0.01)){
 300 //              *dtheta=0.999-*theta;
 301 //          if((*mortar==1)&&(1.-*theta>0.5)){
 302 //              *dtheta=0.900-*theta;
 303 //          }else{
 304                 *dtheta=1.-*theta;
 305 //          }
 306             *dthetaref=*dtheta;
 307             if(iexceed==1)
 308             printf(" the increment size exceeds the remainder of the step and is decreased to %e\n\n",*dtheta**tper);
 309 //          if(*dtheta<=1.e-6){
 310 //              (*ttime)+=(*dtheta**tper);
 311 //              *dtime=0.;
 312 //          }
 313         }
 314     }
 315     else{
 316 
 317         /* no convergence */
 318         
 319         /* check for the amount of iterations */
 320         
 321         if((*iit>ic)&&(*mortar==0)){
 322             printf("\n *ERROR: too many iterations needed\n");
 323             printf(" best solution and residuals are in the frd file\n\n");
 324 
 325             FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,dtime));
 326 
 327             NNEW(fn,double,mt**nk);
 328             NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 329             FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,nk,sti,stn,
 330                 ipkon,inum,kon,lakon,ne,mi,orab,ielorien,co,itg,ntg,vold,
 331                 ielmat,thicke,ielprop,prop));
 332             ++*kode;
 333 
 334             (*ttime)+=(*time);
 335             frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 336                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 337                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 338                         trab,inotr,ntrans,orab,ielorien,norien,description,
 339                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 340                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 341                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 342                         mortar,cdnr,cdni,nmat);
 343 
 344             FORTRAN(stop,());
 345         }       
 346         
 347         /* check for diverging residuals */
 348         
 349         if((*iit>=i0)||(fabs(ram[0])>1.e20)||(fabs(cam[0])>1.e20)||
 350                        (fabs(ram[1])>1.e20)||(fabs(cam[1])>1.e20)||
 351                        (cam[2]>*deltmx)||(qa[2]>0.)){
 352             if((*ithermal!=2)&&(*mortar==0)){
 353                 if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 354                     idivergence=1;
 355             }
 356 
 357 /*          if((*ithermal!=2)&&(*mortar==1)){
 358                 if((ram[0]>ram[6])&&(ram1[0]>ram[6])){
 359                     printf("divergence allowed: residual force not converging\n");
 360                     if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 361                         idivergence=1;
 362                 }
 363                 if(ram[5]<0.5){
 364                     printf("divergence allowed: number of contact elements stabilized\n");
 365                     if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 366                         idivergence=1;
 367                 }
 368                 if(ram[5]>=ram1[3]){
 369                     if(((ram[4]>0.97*ram1[2])&&(ram[4]<1.03*ram1[2]))&&
 370                        ((ram[4]>0.97*ram2[2])&&(ram[4]<1.03*ram2[2]))){
 371                         printf("divergence allowed: repetitive pattern detected\n");
 372                         if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 373                             idivergence=1;
 374                         if(((ram[4]>0.99*ram1[2])&&(ram[4]<1.01*ram1[2]))&&
 375                            ((ram[4]>0.99*ram2[2])&&(ram[4]<1.01*ram2[2])&&(ram[0]>c1[0]*qam[0]))){
 376                             printf("repetitive pattern within 1 per cent -> divergence\n");
 377                             idivergence=1;
 378                         }
 379                         if((ram[4]>(1-1.e-3)*ram1[2])&&(ram[4]<(1+1.e-3)*ram1[2])&&
 380                            (ram[4]>(1-1.e-3)*ram2[2])&&(ram[4]<(1+1.e-3)*ram2[2])&&
 381                            (ram[5]==ram1[3])&&(ram[5]==ram2[3])){
 382                             printf("infinite loop -> divergence\n");
 383                             idivergence=1;
 384                         }
 385                     }
 386                 }   
 387                 }*/
 388 
 389             if((*ithermal!=2)&&(*mortar==1)){
 390 
 391                 if(ram[0]>1.e9){
 392                     printf("divergence allowed: residual force too large\n");
 393                     if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 394                         idivergence=1;
 395                 }
 396 
 397                 /* number of contact elements does not change */
 398 
 399                 if(((ITG)ram[6]==0)&&((ITG)ram1[6]==0)){
 400                     printf("divergence allowed: number of contact elements stabilized\n");
 401                     if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 402                         idivergence=1;
 403                 }
 404             
 405                 /* rate of number of contact elements is increasing */
 406 
 407                 if(((ITG)ram[6]>=(ITG)ram1[6])&&((ITG)ram[6]>=(ITG)ram2[6])){
 408                     
 409                     if(((ram[4]>0.98*ram1[4])&&(ram[4]<1.02*ram1[4]))&&
 410                        ((ram[4]>0.98*ram2[4])&&(ram[4]<1.02*ram2[4]))){
 411                         printf("divergence allowed: repetitive pattern detected\n");
 412                         if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 413                             idivergence=1;
 414                     }
 415                     
 416                     /*      if(((ram[4]>0.99*ram1[4])&&(ram[4]<1.01*ram1[4]))&&
 417                        ((ram[4]>0.99*ram2[4])&&(ram[4]<1.01*ram2[4]))&&
 418                        (ram[0]>c1[0]*qam[0])){
 419                         printf("repetitive pattern within 1 %: divergence\n");
 420                         idivergence=1;
 421                     }
 422 
 423                     if(((ram[4]>(1.-1.e-3)*ram1[4])&&(ram[4]<(1.+1.e-3)*ram1[4]))&&
 424                        ((ram[4]>(1.-1.e-3)*ram2[4])&&(ram[4]<(1.+1.e-3)*ram2[4]))&&
 425                        ((ITG)ram[6]==(ITG)ram1[6])&&((ITG)ram[6]==(ITG)ram2[6])){
 426                         printf("infinite loop -> divergence\n");
 427                         idivergence=1;
 428                         }*/
 429                 }
 430                 
 431                 /*      if(((ITG)ram1[7]==1)&&((ITG)ram2[7]!=-1)){
 432                     printf("divergence if ram>ram1>ram2: contact elements not stabilizing\n");
 433                     if((ram1[0]>ram2[0])&&(ram[0]>ram2[0])&&(ram[0]>c1[0]*qam[0]))
 434                         idivergence=1;
 435                         }*/
 436             }
 437 
 438 
 439             /* for thermal calculations the maximum temperature change
 440                is checked as well */
 441 
 442             if(*ithermal>1){
 443                 if((ram1[1]>ram2[1])&&(ram[1]>ram2[1])&&(ram[1]>c1[1]*qam[1]))
 444                     idivergence=1;
 445                 if(cam[2]>*deltmx) idivergence=2;
 446             }
 447             if(idivergence>0){
 448                 if(*idrct==1) {
 449 
 450                     /* fixed time increments */
 451 
 452                     printf("\n *ERROR: solution seems to diverge; please try \n");
 453                     printf(" automatic incrementation; program stops\n");
 454                     printf(" best solution and residuals are in the frd file\n\n");
 455 
 456                     FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,
 457                                              dtime));
 458 
 459                     NNEW(fn,double,mt**nk);
 460                     NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 461                     FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,nk,
 462                        sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 463                        ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 464                     ++*kode;
 465 
 466                     (*ttime)+=(*time);
 467                     frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 468                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 469                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 470                         trab,inotr,ntrans,orab,ielorien,norien,description,
 471                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 472                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 473                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 474                         mortar,cdnr,cdni,nmat);
 475 
 476                     FORTRAN(stop,());
 477                 }
 478                 else {
 479 
 480                     /* variable time increments */
 481 
 482                     if(qa[2]>0.){
 483                         *dtheta=*dtheta*qa[2];
 484                         printf("increment size decrease requested by a material user routine (through pnewdt)\n\n");
 485                     }else{
 486                         if(idivergence==1){
 487                             *dtheta=*dtheta*df;
 488                         }else{
 489                             *dtheta=*dtheta**deltmx/cam[2]*da;
 490                         }
 491                     }
 492                     *dthetaref=*dtheta;
 493                     printf(" divergence; the increment size is decreased to %e\n",*dtheta**tper);
 494                     printf(" the increment is reattempted\n\n");
 495 
 496                     FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,time,
 497                                              dtime));
 498 
 499                     *istab=0;
 500                     if(*itp==1){
 501                       *itp=0;
 502                       (*inext)--;
 503                     }
 504 
 505                     /* check whether new increment size is smaller than minimum */
 506 
 507                     if(*dtheta<*tmin){
 508                         printf("\n *ERROR: increment size smaller than minimum\n");
 509                         printf(" best solution and residuals are in the frd file\n\n");
 510                         NNEW(fn,double,mt**nk);
 511                         NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 512                         FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
 513                            nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 514                            ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 515                         ++*kode;
 516 
 517                         (*ttime)+=(*time);
 518                         frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 519                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 520                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 521                         trab,inotr,ntrans,orab,ielorien,norien,description,
 522                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 523                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 524                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 525                         mortar,cdnr,cdni,nmat);
 526 
 527                         FORTRAN(stop,());
 528                     }
 529                     *icntrl=1;
 530                     (*icutb)++;
 531 
 532                     /* check whether too many cutbacks */
 533 
 534                     if(*icutb>ia){
 535                         printf("\n *ERROR: too many cutbacks\n");
 536                         printf(" best solution and residuals are in the frd file\n\n");
 537                         NNEW(fn,double,mt**nk);
 538                         NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 539                         FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
 540                            nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 541                            ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 542                         ++*kode;
 543 
 544                         (*ttime)+=(*time);
 545                         frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 546                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 547                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 548                         trab,inotr,ntrans,orab,ielorien,norien,description,
 549                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 550                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 551                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 552                         mortar,cdnr,cdni,nmat);
 553 
 554                         FORTRAN(stop,());
 555                     }
 556                     if(*uncoupled){
 557                       if(*ithermal==1){
 558                         (ctrl[0])/=4;
 559                       }
 560                       *ithermal=3;
 561                     }
 562 
 563                     /* default value for qa[2] */
 564 
 565                     qa[2]=-1.;
 566 
 567                     return;
 568                 }
 569             }
 570         }
 571         
 572         /* check for too slow convergence */
 573         
 574         if(*iit>=ir){
 575             if(*ithermal!=2){
 576                 iest1=(ITG)ceil(*iit+log(ran*qam[0]/(ram[0]))/
 577                                 log(ram[0]/(ram1[0])));
 578             }
 579             if(*ithermal>1){
 580                 iest2=(ITG)ceil(*iit+log(ran*qam[1]/(ram[1]))/
 581                                 log(ram[1]/(ram1[1])));
 582             }
 583             if(iest1>iest2){iest=iest1;}else{iest=iest2;}
 584             if((iest>0)&&(*mortar==0)){
 585             printf(" estimated number of iterations till convergence = %" ITGFORMAT "\n",
 586                    iest);
 587             }
 588             if((((iest>ic)||(*iit==ic))&&(*mortar==0))||((*mortar==1)&&(*iit==60))){
 589                 
 590                 if(*idrct!=1){
 591                     *dtheta=*dtheta*dc;
 592                     *dthetaref=*dtheta;
 593                     printf(" too slow convergence; the increment size is decreased to %e\n",*dtheta**tper);
 594                     printf(" the increment is reattempted\n\n");
 595 
 596                     FORTRAN(writesummarydiv,(istep,iinc,icutb,iit,ttime,
 597                                              time,dtime));
 598 
 599                     *istab=0;
 600 
 601                     /* check whether new increment size is smaller than minimum */
 602 
 603                     if(*dtheta<*tmin){
 604                         printf("\n *ERROR: increment size smaller than minimum\n");
 605                         printf(" best solution and residuals are in the frd file\n\n");
 606                         NNEW(fn,double,mt**nk);
 607                         NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 608                         FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
 609                            nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 610                            ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 611                         ++*kode;
 612 
 613                         (*ttime)+=(*time);
 614                         frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 615                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 616                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 617                         trab,inotr,ntrans,orab,ielorien,norien,description,
 618                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 619                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 620                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 621                         mortar,cdnr,cdni,nmat);
 622 
 623                         FORTRAN(stop,());
 624                     }
 625                     *icntrl=1;
 626                     (*icutb)++;
 627                     if(*icutb>ia){
 628                         printf("\n *ERROR: too many cutbacks\n");
 629                         printf(" best solution and residuals are in the frd file\n\n");
 630                         NNEW(fn,double,mt**nk);
 631                         NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 632                         FORTRAN(storeresidual,(nactdof,b,fn,filab,ithermal,
 633                            nk,sti,stn,ipkon,inum,kon,lakon,ne,mi,orab,
 634                            ielorien,co,itg,ntg,vold,ielmat,thicke,ielprop,prop));
 635                         ++*kode;
 636 
 637                         (*ttime)+=(*time);
 638                         frd(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,nmethod,
 639                         kode,filab,een,t1act,fn,ttime,epn,ielmat,matname,enern,
 640                         xstaten,nstate_,istep,iinc,ithermal,qfn,mode,noddiam,
 641                         trab,inotr,ntrans,orab,ielorien,norien,description,
 642                         ipneigh,neigh,mi,sti,vr,vi,stnr,stni,vmax,stnmax,
 643                         &ngraph,veold,ener,ne,cs,set,nset,istartset,iendset,
 644                         ialset,eenmax,fnr,fni,emn,thicke,jobnamec,output,qfx,cdn,
 645                         mortar,cdnr,cdni,nmat);
 646 
 647                         FORTRAN(stop,());
 648                     }
 649                     if(*uncoupled){
 650                       if(*ithermal==1){
 651                         (ctrl[0])/=4;
 652                       }
 653                       *ithermal=3;
 654                     }
 655                     return;
 656                 }
 657             }
 658         }
 659         
 660         printf(" no convergence\n\n");
 661         
 662         (*iit)++;
 663         
 664     }
 665     
 666     /* default value for qa[2] */
 667 
 668     /*  qa[2]=-1;*/
 669 
 670     return;
 671 }

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