root/src/arpackbu.c

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

DEFINITIONS

This source file includes following definitions.
  1. arpackbu

   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 #ifdef ARPACK
  19 
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include "CalculiX.h"
  24 #ifdef SPOOLES
  25    #include "spooles.h"
  26 #endif
  27 #ifdef SGI
  28    #include "sgi.h"
  29 #endif
  30 #ifdef TAUCS
  31    #include "tau.h"
  32 #endif
  33 #ifdef PARDISO
  34    #include "pardiso.h"
  35 #endif
  36 
  37 void arpackbu(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon,
  38              ITG *ne, 
  39              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  40              ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
  41              ITG *nmpc, 
  42              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  43              ITG *nelemload, char *sideload, double *xload,
  44              ITG *nload, 
  45              ITG *nactdof, 
  46              ITG *icol, ITG *jq, ITG *irow, ITG *neq, ITG *nzl, 
  47              ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, 
  48              ITG *ilboun,
  49              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  50              double *alcon, ITG *nalcon, double *alzero, ITG *ielmat,
  51              ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_,
  52              double *t0, double *t1, double *t1old, 
  53              ITG *ithermal,double *prestr, ITG *iprestr, 
  54              double *vold,ITG *iperturb, double *sti, ITG *nzs,  
  55              ITG *kode, ITG *mei, double *fei,
  56              char *filab, double *eme,
  57              ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon,
  58              ITG *nplkcon,
  59              double *xstate, ITG *npmat_, char *matname, ITG *mi,
  60              ITG *ncmat_, ITG *nstate_, double *ener, char *output, 
  61              char *set, ITG *nset, ITG *istartset,
  62              ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  63              char *prset, ITG *nener, ITG *isolver, double *trab, 
  64              ITG *inotr, ITG *ntrans, double *ttime,double *fmpc,
  65              char *cbody, ITG *ibody,double *xbody, ITG *nbody, 
  66              double *thicke,char *jobnamec,ITG *nmat,ITG *ielprop,
  67              double *prop){
  68   
  69   char bmat[2]="G", which[3]="LM", howmny[2]="A",
  70       description[13]="            ",*tieset=NULL;
  71 
  72   ITG *inum=NULL,k,ido,dz,iparam[11],ipntr[11],lworkl,im,nasym=0,
  73     info,rvec=1,*select=NULL,lfin,j,lint,iout,iconverged=0,ielas,icmd=0,
  74     iinc=1,istep=1,*ncocon=NULL,*nshcon=NULL,nev,ncv,mxiter,jrow,
  75     *ipobody=NULL,inewton=0,coriolis=0,ifreebody,symmetryflag=0,
  76     inputformat=0,ngraph=1,mt=mi[1]+1,mass[2]={0,0}, stiffness=1, buckling=0, 
  77     rhsi=1, intscheme=0, noddiam=-1,*ipneigh=NULL,*neigh=NULL,ne0,
  78     *integerglob=NULL,ntie,icfd=0,*inomat=NULL,mortar=0,*islavnode=NULL,
  79     *islavact=NULL,*nslavnode=NULL,*islavsurf=NULL;
  80 
  81   double *stn=NULL,*v=NULL,*resid=NULL,*z=NULL,*workd=NULL,
  82     *workl=NULL,*d=NULL,sigma,*temp_array=NULL,
  83     *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[3],*fext=NULL,
  84     time=0.,*epn=NULL,*fnr=NULL,*fni=NULL,*emn=NULL,*cdn=NULL,
  85     *xstateini=NULL,*xstiff=NULL,*stiini=NULL,*vini=NULL,*stx=NULL,
  86     *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*cocon=NULL,
  87     *shcon=NULL,*physcon=NULL,*qfx=NULL,*qfn=NULL,tol, *cgr=NULL,
  88     *xloadold=NULL,reltime,*vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,
  89     *vmax=NULL,*stnmax=NULL,*cs=NULL,*springarea=NULL,*eenmax=NULL,
  90     *emeini=NULL,*doubleglob=NULL,*au=NULL,*clearini=NULL,
  91     *ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
  92     *cdnr=NULL,*cdni=NULL;
  93 
  94   /* buckling routine; only for mechanical applications */
  95 
  96   /* dummy arguments for the results call */
  97 
  98   double *veold=NULL,*accold=NULL,bet,gam,dtime;
  99 
 100 #ifdef SGI
 101   ITG token;
 102 #endif
 103  
 104   /* copying the frequency parameters */
 105 
 106   nev=mei[0];
 107   ncv=mei[1];
 108   mxiter=mei[2];
 109   tol=fei[0];
 110 
 111   /* calculating the stresses due to the buckling load; this is a second
 112      order calculation if iperturb != 0 */
 113 
 114   *nmethod=1;
 115   
 116   /* assigning the body forces to the elements */ 
 117 
 118   if(*nbody>0){
 119       ifreebody=*ne+1;
 120       NNEW(ipobody,ITG,2*ifreebody**nbody);
 121       for(k=1;k<=*nbody;k++){
 122           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 123                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 124           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 125       }
 126       RENEW(ipobody,ITG,2*(ifreebody-1));
 127   }
 128 
 129   /* determining the internal forces and the stiffness coefficients */
 130 
 131   NNEW(f,double,neq[0]);
 132 
 133   /* allocating a field for the stiffness matrix */
 134 
 135   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 136 
 137 //  iout=-1;
 138   NNEW(v,double,mt**nk);
 139   NNEW(fn,double,mt**nk);
 140   NNEW(stx,double,6*mi[0]**ne);
 141 
 142   iout=-1;
 143   NNEW(inum,ITG,*nk);
 144   if(*iperturb==0){
 145      results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 146              elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 147              ielorien,norien,orab,ntmat_,t0,t0,ithermal,
 148              prestr,iprestr,filab,eme,emn,een,iperturb,
 149              f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 150              ndirboun,xboun,nboun,ipompc,
 151              nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,
 152              &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 153              xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 154              &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 155              emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 156              iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 157              fmpc,nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,
 158              &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 159              sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 160              &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 161              islavsurf,ielprop,prop);
 162   }else{
 163      results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 164              elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 165              ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
 166              prestr,iprestr,filab,eme,emn,een,iperturb,
 167              f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 168              ndirboun,xboun,nboun,ipompc,
 169              nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,
 170              &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 171              xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 172              &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 173              emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 174              iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 175              fmpc,nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,
 176              &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 177              sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 178              &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 179              islavsurf,ielprop,prop);
 180   }
 181 
 182   SFREE(v);SFREE(fn);SFREE(stx);SFREE(inum);
 183   iout=1;
 184 
 185   /* determining the system matrix and the external forces */
 186 
 187   NNEW(ad,double,neq[0]);
 188   NNEW(au,double,nzs[0]);
 189   NNEW(fext,double,neq[0]);
 190 
 191   if(*iperturb==0){
 192     FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 193               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 194               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 195               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 196               ikmpc,ilmpc,ikboun,ilboun,
 197               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 198               ielorien,norien,orab,ntmat_,
 199               t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 200               &nzs[0],stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 201               xstiff,npmat_,&dtime,matname,mi,
 202               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,physcon,
 203               shcon,nshcon,cocon,ncocon,ttime,&time,&istep,&iinc,&coriolis,
 204               ibody,xloadold,&reltime,veold,springarea,nstate_,
 205               xstateini,xstate,thicke,integerglob,doubleglob,
 206               tieset,istartset,iendset,ialset,&ntie,&nasym,pslavsurf,pmastsurf,
 207               &mortar,clearini,ielprop,prop));
 208   }
 209   else{
 210     FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 211               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 212               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 213               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 214               ikmpc,ilmpc,ikboun,ilboun,
 215               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 216               ielorien,norien,orab,ntmat_,
 217               t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
 218               &nzs[0],stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 219               xstiff,npmat_,&dtime,matname,mi,
 220               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,physcon,
 221               shcon,nshcon,cocon,ncocon,ttime,&time,&istep,&iinc,&coriolis,
 222               ibody,xloadold,&reltime,veold,springarea,nstate_,
 223               xstateini,xstate,thicke,integerglob,doubleglob,
 224               tieset,istartset,iendset,ialset,&ntie,&nasym,pslavsurf,
 225               pmastsurf,&mortar,clearini,ielprop,prop));
 226   }
 227 
 228   /* determining the right hand side */
 229 
 230   NNEW(b,double,neq[0]);
 231   for(k=0;k<neq[0];++k){
 232       b[k]=fext[k]-f[k];
 233   }
 234   SFREE(fext);SFREE(f);
 235 
 236   if(*nmethod==0){
 237 
 238     /* error occurred in mafill: storing the geometry in frd format */
 239 
 240     ++*kode;
 241     NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 242     if(strcmp1(&filab[1044],"ZZS")==0){
 243         NNEW(neigh,ITG,40**ne);
 244         NNEW(ipneigh,ITG,*nk);
 245     }
 246 
 247     frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 248             kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
 249             nstate_,&istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
 250             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 251             mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 252             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 253             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 254     
 255     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 256     SFREE(inum);FORTRAN(stop,());
 257 
 258   }
 259 
 260   *nmethod=3;
 261   buckling=1;rhsi=0;
 262 
 263   sigma=0.;
 264   if(*isolver==0){
 265 #ifdef SPOOLES
 266     spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],&symmetryflag,
 267             &inputformat,&nzs[2]);
 268 #else
 269     printf("*ERROR in arpackbu: the SPOOLES library is not linked\n\n");
 270     FORTRAN(stop,());
 271 #endif
 272   }
 273   else if(*isolver==4){
 274 #ifdef SGI
 275     token=1;
 276     sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
 277 #else
 278     printf("*ERROR in arpackbu: the SGI library is not linked\n\n");
 279     FORTRAN(stop,());
 280 #endif
 281   }
 282   else if(*isolver==5){
 283 #ifdef TAUCS
 284     tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[0],&nzs[0]);
 285 #else
 286     printf("*ERROR in arpackbu: the TAUCS library is not linked\n\n");
 287     FORTRAN(stop,());
 288 #endif
 289   }
 290   else if(*isolver==7){
 291 #ifdef PARDISO
 292     pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
 293                  &symmetryflag,&inputformat,jq,&nzs[2]);
 294 #else
 295     printf("*ERROR in arpackbu: the PARDISO library is not linked\n\n");
 296     FORTRAN(stop,());
 297 #endif
 298   }
 299 
 300   /* calculating the displacements and the stresses and storing */
 301   /* the results in frd format for each valid eigenmode */
 302 
 303   NNEW(v,double,mt**nk);
 304   NNEW(fn,double,mt**nk);
 305   NNEW(stn,double,6**nk);
 306   NNEW(inum,ITG,*nk);
 307   NNEW(stx,double,6*mi[0]**ne);
 308   
 309   if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 310   if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 311   if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 312 
 313   NNEW(eei,double,6*mi[0]**ne);
 314   if(*nener==1){
 315       NNEW(stiini,double,6*mi[0]**ne);
 316       NNEW(enerini,double,mi[0]**ne);}
 317 
 318   if(*iperturb==0){
 319     results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 320             stx,elcon,nelcon,
 321             rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,
 322             ntmat_,t0,t0,ithermal,
 323             prestr,iprestr,filab,eme,emn,een,iperturb,
 324             f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 325             ndirboun,xboun,nboun,ipompc,
 326             nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,&bet,
 327             &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 328             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 329             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 330             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 331             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 332             nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
 333             &ne0,xforc,nforc,thicke,shcon,nshcon,
 334             sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 335             &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 336             islavsurf,ielprop,prop);}
 337   else{
 338     results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 339             stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 340             ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
 341             prestr,iprestr,filab,eme,emn,een,iperturb,
 342             f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 343             ndirboun,xboun,nboun,ipompc,
 344             nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,&bet,
 345             &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 346             xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 347             ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 348             xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 349             ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 350             nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
 351             &ne0,xforc,nforc,thicke,shcon,nshcon,
 352             sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 353             &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 354             islavsurf,ielprop,prop);
 355   }
 356 
 357   for(k=0;k<mt**nk;++k){
 358     vold[k]=v[k];
 359   }
 360 
 361   ++*kode;
 362   if(strcmp1(&filab[1044],"ZZS")==0){
 363       NNEW(neigh,ITG,40**ne);
 364       NNEW(ipneigh,ITG,*nk);
 365   }
 366 
 367   frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 368             kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
 369             nstate_,&istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
 370             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 371             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 372             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 373             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 374 
 375   if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 376   SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);
 377 
 378   if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
 379   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
 380   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
 381 
 382   /* in buckling mode stx and sti are kept */
 383 
 384 
 385   /* calculation of the left hand matrix (ad and au) and the right
 386      hand matrix (adb and aub); stx are the stresses due to the buckling
 387      load, sti due to previous loads, if any */
 388 
 389   NNEW(aub,double,nzs[0]);
 390   NNEW(adb,double,neq[0]);
 391 
 392   NNEW(fext,double,neq[0]);
 393 
 394   if(*iperturb==0){
 395     FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 396               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 397               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 398               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 399               ikmpc,ilmpc,ikboun,ilboun,
 400               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 401               ielorien,norien,orab,ntmat_,
 402               t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 403               &nzs[0],stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 404               xstiff,npmat_,&dtime,matname,mi,
 405               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,physcon,
 406               shcon,nshcon,cocon,ncocon,ttime,&time,&istep,&iinc,&coriolis,
 407               ibody,xloadold,&reltime,veold,springarea,nstate_,
 408               xstateini,xstate,thicke,integerglob,doubleglob,
 409               tieset,istartset,iendset,ialset,&ntie,&nasym,pslavsurf,
 410               pmastsurf,&mortar,clearini,ielprop,prop));
 411   }
 412   else{
 413     FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 414               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 415               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 416               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 417               ikmpc,ilmpc,ikboun,ilboun,
 418               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 419               ielorien,norien,orab,ntmat_,
 420               t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
 421               &nzs[0],stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 422               xstiff,npmat_,&dtime,matname,mi,
 423               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,physcon,
 424               shcon,nshcon,cocon,ncocon,ttime,&time,&istep,&iinc,&coriolis,
 425               ibody,xloadold,&reltime,veold,springarea,nstate_,
 426               xstateini,xstate,thicke,integerglob,doubleglob,
 427               tieset,istartset,iendset,ialset,&ntie,&nasym,pslavsurf,
 428               pmastsurf,&mortar,clearini,ielprop,prop));
 429   }
 430 
 431   SFREE(stx);SFREE(fext);if(*nbody>0) SFREE(ipobody);
 432 
 433   if(*nmethod==1){return;}
 434 
 435   /* loop checking the plausibility of the buckling factor
 436      if (5*sigma<buckling factor<50000*sigma) the solution is accepted,
 437      else sigma is set to buckling factor/500, and a new iteration is
 438      started */
 439 
 440   sigma=1.;
 441 
 442   do{
 443 
 444 
 445   /* LU decomposition of the left hand matrix */
 446 
 447 //  sigma=1.;
 448   if(*isolver==0){
 449 #ifdef SPOOLES
 450     spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[0],&nzs[0],
 451                    &symmetryflag,&inputformat,&nzs[2]);
 452 #endif
 453   }
 454   else if(*isolver==4){
 455 #ifdef SGI
 456     token=2;
 457     sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[0],&nzs[0],token);
 458 #endif
 459   }
 460   else if(*isolver==5){
 461 #ifdef TAUCS
 462     tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[0],&nzs[0]);
 463 #endif
 464   }
 465   else if(*isolver==7){
 466 #ifdef PARDISO
 467     pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[0],&nzs[0],
 468                    &symmetryflag,&inputformat,jq,&nzs[2]);
 469 #endif
 470   }
 471 
 472   /* calculating the bucking factors and buckling modes */
 473 
 474   printf(" Calculating the buckling factors and buckling modes:\n\n");
 475 
 476   ido=0;
 477   dz=neq[0];
 478   iparam[0]=1;
 479   iparam[2]=mxiter;
 480   iparam[3]=1;
 481   iparam[6]=4;
 482 
 483   lworkl=ncv*(8+ncv);
 484   info=0;
 485 
 486   NNEW(resid,double,neq[0]);
 487   NNEW(z,double,ncv*neq[0]);
 488   NNEW(workd,double,3*neq[0]);
 489   NNEW(workl,double,lworkl);
 490 
 491   FORTRAN(dsaupd,(&ido,bmat,&neq[0],which,&nev,&tol,resid,&ncv,z,&dz,iparam,ipntr,workd,
 492           workl,&lworkl,&info));
 493 
 494   NNEW(temp_array,double,neq[0]);
 495 
 496   while((ido==-1)||(ido==1)||(ido==2)){
 497     if(ido==-1){
 498       FORTRAN(op,(&neq[0],&workd[ipntr[0]-1],temp_array,ad,au,jq,irow));
 499     }
 500     if((ido==-1)||(ido==1)){
 501 
 502       /* solve the linear equation system  */
 503 
 504       if(ido==-1){
 505         if(*isolver==0){
 506 #ifdef SPOOLES
 507           spooles_solve(temp_array,&neq[0]);
 508 #endif
 509         }
 510         else if(*isolver==4){
 511 #ifdef SGI
 512           token=2;
 513           sgi_solve(temp_array,token);
 514 #endif
 515         }
 516         else if(*isolver==5){
 517 #ifdef TAUCS
 518           tau_solve(temp_array,&neq[0]);
 519 #endif
 520         }
 521         else if(*isolver==7){
 522 #ifdef PARDISO
 523           pardiso_solve(temp_array,&neq[0],&symmetryflag);
 524 #endif
 525         }
 526         for(jrow=0;jrow<neq[0];jrow++){
 527           workd[ipntr[1]-1+jrow]=temp_array[jrow];
 528         }
 529       }
 530       else if(ido==1){
 531         if(*isolver==0){
 532 #ifdef SPOOLES
 533           spooles_solve(&workd[ipntr[2]-1],&neq[0]);
 534 #endif
 535         }
 536         else if(*isolver==4){
 537 #ifdef SGI
 538           token=2;
 539           sgi_solve(&workd[ipntr[2]-1],token);
 540 #endif
 541         }
 542         else if(*isolver==5){
 543 #ifdef TAUCS
 544           tau_solve(&workd[ipntr[2]-1],&neq[0]);
 545 #endif
 546         }
 547         else if(*isolver==7){
 548 #ifdef PARDISO
 549           pardiso_solve(&workd[ipntr[2]-1],&neq[0],
 550                &symmetryflag);
 551 #endif
 552         }
 553         for(jrow=0;jrow<neq[0];jrow++){
 554           workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
 555         }
 556       }
 557 
 558     }
 559 
 560     if(ido==2){
 561       FORTRAN(op,(&neq[0],&workd[ipntr[0]-1],&workd[ipntr[1]-1],ad,au,jq,irow));
 562     }
 563 
 564     FORTRAN(dsaupd,(&ido,bmat,&neq[0],which,&nev,&tol,resid,&ncv,z,&dz,iparam,ipntr,workd,
 565             workl,&lworkl,&info));
 566   }
 567 
 568 /*--------------------------------------------------------------------*/
 569 /*
 570    -----------
 571    free memory
 572    -----------
 573 */
 574   if(*isolver==0){
 575 #ifdef SPOOLES
 576     spooles_cleanup();
 577 #endif
 578   }
 579   else if(*isolver==4){
 580 #ifdef SGI
 581     token=2;
 582     sgi_cleanup(token);
 583 #endif
 584   }
 585   else if(*isolver==5){
 586 #ifdef TAUCS
 587     tau_cleanup();
 588 #endif
 589   }
 590   else if(*isolver==7){
 591 #ifdef PARDISO
 592     pardiso_cleanup(&neq[0],&symmetryflag);
 593 #endif
 594   }
 595 
 596   if(info!=0){
 597     printf("*ERROR in arpackbu: info=%" ITGFORMAT "\n",info);
 598     printf("       # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
 599   }         
 600 
 601   NNEW(select,ITG,ncv);
 602   NNEW(d,double,nev);
 603 
 604   FORTRAN(dseupd,(&rvec,howmny,select,d,z,&dz,&sigma,bmat,&neq[0],which,&nev,&tol,resid,
 605           &ncv,z,&dz,iparam,ipntr,workd,workl,&lworkl,&info));
 606 
 607   printf("sigma=%f,d[0]=%f\n\n",sigma,d[0]);
 608   if((5.>d[0]/sigma)||(50000.<d[0]/sigma)){
 609     if(iconverged<-4) {
 610       printf("no convergence for the buckling factor; maybe no buckling occurs");
 611       FORTRAN(stop,());
 612     }
 613     sigma=d[0]/500.;
 614     printf("no convergence; new iteration\n\n");
 615     --iconverged;
 616     SFREE(z);SFREE(d);
 617   }
 618   else{iconverged=0;}
 619      
 620   SFREE(resid);SFREE(workd);SFREE(workl);SFREE(select);SFREE(temp_array);
 621 
 622   } while(iconverged<0);
 623 
 624   SFREE(aub);SFREE(adb);SFREE(au);SFREE(ad);SFREE(b);
 625 
 626   FORTRAN(writebv,(d,&nev));
 627 
 628   /* calculating the displacements and the stresses and storing */
 629   /* the results in frd format for each valid eigenmode */
 630 
 631   NNEW(v,double,mt**nk);
 632   NNEW(fn,double,mt**nk);
 633   NNEW(stn,double,6**nk);
 634   NNEW(inum,ITG,*nk);
 635   NNEW(stx,double,6*mi[0]**ne);
 636   
 637   if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 638   if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 639   if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 640 
 641   lfin=0;
 642   for(j=0;j<nev;++j){
 643 
 644     for(k=0;k<6*mi[0]**ne;k++){eme[k]=0.;}
 645 
 646     lint=lfin;
 647     lfin=lfin+neq[0];
 648 
 649     if(*nprint>0) FORTRAN(writehe,(&j));
 650 
 651 //    memset(&v[0],0.,sizeof(double)*mt**nk);
 652     DMEMSET(v,0,mt**nk,0.);
 653     if(*iperturb==0){
 654       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 655               stx,elcon,
 656               nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
 657               norien,orab,ntmat_,t0,t0,ithermal,
 658               prestr,iprestr,filab,eme,emn,een,iperturb,
 659               f,fn,nactdof,&iout,qa,vold,&z[lint],
 660               nodeboun,ndirboun,xboun,nboun,ipompc,
 661               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,&bet,
 662               &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 663               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 664               ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 665               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 666               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 667               nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
 668               &ne0,xforc,nforc,thicke,shcon,nshcon,
 669               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 670               &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 671               islavsurf,ielprop,prop);}
 672     else{
 673       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 674               stx,elcon,
 675               nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
 676               norien,orab,ntmat_,t0,t1old,ithermal,
 677               prestr,iprestr,filab,eme,emn,een,iperturb,
 678               f,fn,nactdof,&iout,qa,vold,&z[lint],
 679               nodeboun,ndirboun,xboun,nboun,ipompc,
 680               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[0],veold,accold,&bet,
 681               &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 682               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 683               ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 684               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 685               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 686               nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
 687               &ne0,xforc,nforc,thicke,shcon,nshcon,
 688               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 689               &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
 690               islavsurf,ielprop,prop);
 691     }
 692 
 693     ++*kode;
 694     if(strcmp1(&filab[1044],"ZZS")==0){
 695         NNEW(neigh,ITG,40**ne);
 696         NNEW(ipneigh,ITG,*nk);
 697     }
 698 
 699     frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 700             kode,filab,een,t1,fn,&d[j],epn,ielmat,matname,enern,xstaten,
 701             nstate_,&istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
 702             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 703             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 704             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 705             thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 706 
 707     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 708   }
 709 
 710   SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);SFREE(z);SFREE(d);SFREE(eei);
 711   SFREE(xstiff);
 712   if(*nener==1){
 713       SFREE(stiini);SFREE(enerini);}
 714 
 715   if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
 716   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
 717   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
 718 
 719   return;
 720 }
 721 
 722 #endif

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