root/src/linstatic.c

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

DEFINITIONS

This source file includes following definitions.
  1. linstatic

   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 void linstatic(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon,
  36              ITG *ne, 
  37              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  38              ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
  39              ITG *nmpc, 
  40              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  41              ITG *nelemload, char *sideload, double *xload,
  42              ITG *nload, ITG *nactdof, 
  43              ITG **icolp, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, 
  44              ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, 
  45              ITG *ilboun,
  46              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  47              double *alcon, ITG *nalcon, double *alzero, ITG *ielmat,
  48              ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_,
  49              double *t0, double *t1, double *t1old,
  50              ITG *ithermal,double *prestr, ITG *iprestr, 
  51              double *vold,ITG *iperturb, double *sti, ITG *nzs,  
  52              ITG *kode, char *filab, double *eme,
  53              ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon,
  54              ITG *nplkcon,
  55              double *xstate, ITG *npmat_, char *matname, ITG *isolver,
  56              ITG *mi, ITG *ncmat_, ITG *nstate_, double *cs, ITG *mcs,
  57              ITG *nkon, double *ener, double *xbounold,
  58              double *xforcold, double *xloadold,
  59              char *amname, double *amta, ITG *namta,
  60              ITG *nam, ITG *iamforc, ITG *iamload,
  61              ITG *iamt1, ITG *iamboun, double *ttime, char *output, 
  62              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, double *timepar,
  67              double *thicke, char *jobnamec,char *tieset,ITG *ntie,
  68                ITG *istep,ITG *nmat,ITG *ielprop,double *prop,char *typeboun){
  69   
  70   char description[13]="            ";
  71 
  72   ITG *inum=NULL,k,*icol=NULL,*irow=NULL,ielas,icmd=0,iinc=1,nasym=0,i,j,ic,ir,
  73       mass[2]={0,0}, stiffness=1, buckling=0, rhsi=1, intscheme=0,*ncocon=NULL,
  74       *nshcon=NULL,mode=-1,noddiam=-1,*ipobody=NULL,inewton=0,coriolis=0,iout,
  75       ifreebody,*itg=NULL,ntg=0,symmetryflag=0,inputformat=0,ngraph=1,im,
  76       mt=mi[1]+1,ne0,*integerglob=NULL,iglob=0,*ipneigh=NULL,*neigh=NULL,
  77       icfd=0,*inomat=NULL,mortar,*islavact=NULL,*islavnode=NULL,*nslavnode=NULL,
  78       *islavsurf=NULL,nretain,*iretain=NULL,*noderetain=NULL,*ndirretain=NULL,
  79       nmethodl;
  80 
  81   double *stn=NULL,*v=NULL,*een=NULL,cam[5],*xstiff=NULL,*stiini=NULL,*tper,
  82          *f=NULL,*fn=NULL,qa[3],*fext=NULL,*epn=NULL,*xstateini=NULL,
  83          *vini=NULL,*stx=NULL,*enern=NULL,*xbounact=NULL,*xforcact=NULL,
  84          *xloadact=NULL,*t1act=NULL,*ampli=NULL,*xstaten=NULL,*eei=NULL,
  85          *enerini=NULL,*cocon=NULL,*shcon=NULL,*physcon=NULL,*qfx=NULL,
  86          *qfn=NULL,sigma=0.,*cgr=NULL,*xbodyact=NULL,*vr=NULL,*vi=NULL,
  87          *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,*springarea=NULL,
  88          *eenmax=NULL,*fnr=NULL,*fni=NULL,*emn=NULL,*clearini=NULL,ptime,
  89          *emeini=NULL,*doubleglob=NULL,*au=NULL,*ad=NULL,*b=NULL,*aub=NULL,
  90          *adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,*cdn=NULL,*cdnr=NULL,
  91          *cdni=NULL,*submatrix=NULL;
  92 
  93 #ifdef SGI
  94   ITG token;
  95 #endif
  96 
  97   /* dummy arguments for the results call */
  98 
  99   double *veold=NULL,*accold=NULL,bet,gam,dtime,time,reltime=1.;
 100 
 101   icol=*icolp;
 102   irow=*irowp;
 103 
 104   tper=&timepar[1];
 105 
 106   time=*tper;
 107   dtime=*tper;
 108 
 109   ne0=*ne;
 110 
 111   /* determining the global values to be used as boundary conditions
 112      for a submodel */
 113 
 114   getglobalresults(jobnamec,&integerglob,&doubleglob,nboun,iamboun,xboun,
 115                    nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
 116                    ithermal,nk,t1,iamt1);
 117 
 118   /* allocating fields for the actual external loading */
 119 
 120   NNEW(xbounact,double,*nboun);
 121   for(k=0;k<*nboun;++k){xbounact[k]=xbounold[k];}
 122   NNEW(xforcact,double,*nforc);
 123   NNEW(xloadact,double,2**nload);
 124   NNEW(xbodyact,double,7**nbody);
 125   /* copying the rotation axis and/or acceleration vector */
 126   for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
 127   if(*ithermal==1){
 128     NNEW(t1act,double,*nk);
 129     for(k=0;k<*nk;++k){t1act[k]=t1old[k];}
 130   }
 131   
 132   /* assigning the body forces to the elements */ 
 133 
 134   if(*nbody>0){
 135       ifreebody=*ne+1;
 136       NNEW(ipobody,ITG,2*ifreebody**nbody);
 137       for(k=1;k<=*nbody;k++){
 138           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 139                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 140           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 141       }
 142       RENEW(ipobody,ITG,2*(ifreebody-1));
 143   }
 144 
 145   /* allocating a field for the instantaneous amplitude */
 146 
 147   NNEW(ampli,double,*nam);
 148 
 149   FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
 150               xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
 151               t1old,t1,t1act,iamt1,nk,amta,
 152               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
 153               xbounold,xboun,xbounact,iamboun,nboun,
 154               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 155               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 156               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 157               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 158 
 159   /* determining the internal forces and the stiffness coefficients */
 160 
 161   NNEW(f,double,*neq);
 162 
 163   /* allocating a field for the stiffness matrix */
 164 
 165   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 166 
 167   iout=-1;
 168   NNEW(v,double,mt**nk);
 169   NNEW(fn,double,mt**nk);
 170   NNEW(stx,double,6*mi[0]**ne);
 171   NNEW(inum,ITG,*nk);
 172   results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 173           elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 174           ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
 175           prestr,iprestr,filab,eme,emn,een,iperturb,
 176           f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 177           ndirboun,xbounact,nboun,ipompc,
 178           nodempc,coefmpc,labmpc,nmpc,nmethod,cam,neq,veold,accold,
 179           &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 180           xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 181           &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 182           emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 183           iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 184           fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 185           &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 186           sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 187           &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 188           islavsurf,ielprop,prop);
 189   SFREE(v);SFREE(fn);SFREE(stx);SFREE(inum);
 190   iout=1;
 191   
 192   /* determining the system matrix and the external forces */
 193 
 194   NNEW(ad,double,*neq);
 195   NNEW(fext,double,*neq);
 196 
 197   if(*nmethod==11){
 198       NNEW(au,double,nzs[2]);
 199       rhsi=0;
 200       nmethodl=2;
 201   }else{
 202       NNEW(au,double,*nzs);
 203       nmethodl=*nmethod;
 204   }
 205       
 206 
 207   FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounact,nboun,
 208             ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
 209             nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
 210             nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,&nmethodl,
 211             ikmpc,ilmpc,ikboun,ilboun,
 212             elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 213             ielorien,norien,orab,ntmat_,
 214             t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
 215             nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 216             xstiff,npmat_,&dtime,matname,mi,
 217             ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,physcon,
 218             shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,&coriolis,
 219             ibody,xloadold,&reltime,veold,springarea,nstate_,
 220             xstateini,xstate,thicke,integerglob,doubleglob,
 221             tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 222             pmastsurf,&mortar,clearini,ielprop,prop));
 223 
 224   /* determining the right hand side */
 225 
 226   NNEW(b,double,*neq);
 227   for(k=0;k<*neq;++k){
 228       b[k]=fext[k]-f[k];
 229   }
 230   SFREE(fext);SFREE(f);
 231 
 232   /* generation of a substructure stiffness matrix */
 233 
 234   if(*nmethod==11){
 235 
 236       /* factorizing the matrix */
 237 
 238       if(*neq>0){
 239           if(*isolver==0){
 240 #ifdef SPOOLES
 241               spooles_factor(ad,au,adb,aub,&sigma,icol,irow,neq,nzs,&symmetryflag,
 242                              &inputformat,&nzs[2]);
 243 #else
 244               printf("*ERROR in linstatic: the SPOOLES library is not linked\n\n");
 245               FORTRAN(stop,());
 246 #endif
 247           }
 248           else if(*isolver==7){
 249 #ifdef PARDISO
 250               pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,neq,nzs,
 251                              &symmetryflag,&inputformat,jq,&nzs[2]);
 252 #else
 253               printf("*ERROR in linstatic: the PARDISO library is not linked\n\n");
 254               FORTRAN(stop,());
 255 #endif
 256           }
 257       }
 258       
 259     /* determining the nodes and the degrees of freedom in those nodes
 260        belonging to the substructure */
 261       
 262     NNEW(iretain,ITG,*nk);
 263     NNEW(noderetain,ITG,*nk);
 264     NNEW(ndirretain,ITG,*nk);
 265     nretain=0;
 266 
 267     for(i=0;i<*nboun;i++){
 268         if(strcmp1(&typeboun[i],"C")==0){
 269             iretain[nretain]=i+1;
 270             noderetain[nretain]=nodeboun[i];
 271             ndirretain[nretain]=ndirboun[i];
 272             nretain++;
 273         }
 274     }
 275 
 276     RENEW(iretain,ITG,nretain);
 277     RENEW(noderetain,ITG,nretain);
 278     RENEW(ndirretain,ITG,nretain);
 279 
 280     /* solving the system of equations with appropriate rhs */
 281 
 282     NNEW(submatrix,double,nretain*nretain);
 283 
 284     for(i=0;i<nretain;i++){
 285         DMEMSET(b,0,*neq,0.);
 286         ic=*neq+iretain[i]-1;
 287         for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
 288             ir=irow[j]-1;
 289             b[ir]-=au[j];
 290         }
 291 
 292         /* solving the system */
 293           
 294         if(*neq>0){
 295             if(*isolver==0){
 296 #ifdef SPOOLES
 297                 spooles_solve(b,neq);
 298 #endif
 299             }
 300             else if(*isolver==7){
 301 #ifdef PARDISO
 302                 pardiso_solve(b,neq,&symmetryflag);
 303 #endif
 304             
 305             }
 306         }
 307 
 308           /* calculating the internal forces */
 309 
 310           NNEW(v,double,mt**nk);
 311           NNEW(fn,double,mt**nk);
 312           NNEW(stn,double,6**nk);
 313           NNEW(inum,ITG,*nk);
 314           NNEW(stx,double,6*mi[0]**ne);
 315           
 316           if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 317           if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 318           if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 319           
 320           NNEW(eei,double,6*mi[0]**ne);
 321           if(*nener==1){
 322               NNEW(stiini,double,6*mi[0]**ne);
 323               NNEW(enerini,double,mi[0]**ne);}
 324           
 325           /* replacing the appropriate boundary value by unity */
 326 
 327           xbounact[iretain[i]-1]=1.;
 328 
 329           results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 330             elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 331             ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
 332             prestr,iprestr,filab,eme,emn,een,iperturb,
 333             f,fn,nactdof,&iout,qa,vold,b,nodeboun,ndirboun,xbounact,nboun,ipompc,
 334             nodempc,coefmpc,labmpc,nmpc,nmethod,cam,neq,veold,accold,&bet,
 335             &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 336             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 337             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 338             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 339             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 340             nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
 341             &ne0,xforc,nforc,thicke,shcon,nshcon,
 342             sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 343             &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 344             islavsurf,ielprop,prop);
 345 
 346           xbounact[iretain[i]-1]=0.;
 347           
 348           SFREE(v);SFREE(stn);SFREE(inum);SFREE(stx);
 349           
 350           if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
 351           if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
 352           if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
 353         
 354           SFREE(eei);if(*nener==1){SFREE(stiini);SFREE(enerini);}
 355 
 356           /* storing the internal forces in the substructure
 357              stiffness matrix */
 358 
 359           for(j=0;j<nretain;j++){
 360               submatrix[i*nretain+j]=fn[mt*(noderetain[j]-1)+ndirretain[j]];
 361           }
 362 
 363           SFREE(fn);
 364 
 365     }
 366 
 367     SFREE(au);SFREE(ad);SFREE(b);
 368     SFREE(iretain);
 369     
 370     SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(t1act);SFREE(ampli);
 371     SFREE(xbodyact);if(*nbody>0) SFREE(ipobody);SFREE(xstiff);
 372     
 373     if(iglob==1){SFREE(integerglob);SFREE(doubleglob);}
 374     
 375     FORTRAN(writesubmatrix,(submatrix,noderetain,ndirretain,&nretain,jobnamec));
 376 
 377     SFREE(submatrix);SFREE(noderetain);SFREE(ndirretain);
 378     
 379     return;
 380 
 381 
 382   }else if(*nmethod!=0){
 383 
 384     if(*isolver==0){
 385 #ifdef SPOOLES
 386       spooles(ad,au,adb,aub,&sigma,b,icol,irow,neq,nzs,&symmetryflag,
 387               &inputformat,&nzs[2]);
 388 #else
 389             printf("*ERROR in linstatic: the SPOOLES library is not linked\n\n");
 390             FORTRAN(stop,());
 391 #endif
 392     }
 393     else if((*isolver==2)||(*isolver==3)){
 394       preiter(ad,&au,b,&icol,&irow,neq,nzs,isolver,iperturb);
 395     }
 396     else if(*isolver==4){
 397 #ifdef SGI
 398       token=1;
 399       sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,neq,nzs,token);
 400 #else
 401             printf("*ERROR in linstatic: the SGI library is not linked\n\n");
 402             FORTRAN(stop,());
 403 #endif
 404     }
 405     else if(*isolver==5){
 406 #ifdef TAUCS
 407       tau(ad,&au,adb,aub,&sigma,b,icol,&irow,neq,nzs);
 408 #else
 409             printf("*ERROR in linstatic: the TAUCS library is not linked\n\n");
 410             FORTRAN(stop,());
 411 #endif
 412     }
 413     else if(*isolver==7){
 414 #ifdef PARDISO
 415       pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,neq,nzs,
 416                    &symmetryflag,&inputformat,jq,&nzs[2]);
 417 #else
 418             printf("*ERROR in linstatic: the PARDISO library is not linked\n\n");
 419             FORTRAN(stop,());
 420 #endif
 421     }
 422 
 423     SFREE(ad);SFREE(au);
 424 
 425     /* calculating the displacements and the stresses and storing */
 426     /* the results in frd format for each valid eigenmode */
 427 
 428     NNEW(v,double,mt**nk);
 429     NNEW(fn,double,mt**nk);
 430     NNEW(stn,double,6**nk);
 431     NNEW(inum,ITG,*nk);
 432     NNEW(stx,double,6*mi[0]**ne);
 433   
 434     if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 435     if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 436     if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 437 
 438     NNEW(eei,double,6*mi[0]**ne);
 439     if(*nener==1){
 440         NNEW(stiini,double,6*mi[0]**ne);
 441         NNEW(enerini,double,mi[0]**ne);}
 442 
 443     results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 444             elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 445             ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
 446             prestr,iprestr,filab,eme,emn,een,iperturb,
 447             f,fn,nactdof,&iout,qa,vold,b,nodeboun,ndirboun,xbounact,nboun,ipompc,
 448             nodempc,coefmpc,labmpc,nmpc,nmethod,cam,neq,veold,accold,&bet,
 449             &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 450             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 451             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 452             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 453             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 454             nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
 455             &ne0,xforc,nforc,thicke,shcon,nshcon,
 456             sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 457             &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 458             islavsurf,ielprop,prop);
 459 
 460     SFREE(eei);
 461     if(*nener==1){
 462         SFREE(stiini);SFREE(enerini);}
 463 
 464     memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
 465     memcpy(&sti[0],&stx[0],sizeof(double)*6*mi[0]**ne);
 466 /*    for(k=0;k<mt**nk;++k){
 467       vold[k]=v[k];
 468     }
 469     for(k=0;k<6*mi[0]**ne;++k){
 470       sti[k]=stx[k];
 471       }*/
 472 
 473     ++*kode;
 474 
 475     /* for cyclic symmetric sectors: duplicating the results */
 476 
 477     if(*mcs>0){
 478         ptime=*ttime+time;
 479       frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,t1act,
 480                    fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
 481                    nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,
 482                    qfn,ialset,istartset,iendset,trab,inotr,ntrans,orab,
 483                    ielorien,norien,sti,veold,&noddiam,set,nset,emn,thicke,
 484                    jobnamec,&ne0,cdn,&mortar,nmat);
 485     }
 486     else{
 487         if(strcmp1(&filab[1044],"ZZS")==0){
 488             NNEW(neigh,ITG,40**ne);
 489             NNEW(ipneigh,ITG,*nk);
 490         }
 491         ptime=*ttime+time;
 492         frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 493             kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
 494             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
 495             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 496             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 497             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 498             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 499         if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 500     }
 501 
 502     SFREE(v);SFREE(stn);SFREE(inum);
 503     SFREE(b);SFREE(stx);SFREE(fn);
 504 
 505     if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
 506     if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
 507     if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
 508 
 509   }
 510   else {
 511 
 512     /* error occurred in mafill: storing the geometry in frd format */
 513 
 514     ++*kode;
 515     NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 516     if(strcmp1(&filab[1044],"ZZS")==0){
 517         NNEW(neigh,ITG,40**ne);
 518         NNEW(ipneigh,ITG,*nk);
 519     }
 520     ptime=*ttime+time;
 521     frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 522             kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
 523             nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
 524             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 525             mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 526             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 527             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 528     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 529     SFREE(inum);FORTRAN(stop,());
 530 
 531   }
 532 
 533   /* updating the loading at the end of the step; 
 534      important in case the amplitude at the end of the step
 535      is not equal to one */
 536 
 537   for(k=0;k<*nboun;++k){xbounold[k]=xbounact[k];}
 538   for(k=0;k<*nforc;++k){xforcold[k]=xforcact[k];}
 539   for(k=0;k<2**nload;++k){xloadold[k]=xloadact[k];}
 540   for(k=0;k<7**nbody;k=k+7){xbodyold[k]=xbodyact[k];}
 541   if(*ithermal==1){
 542     for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
 543     for(k=0;k<*nk;++k){vold[mt*k]=t1act[k];}
 544   }
 545 
 546   SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(t1act);SFREE(ampli);
 547   SFREE(xbodyact);if(*nbody>0) SFREE(ipobody);SFREE(xstiff);
 548 
 549   if(iglob==1){SFREE(integerglob);SFREE(doubleglob);}
 550 
 551   *icolp=icol;
 552   *irowp=irow;
 553 
 554   (*ttime)+=(*tper);
 555  
 556   return;
 557 }

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