root/src/dynacont.c

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

DEFINITIONS

This source file includes following definitions.
  1. dynacont

   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 /*     author: Reinhold Fischer                                          */
  19 
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include <time.h>
  24 #include <sys/time.h>
  25 #include <string.h>
  26 #include "CalculiX.h"
  27 
  28 #ifdef SPOOLES
  29    #include "spooles.h"
  30 #endif
  31 #ifdef SGI
  32    #include "sgi.h"
  33 #endif
  34 #ifdef TAUCS
  35    #include "tau.h"
  36 #endif
  37 
  38 void dynacont(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, 
  39               ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun,
  40               ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
  41               ITG *nmpc, ITG *nodeforc,ITG *ndirforc,double *xforc, 
  42               ITG *nforc,ITG *nelemload, char *sideload,double *xload,
  43               ITG *nload, 
  44               ITG *nactdof,ITG *neq, ITG *nzl,ITG *icol, ITG *irow, 
  45               ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, 
  46               ITG *ilboun,double *elcon, ITG *nelcon, double *rhcon, 
  47               ITG *nrhcon,double *cocon, ITG *ncocon,
  48               double *alcon, ITG *nalcon, double *alzero, 
  49               ITG *ielmat,ITG *ielorien, ITG *norien, double *orab, 
  50               ITG *ntmat_,double *t0, 
  51               double *t1,ITG *ithermal,double *prestr, ITG *iprestr, 
  52               double *vold,ITG *iperturb, double *sti, ITG *nzs, 
  53               double *tinc, double *tper, double *xmodalsteady,
  54               double *veold, char *amname, double *amta,
  55               ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload,
  56               ITG *iamt1,ITG *jout,char *filab,double *eme, double *xforcold, 
  57               double *xloadold,
  58               double *t1old, ITG *iamboun, double *xbounold, ITG *iexpl,
  59               double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
  60               double *xstate, ITG *npmat_, char *matname, ITG *mi,
  61               ITG *ncmat_, ITG *nstate_, double *ener, char *jobnamec,
  62               double *ttime, char *set, ITG *nset, ITG *istartset,
  63               ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  64               char *prset, ITG *nener, double *trab, 
  65               ITG *inotr, ITG *ntrans, double *fmpc, char *cbody, ITG *ibody,
  66               double *xbody, ITG *nbody, double *xbodyold, ITG *istep,
  67               ITG *isolver,ITG *jq, char *output, ITG *mcs, ITG *nkon,
  68               ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset,
  69               ITG *idrct, ITG *jmax, double *tmin, double *tmax,
  70               double *ctrl, ITG *itpamp, double *tietol,ITG *iit,
  71               ITG *ncont,ITG *ne0, double *reltime, double *dtime,
  72               double *bcontini, double *bj, double *aux, ITG *iaux,
  73               double *bcont, ITG *nev, double *v,
  74               ITG *nkon0, double *deltmx, double *dtheta, double *theta,
  75               ITG *iprescribedboundary, ITG *mpcfree, ITG *memmpc_,
  76               ITG *itietri, ITG *koncont, double *cg, double *straight,
  77               ITG *iinc, double *vini,
  78               double *aa, double *bb, double *aanew, double *d, 
  79               double *z, double *zeta,double *b, double *time0,double *time, 
  80               ITG *ipobody,
  81               double *xforcact, double *xloadact, double *t1act, 
  82               double *xbounact, double *xbodyact, double *cd, double *cv,
  83               double *ampli, double *dthetaref, double *bjp, double *bp,
  84               double *cstr,ITG *imddof, ITG *nmddof, 
  85               ITG **ikactcontp, ITG *nactcont,ITG *nactcont_,
  86               double *aamech, double *bprev, ITG *iprev, ITG *inonlinmpc,
  87               ITG **ikactmechp, ITG *nactmech,ITG *imdnode,ITG *nmdnode,
  88               ITG *imdboun,ITG *nmdboun,ITG *imdmpc,ITG *nmdmpc,
  89               ITG *itp, ITG *inext,ITG *imastop,ITG *nslavnode,ITG *islavnode,
  90               ITG *islavsurf,
  91               ITG *itiefac,double *areaslav,ITG *iponoels,ITG *inoels,
  92               double *springarea,ITG *izdof,ITG *nzdof,double *fn,
  93               ITG *imastnode,ITG *nmastnode,double *xmastnor,
  94               double *xstateini, ITG *nslavs,
  95               ITG *cyclicsymmetry,double *xnoels,ITG *ielas,ITG *ielprop,
  96               double *prop){
  97     
  98   char lakonl[9]="        \0";
  99 
 100   ITG i,j,k,l,init,*itg=NULL,ntg=0,maxlenmpc,icascade=0,loop,
 101       konl[20],imat,nope,kodem,indexe,j1,jdof,kmin,kmax,
 102       id,newstep=0,idiscon,*ipiv=NULL,info,nrhs=1,kode,iener=0,
 103       *ikactcont=NULL,*ilactcont=NULL,*ikactcont1=NULL,nactcont1=0,
 104       i1,icutb=0,iconvergence=0,idivergence=0,mt=mi[1]+1,
 105       nactcont1_=100,*ikactmech=NULL,iabsload=0,im,nasym=0,mortar=0;
 106 
 107   long long i2;
 108 
 109   double *adb=NULL,*aub=NULL,*cgr=NULL, *au=NULL,fexp,fcos,fsin,fexm,
 110       physcon[1],zetaj,dj,ddj,h1,h2,h3,h4,h5,h6,sum,aai,bbi,tstart,tend,
 111       *ad=NULL,sigma=0.,alpham,betam,*bact=NULL,*bmin=NULL,*bv=NULL,
 112       xl[27],voldl[mt*9],elas[21],fnl[27],t1l,elconloc[21],veoldl[mt*9],
 113       bbmax,s[3600],*aaa=NULL,*bbb=NULL,func,funcp,*bjbasp=NULL,
 114       *bjbas=NULL, *bjinc=NULL, *dbj=NULL, *lhs=NULL,dbjmax,bjmax,
 115       *bjincp=NULL,sump,h14,*dbjp=NULL,senergy=0.0,*xforcdiff=NULL,
 116       df,i0,ic,ia,dbjmaxOLD1,dbjmaxOLD2,*xloaddiff=NULL,*dbcont=NULL,
 117       zl=0.0,*xbodydiff=NULL,*t1diff=NULL,*xboundiff=NULL,*bdiff=NULL,
 118       *pslavsurf=NULL,*pmastsurf=NULL,*clearini=NULL;
 119 
 120   ikactcont=*ikactcontp;ikactmech=*ikactmechp;
 121 
 122   if(*inonlinmpc==1) iabsload=2;
 123 
 124   if(ithermal[0]<=1){
 125       kmin=1;kmax=3;
 126   }else if(ithermal[0]==2){
 127       kmin=0;kmax=mi[1];if(kmax>2)kmax=2;
 128   }else{
 129       kmin=0;kmax=3;
 130   }
 131 
 132   NNEW(xforcdiff,double,*nforc);
 133   NNEW(xloaddiff,double,2**nload);
 134   NNEW(xbodydiff,double,7**nbody);
 135 
 136   /* copying the rotation axis and/or acceleration vector */
 137 
 138   for(k=0;k<7**nbody;k++){xbodydiff[k]=xbody[k];}
 139   NNEW(xboundiff,double,*nboun);
 140   if(*ithermal==1) NNEW(t1diff,double,*nk);
 141 
 142   /* load the convergence constants from ctrl*/
 143 
 144   i0=ctrl[0];ic=ctrl[3];ia=ctrl[7];df=ctrl[10];
 145 
 146   /* set the convergence parameters*/
 147 
 148   dbjmaxOLD1=0.0;
 149   dbjmaxOLD2=0.0;
 150 
 151   /* calculating the contact forces */
 152         
 153   for(j=0;j<*nactcont;j++){bcont[ikactcont[j]]=0.;}
 154   
 155   *ne=*ne0;*nkon=*nkon0;
 156   
 157   contact(ncont,ntie,tieset,nset,set,istartset,iendset,
 158           ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,
 159           straight,nkon,co,vold,ielmat,cs,elcon,istep,
 160           iinc,iit,ncmat_,ntmat_,ne0,
 161           vini,nmethod,nmpc,mpcfree,memmpc_,
 162           &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 163           iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
 164           itiefac,areaslav,iponoels,inoels,springarea,tietol,reltime,
 165           imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
 166           xnoels,&mortar,pslavsurf,pmastsurf,clearini,theta);
 167 
 168   NNEW(ikactcont1,ITG,nactcont1_);
 169 
 170   for(i=*ne0;i<*ne;i++){
 171       indexe=ipkon[i];
 172       imat=ielmat[mi[2]*i];
 173       kodem=nelcon[2*imat-2];
 174       for(j=0;j<8;j++){lakonl[j]=lakon[8*i+j];}
 175       nope=atoi(&lakonl[7])+1;
 176       for(j=0;j<nope;j++){
 177           konl[j]=kon[indexe+j];
 178           for(j1=0;j1<3;j1++){
 179               xl[j*3+j1]=co[3*(konl[j]-1)+j1];
 180               voldl[mt*j+j1+1]=vold[mt*(konl[j]-1)+j1+1];
 181               veoldl[mt*j+j1+1]=veold[mt*(konl[j]-1)+j1+1];
 182           }
 183       }
 184       konl[nope]=kon[indexe+nope];
 185 
 186       FORTRAN(springforc_n2f,(xl,konl,voldl,&imat,elcon,nelcon,elas,
 187                           fnl,ncmat_,ntmat_,&nope,lakonl,
 188                           &t1l,&kodem,elconloc,plicon,nplicon,npmat_,
 189                           &senergy,&iener,cstr,mi,
 190                           &springarea[2*(konl[nope]-1)],nmethod,ne0,
 191                           nstate_,xstateini,xstate,reltime,ielas));
 192 
 193       storecontactdof(&nope,nactdof,&mt,konl,&ikactcont1,&nactcont1,
 194                       &nactcont1_,bcont,fnl,ikmpc,nmpc,ilmpc,ipompc,nodempc, 
 195                       coefmpc);
 196 
 197   }
 198   RENEW(ikactcont1,ITG,nactcont1);
 199 
 200   /* merging ikactcont with ikactcont1; the result ist
 201      stored in ikactcont */
 202 
 203   for(i=0;i<nactcont1;i++){
 204       jdof=ikactcont1[i];
 205       FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
 206       do{
 207           if(id>0){
 208               if(ikactcont[id-1]==jdof){
 209                   break;
 210               }
 211           }
 212           (*nactcont)++;
 213           if(*nactcont>*nactcont_){
 214               *nactcont_=(ITG)(1.1**nactcont_);
 215               RENEW(ikactcont,ITG,*nactcont_);
 216           }
 217           k=*nactcont-1;
 218           l=k-1;
 219           while(k>id){
 220               ikactcont[k--]=ikactcont[l--];
 221           }
 222           ikactcont[id]=jdof;
 223           break;
 224       }while(1);
 225   }
 226     
 227     /* calculate the change in contact force */
 228     
 229   bbmax=0.;
 230   if(icutb==0){
 231       for(i=0;i<*nactcont;i++){
 232           jdof=ikactcont[i];
 233           if(fabs(bcont[jdof]-bcontini[jdof])>bbmax){
 234               bbmax=fabs(bcont[jdof]-bcontini[jdof]);
 235           }
 236       }
 237   }
 238 
 239   /* removing entries in bcont */
 240 
 241   for(j=0;j<nactcont1;j++){bcont[ikactcont1[j]]=0.;}
 242   SFREE(ikactcont1);
 243   *nactcont=0;
 244   
 245   /* major loop to calculate the correction of bj due to contact */
 246 
 247   NNEW(ilactcont,ITG,*nactcont_);
 248   NNEW(dbcont,double,*nactcont_**nev);
 249 
 250   icutb=0;
 251 
 252   do{
 253 
 254       /* restoring initial values */
 255 
 256       if(*nmdnode>0){
 257           for(i=0;i<*nmdnode;i++){
 258               i1=mt*(imdnode[i]-1);
 259               for(j=kmin;j<=kmax;j++){
 260                   vold[i1+j]=vini[i1+j];
 261               }
 262           }
 263       }else{
 264           memcpy(&vold[0],&vini[0],sizeof(double)*mt**nk);
 265       }
 266 
 267       if(*nstate_!=0){
 268         for(k=0;k<*nstate_*mi[0]*(*ne0+*nslavs);++k){
 269           xstate[k]=xstateini[k];
 270         }       
 271       }
 272     
 273     /* restoring aa[(iinc-1)*nev+i] (before change of *dtime) */
 274 
 275     for(i=0;i<*nev;i++){
 276       aa[i]+=bb[i]*(*time-*dtime);
 277     }
 278     
 279     /* increment size is reduced if:
 280             - the contact force change is too large (only in first iteration)
 281             - or the increment did not converge  */
 282     
 283     if((bbmax>*deltmx || icutb>0)&&(((*itp==1)&&(*dtheta>*tmin))||(*itp==0))){
 284       
 285       /* force increase too big: increment size is decreased */
 286 
 287       if(icutb>0){
 288         *dtheta=*dtheta*df;
 289       }
 290       else{
 291         *dtheta=*dtheta**deltmx/bbmax;
 292       }
 293       *dthetaref=*dtheta;
 294       if(*itp==1){
 295           (*inext)--;
 296           *itp=0;
 297       }
 298       
 299       /* check whether the new increment size is not too small */
 300       
 301       if(*dtheta<*tmin){
 302           *dtheta=*tmin;
 303           *dthetaref=*dtheta;
 304       }
 305       
 306       *reltime=*theta+(*dtheta);
 307       *time=*reltime**tper;
 308       *dtime=*dtheta**tper;
 309       
 310       /* calculating the instantaneous loads (forces, surface loading, 
 311          centrifugal and gravity loading or temperature) */
 312       
 313       FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
 314                 xloadold,xload,xloadact,iamload,nload,ibody,xbody,
 315                 nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
 316                 namta,nam,ampli,time,reltime,ttime,dtime,ithermal,
 317                 nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
 318                 ndirboun,nodeforc,
 319                 ndirforc,istep,iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
 320                 nelemload,sideload,mi,
 321                 xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
 322                 iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
 323                 fn));
 324               
 325       /* calculating the instantaneous loading vector */
 326       
 327       if(iabsload!=2){
 328           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
 329                    ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcdiff,
 330                    nforc,nelemload,sideload,xloaddiff,nload,xbodydiff,
 331                    ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
 332                    ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
 333                    alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
 334                    t0,t1diff,ithermal,iprestr,vold,iperturb,iexpl,plicon,
 335                    nplicon,plkcon,nplkcon,
 336                    npmat_,ttime,time,istep,iinc,dtime,physcon,ibody,
 337                    xbodyold,reltime,veold,matname,mi,ikactmech,nactmech,
 338                    ielprop,prop));
 339       }else{
 340           FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
 341                    ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
 342                    nforc,nelemload,sideload,xloadact,nload,xbodyact,
 343                    ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
 344                    ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
 345                    alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
 346                    t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
 347                    nplicon,plkcon,nplkcon,
 348                    npmat_,ttime,time,istep,iinc,dtime,physcon,ibody,
 349                    xbodyold,reltime,veold,matname,mi,ikactmech,nactmech,
 350                    ielprop,prop));
 351       }
 352       
 353       /* correction for nonzero SPC's */
 354       
 355       if(*iprescribedboundary){
 356           dynboun(amta,namta,nam,ampli,time,ttime,dtime,
 357                   xbounold,xboun,
 358                   xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
 359                   aub,icol,irow,neq,nzs,&sigma,b,isolver,
 360                   &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
 361                   bprev,bdiff,nactmech,&iabsload,iprev);
 362       }
 363 
 364       /* correcting aamech */
 365 
 366       if(!(*cyclicsymmetry)){
 367           for(i=0;i<*nev;i++){
 368               i2=(long long)i*neq[1];
 369               
 370               if(iabsload==2){aamech[i]=0.;}
 371               if(*nactmech<neq[1]/2){
 372                   for(j=0;j<*nactmech;j++){
 373                       aamech[i]+=z[i2+ikactmech[j]]*b[ikactmech[j]];
 374                   }
 375               }else{
 376                   for(j=0;j<neq[1];j++){
 377                       aamech[i]+=z[i2+j]*b[j];
 378                   }
 379               }
 380           }
 381       }else{
 382           for(i=0;i<*nev;i++){
 383               if(iabsload==2){aamech[i]=0.;}
 384           }
 385           for(j=0;j<*nactmech;j++){
 386               FORTRAN(nident,(izdof,&ikactmech[j],nzdof,&id));
 387               if(id!=0){
 388                   if(izdof[id-1]==ikactmech[j]){
 389                       for(i=0;i<*nev;i++){
 390                           aamech[i]+=z[i**nzdof+id-1]*b[ikactmech[j]];
 391                       }
 392                   }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 393               }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 394           }
 395       }
 396       
 397     }
 398 
 399     bbmax=0.;
 400     
 401     /* calculating the linearized force function connecting the
 402        mechanical+contact load at the start of the increment with
 403        the mechanical load at the end of the increment
 404        = base load */
 405 
 406     for(i=0;i<*nev;i++){
 407 
 408       aanew[i]=aamech[i];
 409 
 410       bb[i]=(aanew[i]-aa[i])/(*dtime);
 411       aa[i]=aanew[i]-bb[i]**time;
 412     }
 413     
 414     /* calculating the base response */
 415 
 416     NNEW(bjbas,double,*nev); /* basis response modal decomposition */
 417     NNEW(bjbasp,double,*nev);
 418     for(l=0;l<*nev;l++){
 419       zetaj=zeta[l];
 420       dj=d[l];
 421       
 422       /* zero eigenfrequency: rigid body mode */
 423       
 424       if(fabs(d[l])<=1.e-10){
 425           aai=aa[l];
 426           bbi=bb[l];
 427           tstart=*time0;
 428           tend=*time;
 429           sum=tend*(aai**time+
 430                     tend*((bbi**time-aai)/2.-bbi*tend/3.))-
 431               tstart*(aai**time+
 432                       tstart*((bbi**time-aai)/2.-bbi*tstart/3.));
 433           sump=tend*(aai+bbi*tend/2.)-tstart*(aai+bbi*tstart/2.);
 434           bjbas[l]=sum+cd[l]+*dtime*cv[l];
 435           bjbasp[l]=sump+cv[l];
 436       }
 437       
 438       /*   subcritical damping */
 439       
 440       else if(zetaj<1.-1.e-6){
 441         ddj=dj*sqrt(1.-zetaj*zetaj);
 442         h1=zetaj*dj;
 443         h2=h1*h1+ddj*ddj;
 444         h3=h1*h1-ddj*ddj;
 445         h4=2.*h1*ddj/h2;
 446         h14=zetaj*dj/ddj;
 447         tstart=0;
 448         FORTRAN(fsub,(time,dtime,&aa[l],&bb[l],&ddj,
 449                       &h1,&h2,&h3,&h4,&func,&funcp));
 450         sum=func;sump=funcp;
 451         FORTRAN(fsub,(time,&tstart,&aa[l],&bb[l],&ddj,
 452                       &h1,&h2,&h3,&h4,&func,&funcp));
 453         sum-=func;sump-=funcp;
 454         fexp=exp(-h1**dtime);
 455         fsin=sin(ddj**dtime);
 456         fcos=cos(ddj**dtime);
 457         
 458         bjbas[l]=sum/ddj+fexp*(fcos+zetaj/sqrt(1.-zetaj*zetaj)*fsin)*cd[l]+
 459           fexp*fsin*cv[l]/ddj;
 460         bjbasp[l]=sump/ddj+fexp*((-h1+ddj*h14)*fcos+(-ddj-h1*h14)*fsin)*cd[l]
 461           +fexp*(-h1*fsin+ddj*fcos)*cv[l]/ddj;
 462         
 463       }
 464       
 465       /*      supercritical damping */
 466       
 467       else if(zetaj>1.+1.e-6){
 468         ddj=dj*sqrt(zetaj*zetaj-1.);
 469         h1=ddj-zetaj*dj;
 470         h2=ddj+zetaj*dj;
 471         h3=1./h1;
 472         h4=1./h2;
 473         h5=h3*h3;
 474         h6=h4*h4;
 475         tstart=0;
 476         FORTRAN(fsuper,(time,dtime,&aa[l],&bb[l],
 477                         &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
 478         sum=func;sump=funcp;
 479         FORTRAN(fsuper,(time,&tstart,&aa[l],&bb[l],
 480                         &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
 481         sum-=func;sump-=funcp;
 482         fexm=exp(h1**dtime);
 483         fexp=exp(-h2**dtime);
 484         h14=zetaj*dj/ddj;
 485         
 486         bjbas[l]=sum/(2.*ddj)+(fexm+fexp)*cd[l]/2.+zetaj*(fexm-fexp)/(2.*sqrt(zetaj*zetaj-1.))*cd[l]+(fexm-fexp)*cv[l]/(2.*ddj);
 487         bjbasp[l]=sump/(2.*ddj)+(h1*fexm-h2*fexp)*cd[l]/2.+(h14*cd[l]+cv[l]/ddj)*(h1*fexm+h2*fexp)/2.;
 488       }
 489       
 490       /* critical damping */
 491       
 492       else{
 493         h1=zetaj*dj;
 494         h2=1./h1;
 495         h3=h2*h2;
 496         h4=h2*h3;
 497         tstart=0;
 498         FORTRAN(fcrit,(time,dtime,&aa[l],&bb[l],&zetaj,&dj,
 499                        &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
 500         sum=func;sump=funcp;
 501         FORTRAN(fcrit,(time,&tstart,&aa[l],&bb[l],&zetaj,&dj,
 502                        &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
 503         sum-=func;sump+=funcp;
 504         fexp=exp(-h1**dtime);
 505         bjbas[l]=sum+fexp*((1.+h1**dtime)*cd[l]+*dtime*cv[l]);
 506         bjbasp[l]=sump+fexp*(-h1*h1**dtime*cd[l]+(1.-h1**dtime)*cv[l]);
 507       }
 508     }
 509     
 510     /* calculating the incremental response due to contact */
 511     
 512     aai=-(*time-*dtime)/(*dtime);
 513     bbi=1./(*dtime);
 514     
 515     NNEW(bjinc,double,*nev); /* incremental response modal decomposition */
 516     NNEW(bjincp,double,*nev);
 517     for(l=0;l<*nev;l++){
 518       zetaj=zeta[l];
 519       dj=d[l];
 520       
 521       /* zero eigenfrequency: rigid body mode */
 522       
 523       if(fabs(d[l])<=1.e-10){
 524         tstart=*time0;
 525         tend=*time;
 526         sum=tend*(aai**time+
 527                   tend*((bbi**time-aai)/2.-bbi*tend/3.))-
 528           tstart*(aai**time+
 529                   tstart*((bbi**time-aai)/2.-bbi*tstart/3.));
 530         sump=tend*(aai+bbi*tend/2.)-tstart*(aai+bbi*tstart/2.);
 531         
 532         bjinc[l]=sum;
 533         bjincp[l]=sump;
 534       }
 535       
 536       /*   subcritical damping */
 537       
 538       else if(zetaj<1.-1.e-6){
 539         ddj=dj*sqrt(1.-zetaj*zetaj);
 540         h1=zetaj*dj;
 541         h2=h1*h1+ddj*ddj;
 542         h3=h1*h1-ddj*ddj;
 543         h4=2.*h1*ddj/h2;
 544         tstart=0.;
 545         FORTRAN(fsub,(time,dtime,&aai,&bbi,&ddj,
 546                       &h1,&h2,&h3,&h4,&func,&funcp));
 547         sum=func;sump=funcp;
 548         FORTRAN(fsub,(time,&tstart,&aai,&bbi,&ddj,
 549                       &h1,&h2,&h3,&h4,&func,&funcp));
 550         sum-=func;sump-=funcp;
 551         
 552         bjinc[l]=sum/ddj;
 553         bjincp[l]=sump/ddj;
 554         
 555       }
 556       
 557       /*      supercritical damping */
 558       
 559       else if(zetaj>1.+1.e-6){
 560         ddj=dj*sqrt(zetaj*zetaj-1.);
 561         h1=ddj-zetaj*dj;
 562         h2=ddj+zetaj*dj;
 563         h3=1./h1;
 564         h4=1./h2;
 565         h5=h3*h3;
 566         h6=h4*h4;
 567         tstart=0.;
 568         FORTRAN(fsuper,(time,dtime,&aai,&bbi,
 569                         &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
 570         sum=func;sump=funcp;
 571         FORTRAN(fsuper,(time,&tstart,&aai,&bbi,
 572                         &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
 573         sum-=func;sump-=funcp;
 574         
 575         bjinc[l]=sum/(2.*ddj);
 576         bjincp[l]=sump/(2.*ddj);
 577         
 578       }
 579       
 580       /* critical damping */
 581       
 582       else{
 583         h1=zetaj*dj;
 584         h2=1./h1;
 585         h3=h2*h2;
 586         h4=h2*h3;
 587         tstart=0.;
 588         FORTRAN(fcrit,(time,dtime,&aai,&bbi,&zetaj,&dj,
 589                        &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
 590         sum=func;sump=funcp;
 591         FORTRAN(fcrit,(time,&tstart,&aai,&bbi,&zetaj,&dj,
 592                        &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
 593         sum-=func;sump-=funcp;
 594         
 595         bjinc[l]=sum;
 596         bjincp[l]=sump;
 597         
 598       }
 599     }
 600 
 601     NNEW(aaa,double,*nev);
 602     NNEW(bbb,double,*nev**nev);
 603     NNEW(lhs,double,*nev**nev);
 604     NNEW(ipiv,ITG,*nev);
 605     NNEW(dbj,double,*nev); /* change in bj */
 606     NNEW(dbjp,double,*nev); /* change in djp */
 607     
 608     /* starting solution for the iteration loop = base solution */
 609 
 610     memcpy(&bj[0],&bjbas[0],sizeof(double)**nev);
 611     memcpy(&bjp[0],&bjbasp[0],sizeof(double)**nev);
 612     
 613     /* major iteration loop for the contact response */
 614     
 615     loop=0;
 616     do{
 617       loop++;
 618       
 619       /* composing the response */
 620       
 621       if(*iprescribedboundary){
 622           if(*nmdnode==0){
 623               memcpy(&b[0],&bmin[0],sizeof(double)*neq[1]);
 624               memcpy(&bp[0],&bv[0],sizeof(double)*neq[1]);
 625           }else{
 626               for(i=0;i<*nmddof;i++){
 627                   b[imddof[i]]=bmin[imddof[i]];
 628                   bp[imddof[i]]=bv[imddof[i]];
 629               }
 630           }
 631       }
 632       else{
 633           if(*nmdnode==0){
 634               DMEMSET(b,0,neq[1],0.);
 635               DMEMSET(bp,0,neq[1],0.);
 636           }else{
 637               for(i=0;i<*nmddof;i++){
 638                   b[imddof[i]]=0.;
 639                   bp[imddof[i]]=0.;
 640               }
 641           }
 642       }
 643       
 644       if(!(*cyclicsymmetry)){
 645           if(*nmdnode==0){
 646               for(i=0;i<neq[1];i++){
 647                   for(j=0;j<*nev;j++){
 648                       b[i]+=bj[j]*z[(long long)j*neq[1]+i];
 649                       bp[i]+=bjp[j]*z[(long long)j*neq[1]+i];
 650                   }
 651               }
 652           }else{
 653               for(i=0;i<*nmddof;i++){
 654                   for(j=0;j<*nev;j++){
 655                       b[imddof[i]]+=bj[j]*z[(long long)j*neq[1]+imddof[i]];
 656                       bp[imddof[i]]+=bjp[j]*z[(long long)j*neq[1]+imddof[i]];
 657                   }
 658               }
 659           }
 660       }else{
 661           for(i=0;i<*nmddof;i++){
 662               FORTRAN(nident,(izdof,&imddof[i],nzdof,&id));
 663               if(id!=0){
 664                   if(izdof[id-1]==imddof[i]){
 665                       for(j=0;j<*nev;j++){
 666                           b[imddof[i]]+=bj[j]*z[j**nzdof+id-1];
 667                           bp[imddof[i]]+=bjp[j]*z[j**nzdof+id-1];
 668                       }
 669                   }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 670               }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 671           }
 672       }
 673       
 674       /* update nonlinear MPC-coefficients (e.g. for rigid
 675          body MPC's */
 676 
 677       if(*inonlinmpc==1){
 678           FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
 679                          nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
 680                          &maxlenmpc,ikmpc,ilmpc,&icascade,
 681                          kon,ipkon,lakon,ne,reltime,&newstep,xboun,fmpc,
 682                          iit,&idiscon,ncont,trab,ntrans,ithermal,mi));
 683       }
 684       
 685       /* calculating displacements/temperatures */
 686       
 687       FORTRAN(dynresults,(nk,v,ithermal,nactdof,vold,nodeboun,
 688               ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
 689               b,bp,veold,dtime,mi,imdnode,nmdnode,imdboun,nmdboun,
 690               imdmpc,nmdmpc,nmethod,time));
 691       
 692       if((iconvergence==1)||((*idrct==1)&&(loop>1))){
 693         break;
 694       }
 695       
 696       /* creating contact elements and calculating the contact forces
 697          based on the displacements at the end of the present increment */
 698       
 699           for(j=0;j<*nactcont;j++){bcont[ikactcont[j]]=0.;}
 700 
 701       RENEW(dbcont,double,*nactcont_**nev);
 702       RENEW(ikactcont,ITG,*nactcont_);
 703       RENEW(ilactcont,ITG,*nactcont_);
 704       *nactcont=0;
 705       
 706       DMEMSET(dbcont,0,*nactcont_**nev,0.);
 707       DMEMSET(ikactcont,0,*nactcont_,0.);
 708       
 709       *ne=*ne0;*nkon=*nkon0;
 710       contact(ncont,ntie,tieset,nset,set,istartset,iendset,
 711               ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,
 712               straight,nkon,co,vold,ielmat,cs,elcon,istep,
 713               iinc,iit,ncmat_,ntmat_,ne0,
 714               vini,nmethod,nmpc,mpcfree,memmpc_,
 715               &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 716               iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
 717               itiefac,areaslav,iponoels,inoels,springarea,tietol,reltime,
 718               imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
 719               xnoels,&mortar,pslavsurf,pmastsurf,clearini,theta);
 720 
 721       for(i=*ne0;i<*ne;i++){
 722         indexe=ipkon[i];
 723         imat=ielmat[mi[2]*i];
 724         kodem=nelcon[2*imat-2];
 725         for(j=0;j<8;j++){lakonl[j]=lakon[8*i+j];}
 726         nope=atoi(&lakonl[7])+1;
 727         for(j=0;j<nope;j++){
 728           konl[j]=kon[indexe+j];
 729           for(j1=0;j1<3;j1++){
 730             xl[j*3+j1]=co[3*(konl[j]-1)+j1];
 731             voldl[mt*j+j1+1]=vold[mt*(konl[j]-1)+j1+1];
 732             veoldl[mt*j+j1+1]=veold[mt*(konl[j]-1)+j1+1];
 733           }
 734         }
 735         konl[nope]=kon[indexe+nope];
 736 
 737         FORTRAN(springforc_n2f,(xl,konl,voldl,&imat,elcon,nelcon,elas,
 738                             fnl,ncmat_,ntmat_,&nope,lakonl,
 739                             &t1l,&kodem,elconloc,plicon,nplicon,npmat_,
 740                             &senergy,&iener,cstr,mi,
 741                             &springarea[2*(konl[nope]-1)],nmethod,ne0,
 742                             nstate_,xstateini,xstate,reltime,ielas));
 743         
 744         FORTRAN(springstiff_n2f,(xl,elas,konl,voldl,s,&imat,elcon,nelcon,
 745                 ncmat_,ntmat_,&nope,lakonl,&t1l,&kode,elconloc,
 746                 plicon,nplicon,npmat_,iperturb,&springarea[2*(konl[nope]-1)],
 747                 nmethod,mi,ne0,nstate_,xstateini,xstate,reltime,&nasym));
 748 
 749         dfdbj(bcont,&dbcont,&neq[1],&nope,konl,nactdof,
 750               s,z,ikmpc,ilmpc,ipompc,nodempc,nmpc,coefmpc,
 751               fnl,nev,&ikactcont,&ilactcont,nactcont,nactcont_,mi,
 752               cyclicsymmetry,izdof,nzdof);
 753       }
 754 
 755       if(*nactcont>100){*nactcont_=*nactcont;}else{*nactcont_=100;}
 756       RENEW(ikactcont,ITG,*nactcont_);
 757       RENEW(ilactcont,ITG,*nactcont_);
 758       RENEW(dbcont,double,*nactcont_**nev);
 759 
 760       /* aaa(i) is the internal product of the contact force at the end of the
 761          increment with eigenmode i
 762          bbb(i,j) is the internal product of the change of the contact force with
 763          respect to modal coordinate j with the eigenmode i */
 764 
 765       DMEMSET(bbb,0,*nev**nev,0.);
 766       DMEMSET(aaa,0,*nev,0.);
 767 
 768       if(!(*cyclicsymmetry)){
 769           for(k=0; k<*nactcont; k++){
 770               i1=ikactcont[k];
 771               i2=(ilactcont[k]-1)**nev;
 772               for(j=0; j<*nev; j++){
 773                   zl=z[(long long)j*neq[1]+i1];
 774                   aaa[j]+=zl*bcont[i1];
 775                   for(l=0; l<*nev; l++){
 776                       bbb[l**nev+j]+=zl*dbcont[i2+l];
 777                   }
 778               }
 779           }
 780       }else{
 781           for(k=0; k<*nactcont; k++){
 782               i1=ikactcont[k];
 783               i2=(ilactcont[k]-1)**nev;
 784               FORTRAN(nident,(izdof,&i1,nzdof,&id));
 785               if(id!=0){
 786                   if(izdof[id-1]==i1){
 787                       for(j=0; j<*nev; j++){
 788                           zl=z[j**nzdof+id-1];
 789                           aaa[j]+=zl*bcont[i1];
 790                           for(l=0; l<*nev; l++){
 791                               bbb[l**nev+j]+=zl*dbcont[i2+l];
 792                           }
 793                       }
 794                   }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 795               }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 796           }
 797       }
 798       
 799       for(l=0;l<*nev;l++){
 800         i1=l**nev;
 801         for(j=0;j<*nev;j++){
 802           if(j==l){lhs[i1+j]=1.;}else{lhs[i1+j]=0.;}
 803           lhs[i1+j]-=bjinc[j]*bbb[i1+j];
 804         }
 805         dbj[l]=bjbas[l]+bjinc[l]*aaa[l]-bj[l];
 806       }
 807 
 808       /* solve the system of equations; determine dbj */
 809       
 810       FORTRAN(dgesv,(nev,&nrhs,lhs,nev,ipiv,dbj,nev,&info));
 811       
 812       /* check the size of dbj */
 813       
 814       bjmax=0.;
 815       dbjmaxOLD2=dbjmaxOLD1;
 816       dbjmaxOLD1=dbjmax;
 817       dbjmax=0.;
 818       for(i=0;i<*nev;i++){
 819         if(fabs(bj[i])>bjmax) bjmax=fabs(bj[i]);
 820         if(fabs(dbj[i])>dbjmax) dbjmax=fabs(dbj[i]);
 821       }
 822 
 823       iconvergence=0;
 824       idivergence=0;
 825 
 826       if(dbjmax<=0.005*bjmax){
 827         
 828         //calculate bjp: the derivative of bj w.r.t. time
 829         
 830         for(j=0; j<*nev; j++){
 831           bjp[j]=bjbasp[j]+bjincp[j]*aaa[j];
 832         }
 833         FORTRAN(dgetrs,("No transpose",nev,&nrhs,lhs,nev,ipiv,bjp,nev,&info));
 834         iconvergence=1;    
 835       }
 836       else{
 837           if(loop>=i0 && loop<=ic){
 838 
 839           /* check for divergence */
 840 
 841           if((dbjmax>dbjmaxOLD1) && (dbjmax>dbjmaxOLD2)){
 842 
 843             /* divergence --> cutback */ 
 844 
 845             idivergence=1;
 846             icutb++;
 847             break;
 848           }
 849         }
 850         else{
 851           if(loop>ic){
 852 
 853             /* cutback after ic iterations*/
 854 
 855             idivergence=1;
 856             icutb++;
 857             break;
 858           }
 859         }
 860       }
 861 
 862       /* add dbj to db */
 863       
 864       for(j=0;j<*nev;j++){
 865         bj[j]+=dbj[j];
 866       }
 867       
 868     }while(1);
 869   }while(idivergence==1 && icutb<10 && *idrct==0);
 870 
 871   if(icutb>=10){
 872 
 873     //no convergence, stop all
 874 
 875     printf("*ERROR: Contact did not converge.\n");
 876     FORTRAN(stop,());
 877   }
 878   
 879   /* convergence has been reached */
 880   /* restoring aa[(*iinc-1)*nev+i] */
 881 
 882   for(i=0;i<*nev;i++){
 883     aa[i]+=bb[i]*(*time-*dtime);
 884   }
 885   
 886   /* calculating the linearized force function connecting the
 887      mechanical+contact load at the start of the increment with
 888      the mechanical+contact load at the end of the increment */
 889   
 890   if(!(*cyclicsymmetry)){
 891       for(i=0;i<*nev;i++){
 892           i2=(long long)i*neq[1];
 893           
 894           aanew[i]=aamech[i];
 895           for(j=0;j<*nactcont;j++){
 896               aanew[i]+=z[i2+ikactcont[j]]*bcont[ikactcont[j]];
 897           }
 898           
 899           bb[i]=(aanew[i]-aa[i])/(*dtime);
 900           aa[i]=aanew[i]-bb[i]**time;
 901       }
 902   }else{
 903       memcpy(&aanew[0],&aamech[0],sizeof(double)**nev);
 904       for(j=0;j<*nactcont;j++){
 905           FORTRAN(nident,(izdof,&ikactcont[j],nzdof,&id));
 906           if(id!=0){
 907               if(izdof[id-1]==ikactcont[j]){
 908                   for(i=0;i<*nev;i++){
 909                       aanew[i]+=z[i**nzdof+id-1]*bcont[ikactcont[j]];
 910                   }
 911               }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 912           }else{printf("*ERROR in dynacont\n");FORTRAN(stop,());}
 913       }
 914       for(i=0;i<*nev;i++){
 915           bb[i]=(aanew[i]-aa[i])/(*dtime);
 916           aa[i]=aanew[i]-bb[i]**time;
 917       }
 918   }
 919   
 920   SFREE(aaa);SFREE(bbb);SFREE(bjbas);SFREE(bjinc);SFREE(dbj);SFREE(lhs);
 921   SFREE(ipiv);SFREE(bjbasp);SFREE(bjincp);SFREE(dbjp);SFREE(ilactcont);
 922   SFREE(dbcont);
 923   SFREE(xforcdiff);SFREE(xloaddiff);SFREE(xboundiff),SFREE(xbodydiff);
 924 
 925   if(*ithermal==1) SFREE(t1diff);
 926   
 927   *ikactcontp=ikactcont;*ikactmechp=ikactmech;
 928   
 929   return;
 930 }
 931 
 932 

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