root/src/arpackcs.c

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

DEFINITIONS

This source file includes following definitions.
  1. arpackcs

   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 <string.h>
  24 #include "CalculiX.h"
  25 #ifdef SPOOLES
  26    #include "spooles.h"
  27 #endif
  28 #ifdef SGI
  29    #include "sgi.h"
  30 #endif
  31 #ifdef TAUCS
  32    #include "tau.h"
  33 #endif
  34 #ifdef MATRIXSTORAGE
  35    #include "matrixstorage.h"
  36 #endif
  37 #ifdef PARDISO
  38    #include "pardiso.h"
  39 #endif
  40 
  41 void arpackcs(double *co, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp,
  42              ITG *ne, 
  43              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  44              ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
  45              ITG *nmpc, 
  46              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  47              ITG *nelemload, char *sideload, double *xload,
  48              ITG *nload, ITG *nactdof, 
  49              ITG *icol, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, 
  50              ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, 
  51              ITG *ilboun,
  52              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  53              double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
  54              ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
  55              double *t0, double *t1, double *t1old, 
  56              ITG *ithermal,double *prestr, ITG *iprestr, 
  57              double *vold,ITG *iperturb, double *sti, ITG *nzs,  
  58              ITG *kode, ITG *mei, double *fei,
  59              char *filab, double *eme,
  60              ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon,
  61              ITG *nplkcon,
  62              double **xstatep, ITG *npmat_, char *matname, ITG *mi,
  63              ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_,
  64              ITG *nstate_, ITG *mcs, ITG *nkon, double **enerp,
  65              char *jobnamec, char *output, char *set, ITG *nset, 
  66              ITG *istartset,
  67              ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  68              char *prset, ITG *nener, ITG *isolver, double *trab, 
  69              ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
  70              char *cbody, ITG *ibody, double *xbody, ITG *nbody,
  71              ITG *nevtot, double *thicke, ITG *nslavs, double *tietol, 
  72              ITG *mpcinfo,ITG *ntie,ITG *istep,
  73              char *tieset,ITG *nintpoint,ITG *mortar,ITG *ifacecount,
  74              ITG **islavsurfp,double **pslavsurfp,double **clearinip,
  75              ITG *nmat,char *typeboun,ITG *ielprop,double *prop){
  76 
  77   /* calls the Arnoldi Package (ARPACK) for cyclic symmetry calculations */
  78   
  79   char bmat[2]="G", which[3]="LM", howmny[2]="A",*lakont=NULL,
  80       description[13]="            ",fneig[132]="",
  81       lakonl[2]=" \0",*lakon=NULL;
  82 
  83   ITG *inum=NULL,k,ido,ldz,iparam[11],ipntr[14],lworkl,idir,nherm=1,
  84     info,rvec=1,*select=NULL,lfin,j,lint,iout=1,nm,index,inode,id,i,idof,
  85     ielas=0,icmd=0,kk,l,nkt,icntrl,*kont=NULL,*ipkont=NULL,*inumt=NULL,
  86     *ielmatt=NULL,net,imag,icomplex,kkv,kk6,iinc=1,nev,ncv,
  87     mxiter,lprev,ilength,ij,i1,i2,iel,ielset,node,indexe,nope,ml1,
  88     *inocs=NULL,*ielcs=NULL,jj,l1,l2,ngraph,is,jrow,*ipobody=NULL,
  89     *inotrt=NULL,symmetryflag=0,inputformat=0,inewton=0,ifreebody,
  90     mass=1, stiffness=1, buckling=0, rhsi=0, intscheme=0,*ncocon=NULL,
  91     coriolis=0,iworsttime,l3,iray,mt,kkx,im,ne0,*integerglob=NULL,
  92     *nshcon=NULL,one=1,irenewxstate,ncont=0,*itietri=NULL,neq2,
  93     *koncont=NULL,ismallsliding=0,*itiefac=NULL,*islavsurf=NULL,
  94     *islavnode=NULL,*imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,
  95     *imastop=NULL,*iponoels=NULL,*inoels=NULL,*ipe=NULL,*ime=NULL,
  96     mpcfree,memmpc_,icascade,maxlenmpc,nkon0,iit=-1,*irow=NULL,nasym=0,
  97     kmax1,kmax2,icfd=0,*inomat=NULL,*ipkon=NULL,*kon=NULL,*ielmat=NULL,
  98     *ielorien=NULL,*islavact=NULL,*islavsurfold=NULL,nslavs_prev_step,
  99     maxprevcontel,iex,iflagact=0,*nmc=NULL;
 100 
 101   double *stn=NULL,*v=NULL,*resid=NULL,*z=NULL,*workd=NULL,*vr=NULL,
 102       *workl=NULL,*d=NULL,sigma,*temp_array=NULL,*vini=NULL,
 103     *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[3],*fext=NULL,*emn=NULL,
 104     *epn=NULL,*stiini=NULL,*fnr=NULL,*fni=NULL,fnreal,fnimag,*emeini=NULL,
 105     *xstateini=NULL,theta=0,pi,*coefmpcnew=NULL,*xstiff=NULL,*vi=NULL,
 106     *vt=NULL,*fnt=NULL,*stnt=NULL,*eent=NULL,*cot=NULL,t[3],ctl,stl,
 107     *t1t=NULL,freq,*stx=NULL,*enern=NULL,*enernt=NULL,*xstaten=NULL,
 108     *eei=NULL,*enerini=NULL,*cocon=NULL,*qfx=NULL,*qfn=NULL,*qfnt=NULL,
 109     tol,fmin,fmax,xreal,ximag,*cgr=NULL,*xloadold=NULL,reltime=1.,constant,
 110     vreal,vimag,*stnr=NULL,*stni=NULL,stnreal,stnimag,*vmax=NULL,
 111     *stnmax=NULL,vl[4],stnl[6],dd,v1,v2,v3,bb,cc,al[3],cm,cn,tt,
 112     worstpsmax,vray[3],worstumax,p1[3],p2[3],q[3],tan[3],*springarea=NULL,
 113     *stxt=NULL,*eenmax=NULL,eenl[6],*emnt=NULL,*clearini=NULL,
 114     *doubleglob=NULL,*shcon=NULL,*cg=NULL,*straight=NULL,*cdn=NULL,
 115     *xmastnor=NULL,*areaslav=NULL,*dc=NULL,*di=NULL,*xnoels=NULL,
 116     *workev=NULL,*temp_array2=NULL,*ener=NULL,*xstate=NULL,cdnimag,
 117     sigmai=0,amp,ampmax,*zstorage=NULL,*au=NULL,*ad=NULL,cdnreal,
 118     *b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
 119     *cdnt=NULL,*cdnr=NULL,*cdni=NULL,
 120     *pslavsurfold=NULL;
 121 
 122   FILE *f1;
 123 
 124   /* dummy arguments for the results call */
 125 
 126   double *veold=NULL,*accold=NULL,bet,gam,dtime,time=1.;
 127 
 128   ITG *ipneigh=NULL,*neigh=NULL;
 129   
 130 #ifdef SGI
 131   ITG token;
 132 #endif
 133   
 134   irow=*irowp;ener=*enerp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
 135   kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;
 136 
 137   islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
 138   
 139   if(*mortar==0){
 140       maxprevcontel=*nslavs;
 141   }else if(*mortar==1){
 142       maxprevcontel=*nintpoint;
 143   }
 144   nslavs_prev_step=*nslavs;
 145   
 146   mt=mi[1]+1;
 147   pi=4.*atan(1.);
 148   constant=180./pi;
 149   
 150   /* copying the frequency parameters */
 151   
 152   nev=mei[0];
 153   ncv=mei[1];
 154   mxiter=mei[2];
 155   tol=fei[0];
 156   fmin=2*pi*fei[1];
 157   fmax=2*pi*fei[2];
 158   
 159   /* assigning the body forces to the elements */ 
 160   
 161   if(*nbody>0){
 162       ifreebody=*ne+1;
 163       NNEW(ipobody,ITG,2*ifreebody**nbody);
 164       for(k=1;k<=*nbody;k++){
 165           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 166                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 167           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 168       }
 169       RENEW(ipobody,ITG,2*(ifreebody-1));
 170       if(inewton==1){
 171           printf("*ERROR in arpackcs: generalized gravity loading is not allowed in frequency calculations");
 172           FORTRAN(stop,());
 173       }
 174   }
 175   
 176   ne0=*ne;nkon0=*nkon;
 177   
 178   /* contact conditions */
 179   
 180   if(*iperturb!=0){
 181       
 182       memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
 183       maxlenmpc=mpcinfo[3];
 184       
 185       if(*nslavs==0){irenewxstate=1;}else{irenewxstate=0;}
 186       inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
 187               lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
 188               &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
 189               mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
 190               nmpc,&mpcfree,&memmpc_,
 191               &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 192               iperturb,ikboun,nboun,co,istep,&xnoels);
 193       
 194       if(ncont!=0){
 195 
 196           NNEW(cg,double,3*ncont);
 197           NNEW(straight,double,16*ncont);
 198           
 199           /* 11 instead of 10: last position is reserved for the
 200              local contact spring element number; needed as
 201              pointer into springarea */
 202           
 203           if(*mortar==0){
 204               RENEW(kon,ITG,*nkon+11**nslavs);
 205               NNEW(springarea,double,2**nslavs);
 206               if(*nener==1){
 207                   RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
 208               }
 209               RENEW(ipkon,ITG,*ne+*nslavs);
 210               RENEW(lakon,char,8*(*ne+*nslavs));
 211               
 212               if(*norien>0){
 213                   RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
 214                   for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
 215               }
 216               
 217               RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
 218               for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
 219               
 220               if((maxprevcontel==0)&&(*nslavs!=0)){
 221                   RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
 222                   for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
 223                       xstate[k]=0.;
 224                   }
 225               }
 226               maxprevcontel=*nslavs;
 227               
 228               NNEW(areaslav,double,*ifacecount);
 229               NNEW(xmastnor,double,3*nmastnode[*ntie]);
 230           }else if(*mortar==1){
 231               NNEW(islavact,ITG,nslavnode[*ntie]);
 232               DMEMSET(islavact,0,nslavnode[*ntie],1);
 233               if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
 234               NNEW(xmastnor,double,3*nmastnode[*ntie]);
 235 
 236               if(*nstate_!=0){
 237                   if(maxprevcontel!=0){
 238                       NNEW(islavsurfold,ITG,2**ifacecount+2);
 239                       NNEW(pslavsurfold,double,3**nintpoint);
 240                       memcpy(&islavsurfold[0],&islavsurf[0],
 241                              sizeof(ITG)*(2**ifacecount+2));
 242                       memcpy(&pslavsurfold[0],&pslavsurf[0],
 243                              sizeof(double)*(3**nintpoint));
 244                   }
 245               }
 246 
 247               *nintpoint=0;
 248               
 249               precontact(&ncont,ntie,tieset,nset,set,istartset,
 250                          iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
 251                          cg,straight,co,vold,istep,&iinc,&iit,itiefac,
 252                          islavsurf,islavnode,imastnode,nslavnode,nmastnode,
 253                          imastop,mi,ipe,ime,tietol,&iflagact,
 254                          nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
 255                          nslavs);
 256               
 257               /* changing the dimension of element-related fields */
 258               
 259               RENEW(kon,ITG,*nkon+22**nintpoint);
 260               RENEW(springarea,double,2**nintpoint);
 261               RENEW(pmastsurf,double,6**nintpoint);
 262               
 263               if(*nener==1){
 264                   RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
 265               }
 266               RENEW(ipkon,ITG,*ne+*nintpoint);
 267               RENEW(lakon,char,8*(*ne+*nintpoint));
 268               
 269               if(*norien>0){
 270                   RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
 271                   for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
 272               }
 273               RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
 274               for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
 275           }
 276           
 277           /* generating contact spring elements */
 278           
 279           contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
 280                   ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
 281                   co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
 282                   &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
 283                   &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 284                   iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
 285                   itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
 286                   imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
 287                   xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta);
 288           
 289           printf("number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
 290                   
 291           /* interpolating the state variables */
 292 
 293           if(*mortar==1){
 294               if(*nstate_!=0){
 295                   if(maxprevcontel!=0){
 296                       RENEW(xstateini,double,
 297                             *nstate_*mi[0]*(ne0+maxprevcontel));
 298                       for(k=*nstate_*mi[0]*ne0;
 299                           k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
 300                           xstateini[k]=xstate[k];
 301                       }
 302                   }
 303                   
 304                   RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
 305                   for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
 306                       xstate[k]=0.;
 307                   }
 308                   
 309                   if((*nintpoint>0)&&(maxprevcontel>0)){
 310                       iex=2;
 311                       
 312                       /* interpolation of xstate */
 313                       
 314                       FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
 315                                &ne0,mi,xstate,pslavsurf,nstate_,
 316                                xstateini,islavsurf,islavsurfold,
 317                                pslavsurfold));
 318                       
 319                   }
 320                   
 321                   if(maxprevcontel!=0){
 322                       SFREE(islavsurfold);SFREE(pslavsurfold);
 323                   }
 324 
 325                   maxprevcontel=*nintpoint;
 326                   
 327                   RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
 328                   for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
 329                       xstateini[k]=xstate[k];
 330                   }
 331               }
 332           }
 333           
 334           /* determining the structure of the stiffness/mass matrix */
 335           
 336           remastructar(ipompc,&coefmpc,&nodempc,nmpc,
 337                        &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
 338                        labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
 339                        kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
 340                        neq,nzs,nmethod,ithermal,iperturb,&mass,mi,ics,cs,
 341                        mcs,mortar,typeboun);
 342       }
 343   }
 344   
 345   /* field for initial values of state variables (needed if
 346      previous static step was viscoplastic and for contact */
 347   
 348   if((*nstate_!=0)&&((*mortar==0)||(ncont==0))){
 349       NNEW(xstateini,double,*nstate_*mi[0]*(ne0+*nslavs));
 350       for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
 351           xstateini[k]=xstate[k];
 352       }
 353   }
 354   
 355   /* determining the internal forces and the stiffness coefficients */
 356   
 357   NNEW(f,double,neq[1]);
 358   
 359   /* allocating a field for the stiffness matrix */
 360   
 361   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 362   
 363   iout=-1;
 364   NNEW(v,double,mt**nk);
 365   memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
 366   NNEW(fn,double,mt**nk);
 367   NNEW(stx,double,6*mi[0]**ne);
 368   NNEW(inum,ITG,*nk);
 369   if(*iperturb==0){
 370       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 371               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 372               ielorien,norien,orab,ntmat_,t0,t0,ithermal,
 373               prestr,iprestr,filab,eme,emn,een,iperturb,
 374               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 375               ndirboun,xboun,nboun,ipompc,
 376               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 377               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 378               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 379               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 380               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 381               iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 382               fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 383               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 384               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 385               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 386               islavsurf,ielprop,prop);
 387   }else{
 388       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 389               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 390               ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
 391               prestr,iprestr,filab,eme,emn,een,iperturb,
 392               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 393               ndirboun,xboun,nboun,ipompc,
 394               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 395               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 396               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 397               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 398               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 399               iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 400               fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 401               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 402               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 403               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 404               islavsurf,ielprop,prop);
 405   }
 406   SFREE(f);SFREE(v);SFREE(fn);SFREE(stx);SFREE(inum);
 407   iout=1;
 408   
 409   /* for the frequency analysis linear strain and elastic properties
 410      are used */
 411   
 412   iperturb[1]=0;ielas=1;
 413   
 414   /* determining the maximum number of sectors to be plotted */
 415   
 416   ngraph=1;
 417   for(j=0;j<*mcs;j++){
 418       if(cs[17*j+4]>ngraph) ngraph=cs[17*j+4];
 419   }
 420   
 421   /* assigning nodes and elements to sectors */
 422   
 423   NNEW(inocs,ITG,*nk);
 424   NNEW(ielcs,ITG,*ne);
 425   ielset=cs[12];
 426   if((*mcs!=1)||(ielset!=0)){
 427       for(i=0;i<*nk;i++) inocs[i]=-1;
 428       for(i=0;i<*ne;i++) ielcs[i]=-1;
 429   }
 430   
 431   for(i=0;i<*mcs;i++){
 432       is=cs[17*i+4];
 433       if((is==1)&&(*mcs==1)) continue;
 434       ielset=cs[17*i+12];
 435       if(ielset==0) continue;
 436       for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
 437           if(ialset[i1]>0){
 438               iel=ialset[i1]-1;
 439               if(ipkon[iel]<0) continue;
 440               ielcs[iel]=i;
 441               indexe=ipkon[iel];
 442               if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 443               else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 444               else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 445               else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 446               else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 447               else if (strcmp1(&lakon[8*iel+3],"6")==0)nope=6;
 448               else if (strcmp1(&lakon[8*iel],"ES")==0){
 449                   lakonl[0]=lakon[8*iel+7];
 450                   nope=atoi(lakonl)+1;}
 451               else continue;
 452               
 453               for(i2=0;i2<nope;++i2){
 454                   node=kon[indexe+i2]-1;
 455                   inocs[node]=i;
 456               }
 457           }
 458           else{
 459               iel=ialset[i1-2]-1;
 460               do{
 461                   iel=iel-ialset[i1];
 462                   if(iel>=ialset[i1-1]-1) break;
 463                   if(ipkon[iel]<0) continue;
 464                   ielcs[iel]=i;
 465                   indexe=ipkon[iel];
 466                   if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 467                   else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 468                   else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 469                   else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 470                   else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 471                   else {nope=6;}
 472                   for(i2=0;i2<nope;++i2){
 473                       node=kon[indexe+i2]-1;
 474                       inocs[node]=i;
 475                   }
 476               }while(1);
 477           }
 478       } 
 479   }
 480   
 481   /* loop over the nodal diameters */
 482   
 483   for(nm=cs[1];nm<=cs[2];++nm){
 484       
 485       NNEW(ad,double,neq[1]);
 486       NNEW(au,double,nzs[1]);
 487       
 488       NNEW(adb,double,neq[1]);
 489       NNEW(aub,double,nzs[1]);
 490       
 491       NNEW(fext,double,neq[1]);
 492       if(*iperturb==0){
 493           FORTRAN(mafillsmcs,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 494                 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 495                 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 496                 ad,au,fext,nactdof,icol,jq,irow,&neq[1],nzl,nmethod,
 497                 ikmpc,ilmpc,ikboun,ilboun,
 498                 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 499                 ielorien,norien,orab,ntmat_,
 500                 t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 501                 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 502                 xstiff,npmat_,&dtime,matname,mi,
 503                 ics,cs,&nm,ncmat_,labmpc,&mass,&stiffness,&buckling,&rhsi,
 504                 &intscheme,mcs,&coriolis,ibody,xloadold,&reltime,ielcs,veold,
 505                 springarea,thicke,integerglob,doubleglob,
 506                 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 507                 pmastsurf,mortar,clearini,ielprop,prop));
 508       }
 509       else{
 510           FORTRAN(mafillsmcs,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 511                 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 512                 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 513                 ad,au,fext,nactdof,icol,jq,irow,&neq[1],nzl,nmethod,
 514                 ikmpc,ilmpc,ikboun,ilboun,
 515                 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 516                 ielorien,norien,orab,ntmat_,
 517                 t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
 518                 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 519                 xstiff,npmat_,&dtime,matname,mi,
 520                 ics,cs,&nm,ncmat_,labmpc,&mass,&stiffness,&buckling,&rhsi,
 521                 &intscheme,mcs,&coriolis,ibody,xloadold,&reltime,ielcs,veold,
 522                 springarea,thicke,integerglob,doubleglob,
 523                 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 524                 pmastsurf,mortar,clearini,ielprop,prop));
 525           
 526           if(nasym==1){
 527               RENEW(au,double,nzs[2]+nzs[1]);
 528               RENEW(aub,double,nzs[2]+nzs[1]);
 529               symmetryflag=2;
 530               inputformat=1;
 531 
 532               FORTRAN(mafillsmcsas,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
 533                 xboun,nboun,
 534                 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 535                 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 536                 ad,au,fext,nactdof,icol,jq,irow,&neq[1],nzl,nmethod,
 537                 ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
 538                 nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
 539                 t0,t1,ithermal,prestr,
 540                 iprestr,vold,iperturb,sti,nzs,stx,adb,aub,iexpl,plicon,
 541                 nplicon,plkcon,nplkcon,xstiff,npmat_,&dtime,
 542                 matname,mi,ics,cs,&nm,ncmat_,labmpc,&mass,&stiffness,&buckling,
 543                 &rhsi,&intscheme,mcs,&coriolis,ibody,xloadold,&reltime,ielcs,
 544                 veold,springarea,thicke,integerglob,doubleglob,
 545                 tieset,istartset,iendset,ialset,ntie,&nasym,nstate_,xstateini,
 546                 xstate,pslavsurf,pmastsurf,mortar,clearini,ielprop,prop));
 547               
 548           }
 549       }
 550       
 551       SFREE(fext);
 552       
 553       if(*nmethod==0){
 554           
 555           /* error occurred in mafill: storing the geometry in frd format */
 556           
 557           ++*kode;
 558           NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 559           if(strcmp1(&filab[1044],"ZZS")==0){
 560               NNEW(neigh,ITG,40**ne);
 561               NNEW(ipneigh,ITG,*nk);
 562           }
 563           
 564           frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 565               kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
 566               nstate_,istep,&iinc,ithermal,qfn,&j,&nm,trab,inotr,
 567               ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 568               mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 569               cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 570               thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
 571           
 572           if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 573           SFREE(inum);FORTRAN(stop,());
 574           
 575       }
 576       
 577       /* LU decomposition of the left hand matrix */
 578       
 579       if(nasym==1){sigma=0.;}else{sigma=1.;}
 580       
 581       if(*isolver==0){
 582 #ifdef SPOOLES
 583           spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,&symmetryflag,
 584                          &inputformat,&nzs[2]);
 585 #else
 586           printf("*ERROR in arpackcs: the SPOOLES library is not linked\n\n");
 587           FORTRAN(stop,());
 588 #endif
 589       }
 590       else if(*isolver==4){
 591 #ifdef SGI
 592           token=1;
 593           sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,token);
 594 #else
 595           printf("*ERROR in arpackcs: the SGI library is not linked\n\n");
 596           FORTRAN(stop,());
 597 #endif
 598       }
 599       else if(*isolver==5){
 600 #ifdef TAUCS
 601           tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],nzs);
 602 #else
 603           printf("*ERROR in arpackcs: the TAUCS library is not linked\n\n");
 604           FORTRAN(stop,());
 605 #endif
 606       }
 607       else if(*isolver==6){
 608 #ifdef MATRIXSTORAGE
 609           matrixstorage(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1],
 610                         ntrans,inotr,trab,co,nk,nactdof,jobnamec,mi,ipkon,
 611                         lakon,kon,ne,mei,nboun,nmpc,cs,mcs);
 612 #else
 613           printf("*ERROR in arpack: the MATRIXSTORAGE library is not linked\n\n");
 614           FORTRAN(stop,());
 615 #endif
 616       }
 617       if(*isolver==7){
 618 #ifdef PARDISO
 619           pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,
 620                          &symmetryflag,&inputformat,jq,&nzs[2]);
 621 #else
 622           printf("*ERROR in arpackcs: the PARDISO library is not linked\n\n");
 623           FORTRAN(stop,());
 624 #endif
 625       }
 626       
 627       SFREE(au);SFREE(ad);
 628       
 629       /* calculating the eigenvalues and eigenmodes */
 630       
 631       printf(" Calculating the eigenvalues and the eigenmodes\n");
 632       
 633       ido=0;
 634       ldz=neq[1];
 635       for(k=0;k<11;k++)iparam[k]=0;
 636       iparam[0]=1;
 637       iparam[2]=mxiter;
 638       iparam[3]=1;
 639       iparam[6]=3;
 640       
 641       info=0;
 642       
 643       NNEW(resid,double,neq[1]);
 644       NNEW(z,double,(long long)ncv*neq[1]);
 645       NNEW(workd,double,3*neq[1]);
 646       
 647       if(nasym==1){
 648           lworkl=3*ncv*(2+ncv);
 649           NNEW(workl,double,lworkl);
 650           FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
 651                           workl,&lworkl,&info));
 652           /*      neq2=neq[1]/2;
 653                   lworkl=6*ncv*ncv+10*ncv;
 654                   NNEW(workl,double,lworkl);
 655                   NNEW(rwork,double,ncv);
 656                   FORTRAN(znaupd,(&ido,bmat,&neq2,which,&nev,&tol,resid,&ncv,z,&neq2,
 657                   iparam,ipntr,workd,workl,&lworkl,rwork,&info));
 658                   NNEW(temp_array2,double,neq[1]);*/
 659       }else{
 660           lworkl=ncv*(8+ncv);
 661           NNEW(workl,double,lworkl);
 662           FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
 663                           iparam,ipntr,workd,workl,&lworkl,&info));
 664       }
 665       
 666       NNEW(temp_array,double,neq[1]);
 667       
 668       while((ido==-1)||(ido==1)||(ido==2)){
 669           if(ido==-1){
 670               if(nasym==1){
 671                   FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow,nzs));
 672                   /*        for(jrow=0;jrow<neq2;jrow++){
 673                             temp_array2[jrow]=workd[ipntr[0]-1+2*jrow];
 674                             temp_array2[neq2+jrow]=workd[ipntr[0]+2*jrow];
 675                             }
 676                             FORTRAN(opas,(&neq[1],aux,temp_array2,temp_array,adb,aub,
 677                             icol,irow,nzl,nzs));*/
 678               }else{
 679                   FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,
 680                               jq,irow));
 681               }
 682           }
 683           
 684           /* solve the linear equation system  */
 685           
 686           if((ido==-1)||(ido==1)){
 687               
 688               if(ido==-1){
 689                   if(*isolver==0){
 690 #ifdef SPOOLES
 691                       spooles_solve(temp_array,&neq[1]);
 692 #endif
 693                   }
 694                   else if(*isolver==4){
 695 #ifdef SGI
 696                       sgi_solve(temp_array,token);
 697 #endif
 698                   }
 699                   else if(*isolver==5){
 700 #ifdef TAUCS
 701                       tau_solve(temp_array,&neq[1]);
 702 #endif
 703                   }
 704                   else if(*isolver==7){
 705 #ifdef PARDISO
 706                       pardiso_solve(temp_array,&neq[1],&symmetryflag);
 707 #endif
 708                   }
 709                   /*        if(nasym==1){
 710                             for(jrow=0;jrow<neq2;jrow++){
 711                             workd[ipntr[1]-1+2*jrow]=temp_array[jrow];
 712                             workd[ipntr[1]+2*jrow]=temp_array[neq2+jrow];
 713                             }
 714                             }else{*/
 715                   for(jrow=0;jrow<neq[1];jrow++){
 716                       workd[ipntr[1]-1+jrow]=temp_array[jrow];
 717                       //                }
 718                   }
 719               }
 720               else if(ido==1){
 721                   /*        if(nasym==1){
 722                             for(jrow=0;jrow<neq2;jrow++){
 723                             temp_array2[jrow]=workd[ipntr[2]-1+2*jrow];
 724                             temp_array2[neq2+jrow]=workd[ipntr[2]+2*jrow];
 725                             }
 726                             if(*isolver==0){
 727                             #ifdef SPOOLES
 728                             spooles_solve(temp_array2,&neq[1]);
 729                             #endif
 730                             }
 731                             else if(*isolver==4){
 732                             #ifdef SGI
 733                             sgi_solve(temp_array2,token);
 734                             #endif
 735                             }
 736                             else if(*isolver==5){
 737                             #ifdef TAUCS
 738                             tau_solve(temp_array2,&neq[1]);
 739                             #endif
 740                             }
 741                             else if(*isolver==7){
 742                             #ifdef PARDISO
 743                             pardiso_solve(temp_array2,&neq[1],&symmetryflag);
 744                             #endif
 745                             }
 746                             for(jrow=0;jrow<neq2;jrow++){
 747                             workd[ipntr[1]-1+2*jrow]=temp_array2[jrow];
 748                             workd[ipntr[1]+2*jrow]=temp_array2[neq2+jrow];
 749                             }
 750                             }else{*/
 751                   if(*isolver==0){
 752 #ifdef SPOOLES
 753                       spooles_solve(&workd[ipntr[2]-1],&neq[1]);
 754 #endif
 755                   }
 756                   else if(*isolver==4){
 757 #ifdef SGI
 758                       sgi_solve(&workd[ipntr[2]-1],token);
 759 #endif
 760                   }
 761                   else if(*isolver==5){
 762 #ifdef TAUCS
 763                       tau_solve(&workd[ipntr[2]-1],&neq[1]);
 764 #endif
 765                   }
 766                   else if(*isolver==7){
 767 #ifdef PARDISO
 768                       pardiso_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag);
 769 #endif
 770                   }
 771                   for(jrow=0;jrow<neq[1];jrow++){
 772                       workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
 773                   }
 774                   //        }
 775                   
 776               }
 777           }
 778           
 779           if(ido==2){
 780               if(nasym==1){
 781                   FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
 782                                 adb,aub,jq,irow,nzs));
 783                   /*        for(jrow=0;jrow<neq2;jrow++){
 784                             temp_array2[jrow]=workd[ipntr[0]-1+2*jrow];
 785                             temp_array2[neq2+jrow]=workd[ipntr[0]+2*jrow];
 786                             }
 787                             FORTRAN(opas,(&neq[1],aux,temp_array2,temp_array,
 788                             adb,aub,icol,irow,nzl,nzs));
 789                             for(jrow=0;jrow<neq2;jrow++){
 790                             workd[ipntr[1]-1+2*jrow]=temp_array[jrow];
 791                             workd[ipntr[1]+2*jrow]=temp_array[neq2+jrow];
 792                             }*/
 793               }else{
 794                   FORTRAN(op,(neq,&workd[ipntr[0]-1],&workd[ipntr[1]-1],
 795                               adb,aub,jq,irow));
 796               }
 797           }
 798           
 799           if(nasym==1){
 800               FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
 801                               iparam,ipntr,workd,workl,&lworkl,&info));
 802               /*        FORTRAN(znaupd,(&ido,bmat,&neq2,which,&nev,&tol,resid,&ncv,z,&neq2,
 803                         iparam,ipntr,workd,workl,&lworkl,rwork,&info));*/
 804           }else{
 805               FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,
 806                               z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
 807           }
 808       }
 809       
 810 /*--------------------------------------------------------------------*/
 811 /*
 812   -----------
 813   free memory
 814   -----------
 815 */
 816       SFREE(temp_array);SFREE(temp_array2);
 817       if(*isolver==0){
 818 #ifdef SPOOLES
 819           spooles_cleanup();
 820 #endif
 821       }
 822       else if(*isolver==4){
 823 #ifdef SGI
 824           sgi_cleanup(token);
 825 #endif
 826       }
 827       else if(*isolver==5){
 828 #ifdef TAUCS
 829           tau_cleanup();
 830 #endif
 831       }
 832       else if(*isolver==7){
 833 #ifdef PARDISO
 834           pardiso_cleanup(&neq[1],&symmetryflag);
 835 #endif
 836       }
 837       
 838       if(info!=0){
 839           printf("*ERROR in arpackcs: info=%" ITGFORMAT "\n",info);
 840           printf("       # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
 841       }         
 842       
 843       NNEW(select,ITG,ncv);
 844       
 845       if(nasym==1){
 846           NNEW(d,double,nev+1);
 847           NNEW(di,double,nev+1);
 848           NNEW(workev,double,3*ncv);
 849           FORTRAN(dneupd,(&rvec,howmny,select,d,di,z,&ldz,&sigma,&sigmai,
 850                           workev,bmat,&neq[1],which,&nev,&tol,resid,
 851                           &ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
 852           SFREE(workev);
 853           NNEW(dc,double,2*nev);
 854           NNEW(nmc,ITG,nev);
 855           
 856           /* storing as complex number and taking the square root */
 857           
 858           for(j=0;j<nev;j++){
 859               dc[2*j]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])+d[j])/sqrt(2.);
 860               dc[2*j+1]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])-d[j])/sqrt(2.);
 861               if(di[j]<0.) dc[2*j+1]=-dc[2*j+1];
 862               nmc[j]=nm;
 863           }
 864 //        FORTRAN(writeevcomplex,(dc,&nev,&fmin,&fmax));
 865           FORTRAN(writeevcscomplex,(dc,&nev,nmc,&fmin,&fmax));
 866           SFREE(di);SFREE(dc);SFREE(nmc);
 867           /*      dNNEW(      d,double,2*nev+2);
 868                   NNEW(workev,double,4*ncv);
 869                   FORTRAN(zneupd,(&rvec,howmny,select,d,z,&neq2,sigmaz,workev,bmat,&neq2,
 870                   which,&nev,&tol,resid,&ncv,
 871                   z,&neq2,iparam,ipntr,workd,workl,&lworkl,rwork,&info));
 872                   SFREE(workev);SFREE(rwork);*/
 873           
 874           /* taking the square root */
 875           
 876           /*      dcNNEW(     dc,double,2*nev);
 877                   NNEW(nmc,ITG,nev);
 878                   for(j=0;j<nev;j++){
 879                   dc[2*j]=sqrt(sqrt(d[2*j]*d[2*j]+d[2*j+1]*d[2*j+1])+d[2*j])/sqrt(2.);
 880                   dc[2*j+1]=sqrt(sqrt(d[2*j]*d[2*j]+d[2*j+1]*d[2*j+1])-d[2*j])/sqrt(2.);
 881                   if(d[2*j+1]<0.) dc[2*j+1]=-dc[2*j+1];
 882                   nmc[j]=nm;
 883                   }
 884                   FORTRAN(writeevcscomplex,(d,&nev,nmc,&fmin,&fmax));
 885                   FORTRAN(writeevcscomplex,(dc,&nev,nmc,&fmin,&fmax));SFREE(dc);SFREE(nmc);*/
 886           
 887           /* rearranging the real part of the eigenvalues */
 888           
 889           /*      for(j=0;j<nev;j++)d[j]=d[2*j];
 890                   RENEW(d,double,nev);*/
 891           
 892           /* rearranging the eigenvectors 
 893              (first all real parts, then all imaginary parts */
 894           
 895           /*      vNNEW(      v,double,neq[1]);
 896                   for(j=0;j<nev;j++){
 897                   for(k=0;k<neq[1];k++){
 898                   v[k]=z[j*neq[1]+k];
 899                   }
 900                   for(k=0;k<neq2;k++){
 901                   z[j*neq[1]+k]=v[2*k];
 902                   z[j*neq[1]+neq2+k]=v[2*k+1];
 903                   }
 904                   }
 905                   SFREE(v);*/
 906       }else{
 907           NNEW(d,double,nev);
 908           FORTRAN(dseupd,(&rvec,howmny,select,d,z,&ldz,&sigma,bmat,&neq[1],which,&nev,
 909                           &tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
 910           FORTRAN(writeevcs,(d,&nev,&nm,&fmin,&fmax));
 911       }
 912       SFREE(select);SFREE(resid);SFREE(workd);SFREE(workl);
 913       
 914       /* for double eigenmodes: the eigenmode for which the largest
 915          amplitude has the lowest dof comes first */
 916       
 917       neq2=neq[1]/2;
 918       for (j=0;j<nev;j+=2){
 919           
 920           ampmax=0.;kmax1=0;
 921           for(k=0;k<neq2;k++){
 922               amp=z[2*j*neq[1]+k]*z[2*j*neq[1]+k]+
 923                   z[2*j*neq[1]+neq2+k]*z[2*j*neq[1]+neq2+k];
 924               if(amp>ampmax){ampmax=amp;kmax1=k;}
 925           }
 926           
 927           ampmax=0.;kmax2=0;
 928           for(k=0;k<neq2;k++){
 929               amp=z[(2*j+1)*neq[1]+k]*z[(2*j+1)*neq[1]+k]+
 930                   z[(2*j+1)*neq[1]+neq2+k]*z[(2*j+1)*neq[1]+neq2+k];
 931               if(amp>ampmax){ampmax=amp;kmax2=k;}
 932           }
 933           
 934           if(kmax2<kmax1){
 935               printf("exchange!\n");
 936               NNEW(zstorage,double,neq[1]);
 937               memcpy(zstorage,&z[2*j*neq[1]],sizeof(double)*neq[1]);
 938               memcpy(&z[2*j*neq[1]],&z[(2*j+1)*neq[1]],sizeof(double)*neq[1]);
 939               memcpy(&z[(2*j+1)*neq[1]],zstorage,sizeof(double)*neq[1]);
 940               SFREE(zstorage);
 941           }
 942       }
 943       
 944       /* writing the eigenvalues and mass matrix to a binary file */
 945       
 946       if(mei[3]==1){
 947           
 948           strcpy(fneig,jobnamec);
 949           strcat(fneig,".eig");
 950           
 951           /* the first time the file is erased before writing, all subsequent
 952              times the data is appended */
 953           
 954           if(*nevtot==0){
 955               if((f1=fopen(fneig,"wb"))==NULL){
 956                   printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
 957                   exit(0);
 958               }
 959               
 960               /* storing a one as indication that this was a
 961                  cyclic symmetry calculation */
 962               
 963               if(fwrite(&one,sizeof(ITG),1,f1)!=1){
 964                   printf("*ERROR saving the cyclic symmetry flag to the eigenvalue file...");
 965                   exit(0);
 966               }
 967               
 968               /* Hermitian */
 969               
 970               if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
 971                   printf("*ERROR saving the Hermitian flag to the eigenvalue file...");
 972                   exit(0);
 973               }
 974               
 975           }else{
 976               if((f1=fopen(fneig,"ab"))==NULL){
 977                   printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
 978                   exit(0);
 979               }
 980           }
 981           
 982           /* nodal diameter */
 983           
 984           if(fwrite(&nm,sizeof(ITG),1,f1)!=1){
 985               printf("*ERROR saving the nodal diameter to the eigenvalue file...");
 986               exit(0);
 987           }
 988           
 989           /* number of eigenfrequencies */
 990           
 991           if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
 992               printf("*ERROR saving the number of eigenfrequencies to the eigenvalue file...");
 993               exit(0);
 994           }
 995           
 996           /* the eigenfrequencies are stored as radians/time */
 997           
 998           if(fwrite(d,sizeof(double),nev,f1)!=nev){
 999               printf("*ERROR saving the eigenfrequencies to the eigenvalue file...");
1000               exit(0);
1001           }
1002           if(*nevtot==0){
1003               
1004               /* mass matrix */
1005               
1006               if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
1007                   printf("*ERROR saving the diagonal of the mass matrix to the eigenvalue file...");
1008                   exit(0);
1009               }
1010               if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
1011                   printf("*ERROR saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
1012                   exit(0);
1013               }
1014           }
1015       }
1016 
1017       /* calculating the participation factors and the relative effective
1018          modal mass */
1019 
1020       if(*ithermal!=2){
1021           FORTRAN(effectivemodalmass,(neq,nactdof,mi,adb,aub,jq,irow,&nev,z,co,nk));
1022       }
1023       
1024       /* calculating the displacements and the stresses and storing */
1025       /* the results in frd format for each valid eigenmode */
1026       
1027       NNEW(v,double,2*mt**nk);
1028       NNEW(fn,double,2*mt**nk);
1029       if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1653],"MAXS")==0)|| 
1030          (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1031          (strcmp1(&filab[1044],"ERR ")==0)) 
1032           NNEW(stn,double,12**nk);
1033       
1034       if((strcmp1(&filab[261],"E   ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) 
1035           NNEW(een,double,12**nk);
1036       if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,2**nk);
1037       if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,12**nk);
1038       if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1039          &&(*mortar==1)) NNEW(cdn,double,12**nk);
1040       
1041       NNEW(inum,ITG,*nk);
1042       NNEW(stx,double,2*6*mi[0]**ne);
1043       
1044       NNEW(temp_array,double,neq[1]);
1045       NNEW(coefmpcnew,double,*mpcend);
1046       
1047       NNEW(cot,double,3**nk*ngraph);
1048       if(*ntrans>0){NNEW(inotrt,ITG,2**nk*ngraph);}
1049       if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0))
1050           
1051 // real and imaginary part of the displacements
1052           
1053           NNEW(vt,double,2*mt**nk*ngraph);
1054       if(strcmp1(&filab[87],"NT  ")==0)
1055           NNEW(t1t,double,*nk*ngraph);
1056       if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1057          (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1058           
1059 // real and imaginary part of the stresses
1060           
1061           NNEW(stnt,double,2*6**nk*ngraph);
1062       if(strcmp1(&filab[261],"E   ")==0) NNEW(eent,double,2*6**nk*ngraph);
1063       if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0))
1064           
1065 // real and imaginary part of the forces
1066           
1067           NNEW(fnt,double,2*mt**nk*ngraph);
1068       if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emnt,double,2*6**nk*ngraph);
1069       if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1070          &&(*mortar==1)) NNEW(cdnt,double,2*6**nk*ngraph);
1071       if(strcmp1(&filab[522],"ENER")==0)
1072           NNEW(enernt,double,*nk*ngraph);
1073       if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
1074          ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0)))
1075           NNEW(stxt,double,2*6*mi[0]**ne*ngraph);
1076       
1077       NNEW(kont,ITG,*nkon*ngraph);
1078       NNEW(ipkont,ITG,*ne*ngraph);
1079       for(l=0;l<*ne*ngraph;l++)ipkont[l]=-1;
1080       NNEW(lakont,char,8**ne*ngraph);
1081       NNEW(inumt,ITG,*nk*ngraph);
1082       NNEW(ielmatt,ITG,mi[2]**ne*ngraph);
1083       
1084       nkt=ngraph**nk;
1085       net=ngraph**ne;
1086       
1087       /* copying the coordinates of the first sector */
1088       
1089       for(l=0;l<3**nk;l++){cot[l]=co[l];}
1090       if(*ntrans>0){for(l=0;l<*nk;l++){inotrt[2*l]=inotr[2*l];}}
1091       for(l=0;l<*nkon;l++){kont[l]=kon[l];}
1092       for(l=0;l<*ne;l++){ipkont[l]=ipkon[l];}
1093       for(l=0;l<8**ne;l++){lakont[l]=lakon[l];}
1094       for(l=0;l<*ne;l++){ielmatt[mi[2]*l]=ielmat[mi[2]*l];}
1095       
1096       /* generating the coordinates for the other sectors */
1097       
1098       icntrl=1;
1099       
1100       FORTRAN(rectcyl,(cot,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1101       
1102       for(jj=0;jj<*mcs;jj++){
1103           is=cs[17*jj+4];
1104           for(i=1;i<is;i++){
1105               
1106               theta=i*2.*pi/cs[17*jj];
1107               
1108               for(l=0;l<*nk;l++){
1109                   if(inocs[l]==jj){
1110                       cot[3*l+i*3**nk]=cot[3*l];
1111                       cot[1+3*l+i*3**nk]=cot[1+3*l]+theta;
1112                       cot[2+3*l+i*3**nk]=cot[2+3*l];
1113                       if(*ntrans>0){inotrt[2*l+i*2**nk]=inotrt[2*l];}
1114                   }
1115               }
1116               for(l=0;l<*nkon;l++){kont[l+i**nkon]=kon[l]+i**nk;}
1117               for(l=0;l<*ne;l++){
1118                   if(ielcs[l]==jj){
1119                       if(ipkon[l]>=0){
1120                           ipkont[l+i**ne]=ipkon[l]+i**nkon;
1121                           ielmatt[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
1122                           for(l1=0;l1<8;l1++){
1123                               l2=8*l+l1;
1124                               lakont[l2+i*8**ne]=lakon[l2];
1125                           }
1126                       }
1127                   }
1128               }
1129           }
1130       }
1131       
1132       icntrl=-1;
1133       
1134       FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1135                        &imag,mi,emn));
1136       
1137       /* check that the tensor fields which are extrapolated from the
1138          integration points are requested in global coordinates */
1139       
1140       if(strcmp1(&filab[174],"S   ")==0){
1141           if((strcmp1(&filab[179],"L")==0)&&(*norien>0)){
1142               printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculations\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1143               strcpy1(&filab[179],"G",1);
1144           }
1145       }
1146       
1147       if(strcmp1(&filab[261],"E   ")==0){
1148           if((strcmp1(&filab[266],"L")==0)&&(*norien>0)){
1149               printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1150               strcpy1(&filab[266],"G",1);
1151           }
1152       }
1153       
1154       if(strcmp1(&filab[1479],"PHS ")==0){
1155           if((strcmp1(&filab[1484],"L")==0)&&(*norien>0)){
1156               printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1157               strcpy1(&filab[1484],"G",1);
1158           }
1159       }
1160       
1161       if(strcmp1(&filab[1653],"MAXS")==0){
1162           if((strcmp1(&filab[1658],"L")==0)&&(*norien>0)){
1163               printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1164               strcpy1(&filab[1658],"G",1);
1165           }
1166       }   
1167       
1168       if(strcmp1(&filab[2523],"MAXE")==0){
1169           if((strcmp1(&filab[2528],"L")==0)&&(*norien>0)){
1170               printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1171               strcpy1(&filab[1658],"G",1);
1172           }
1173       }   
1174       
1175       /* allocating fields for magnitude and phase information of
1176          displacements and stresses */
1177       
1178       if(strcmp1(&filab[870],"PU")==0){
1179           NNEW(vr,double,mt*nkt);
1180           NNEW(vi,double,mt*nkt);
1181       }
1182       
1183       if(strcmp1(&filab[1479],"PHS")==0){
1184           NNEW(stnr,double,6*nkt);
1185           NNEW(stni,double,6*nkt);
1186       }
1187       
1188       if(strcmp1(&filab[2610],"PRF")==0){
1189           NNEW(fnr,double,mt*nkt);
1190           NNEW(fni,double,mt*nkt);
1191       }
1192       
1193       if(strcmp1(&filab[3915],"PCON")==0){
1194           NNEW(cdnr,double,6*nkt);
1195           NNEW(cdni,double,6*nkt);
1196       }
1197       
1198       if(strcmp1(&filab[1566],"MAXU")==0){
1199           NNEW(vmax,double,4*nkt);
1200       }
1201       
1202       if(strcmp1(&filab[1653],"MAXS")==0){
1203           NNEW(stnmax,double,7*nkt);
1204       }
1205       
1206       if(strcmp1(&filab[2523],"MAXE")==0){
1207           NNEW(eenmax,double,7*nkt);
1208       }
1209       
1210       /* start of output calculations */
1211       
1212       lfin=0;
1213       for(j=0;j<nev;++j){
1214           lint=lfin;
1215           lfin=lfin+neq[1];
1216           
1217           if(mei[3]==1){
1218               if(fwrite(&z[lint],sizeof(double),neq[1],f1)!=neq[1]){
1219                   printf("*ERROR saving data to the eigenvalue file...");
1220                   exit(0);
1221               }
1222           }
1223           
1224           /* check whether the frequency belongs to the requested
1225              interval */
1226           
1227           if(fmin>-0.5){
1228               if(fmin*fmin>d[j]) continue;
1229           }
1230           if(fmax>-0.5){
1231               if(fmax*fmax<d[j]) continue;
1232           }
1233           
1234           if(*nprint>0)FORTRAN(writehe,(&j));
1235           
1236           NNEW(eei,double,6*mi[0]*ne0);
1237           if(*nener==1){
1238               NNEW(stiini,double,6*mi[0]*ne0);
1239               NNEW(enerini,double,mi[0]*ne0);}
1240           
1241 //    memset(&v[0],0.,sizeof(double)*2*mt**nk);
1242           DMEMSET(v,0,2*mt**nk,0.);
1243           
1244           for(k=0;k<neq[1];k+=neq[1]/2){
1245               
1246               for(i=0;i<6*mi[0]*ne0;i++){eme[i]=0.;}
1247               
1248               if(k==0) {kk=0;kkv=0;kk6=0;kkx=0;if(*nprint>0)FORTRAN(writere,());}
1249               else {kk=*nk;kkv=mt**nk;kk6=6**nk;kkx=6*mi[0]**ne;
1250                   if(*nprint>0)FORTRAN(writeim,());}
1251               
1252               /* generating the cyclic MPC's (needed for nodal diameters
1253                  different from 0 */
1254               
1255               for(i=0;i<*nmpc;i++){
1256                   index=ipompc[i]-1;
1257                   /* check whether thermal mpc */
1258                   if(nodempc[3*index+1]==0) continue;
1259                   coefmpcnew[index]=coefmpc[index];
1260                   while(1){
1261                       index=nodempc[3*index+2];
1262                       if(index==0) break;
1263                       index--;
1264                       
1265                       icomplex=0;
1266                       inode=nodempc[3*index];
1267                       if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
1268                           icomplex=atoi(&labmpc[20*i+6]);}
1269                       else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
1270                           for(ij=0;ij<*mcs;ij++){
1271                               lprev=cs[ij*17+13];
1272                               ilength=cs[ij*17+3];
1273                               FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
1274                               if(id!=0){
1275                                   if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
1276                               }
1277                           }
1278                       }
1279                       
1280                       if(icomplex!=0){
1281                           idir=nodempc[3*index+1];
1282                           idof=nactdof[mt*(inode-1)+idir]-1;
1283                           if(idof==-1){xreal=1.;ximag=1.;}
1284                           else{xreal=z[lint+idof];ximag=z[lint+idof+neq[1]/2];}
1285                           if(k==0) {
1286                               if(fabs(xreal)<1.e-30)xreal=1.e-30;
1287                               coefmpcnew[index]=coefmpc[index]*
1288                                   (cs[17*(icomplex-1)+14]+ximag/xreal*cs[17*(icomplex-1)+15]);}
1289                           else {
1290                               if(fabs(ximag)<1.e-30)ximag=1.e-30;
1291                               coefmpcnew[index]=coefmpc[index]*
1292                                   (cs[17*(icomplex-1)+14]-xreal/ximag*cs[17*(icomplex-1)+15]);}
1293                       }
1294                       else{coefmpcnew[index]=coefmpc[index];}
1295                   }
1296               }
1297               
1298               if(*iperturb==0){
1299                   results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1300                     &stx[kkx],elcon,
1301                     nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1302                     norien,orab,ntmat_,t0,t0,ithermal,
1303                     prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1304                     f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
1305                     nodeboun,ndirboun,xboun,nboun,ipompc,
1306                     nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1307                     &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1308                     xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1309                     ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1310                     xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1311                     ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1312                     nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1313                     &ne0,xforc,nforc,thicke,shcon,nshcon,
1314                     sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1315                     mortar,islavact,&cdn[kk6],islavnode,nslavnode,ntie,clearini,
1316                     islavsurf,ielprop,prop);}
1317               else{
1318                   results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1319                     &stx[kkx],elcon,
1320                     nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1321                     norien,orab,ntmat_,t0,t1old,ithermal,
1322                     prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1323                     f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
1324                     nodeboun,ndirboun,xboun,nboun,ipompc,
1325                     nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1326                     &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1327                     xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1328                     ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1329                     xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1330                     ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1331                     nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1332                     &ne0,xforc,nforc,thicke,shcon,nshcon,
1333                     sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1334                     mortar,islavact,&cdn[kk6],islavnode,nslavnode,ntie,clearini,
1335                     islavsurf,ielprop,prop);
1336               }
1337               
1338           }
1339           SFREE(eei);
1340           if(*nener==1){SFREE(stiini);SFREE(enerini);}
1341           
1342           if(strcmp1(&filab[1566],"MAXU")==0){
1343               
1344               /* determining the ray vector; the components of the
1345                  ray vector are the coordinates of the node in node set
1346                  RAY */
1347               
1348               iray=0;
1349               for(i=0;i<*nset;i++){
1350                   if(strcmp1(&set[81*i],"RAYN")==0){
1351                       iray=ialset[istartset[i]-1];
1352                       vray[0]=co[3*iray-3];
1353                       vray[1]=co[3*iray-2];
1354                       vray[2]=co[3*iray-1];
1355                       break;
1356                   }
1357               }
1358               if(iray==0){
1359                   printf("/n*ERROR in arpackcs: no light ray vector/n/n");
1360                   FORTRAN(stop,());
1361               }
1362               
1363               /* initialization */
1364               
1365               for(l1=0;l1<4**nk;l1++){vmax[l1]=0.;}
1366               
1367               /* vector p1 is a point on the rotation axis
1368                  vector p2 is a unit vector along the axis */
1369               
1370               for(l2=0;l2<3;l2++){p1[l2]=cs[5+l2];}
1371               for(l2=0;l2<3;l2++){p2[l2]=cs[8+l2]-p1[l2];}
1372               dd=sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]);
1373               for(l2=0;l2<3;l2++){p2[l2]/=dd;}
1374               
1375               /* determine the time for the worst displacement
1376                  orthogonal to a give light ray vector ; */
1377               
1378               for(l1=0;l1<*nk;l1++){
1379                   
1380                   /*  determining a vector through node (l1+1) and
1381                       orthogonal to the rotation axis */
1382                   
1383                   for(l2=0;l2<3;l2++){q[l2]=co[3*l1+l2]-p1[l2];}
1384                   dd=q[0]*p2[0]+q[1]*p2[1]+q[2]*p2[2];
1385                   for(l2=0;l2<3;l2++){q[l2]-=dd*p2[l2];}
1386                   
1387                   /* determining a vector tan orthogonal to vector q
1388                      and the ray vector */
1389                   
1390                   tan[0]=q[1]*vray[2]-q[2]*vray[1];
1391                   tan[1]=q[2]*vray[0]-q[0]*vray[2];
1392                   tan[2]=q[0]*vray[1]-q[1]*vray[0];
1393                   
1394                   printf("tangent= %" ITGFORMAT ",%e,%e,%e\n",l1,tan[0],tan[1],tan[2]);
1395                   
1396                   worstumax=0.;
1397                   iworsttime=0;
1398                   for(l3=0;l3<360;l3++){
1399                       ctl=cos(l3/constant);
1400                       stl=sin(l3/constant);
1401                       for(l2=1;l2<4;l2++){
1402                           l=mt*l1+l2;
1403                           vl[l2]=ctl*v[l]-stl*v[l+mt**nk];
1404                       }
1405                       
1406                       /* displacement component along the tangent vector
1407                          (no absolute value!) */
1408                       
1409                       dd=vl[1]*tan[0]+vl[2]*tan[1]+vl[3]*tan[2];
1410                       if(dd>worstumax){
1411                           worstumax=dd;
1412                           iworsttime=l3;
1413                       }
1414                   }
1415                   ctl=cos(iworsttime/constant);
1416                   stl=sin(iworsttime/constant);
1417                   for(l2=1;l2<4;l2++){
1418                       l=mt*l1+l2;
1419                       vl[l2]=ctl*v[l]-stl*v[l+mt**nk];
1420                   }
1421                   vmax[4*l1]=1.*iworsttime;
1422                   vmax[4*l1+1]=vl[1];
1423                   vmax[4*l1+2]=vl[2];
1424                   vmax[4*l1+3]=vl[3];
1425                   
1426               }
1427           }
1428           
1429           /* determine the worst principal stress anywhere
1430              in the structure as a function of time; 
1431              the worst principal stress is the maximum
1432              of the absolute value of the principal stresses */
1433           
1434           if(strcmp1(&filab[1653],"MAXS")==0){
1435               
1436               /* determining the set of nodes for the 
1437                  worst principal stress calculation */
1438               
1439               ielset=0;
1440               for(i=0;i<*nset;i++){
1441                   if(strcmp1(&set[81*i],"STRESSDOMAINN")==0){
1442                       ielset=i+1;
1443                       break;
1444                   }
1445               }
1446               if(ielset==0){
1447                   printf("\n*ERROR in arpackcs: no node set for MAXS\n");
1448                   printf("       (must have the name STRESSDOMAIN)\n\n");
1449                   FORTRAN(stop,());
1450               }
1451               
1452               for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
1453                   if(ialset[i1]>0){
1454                       l1=ialset[i1]-1;
1455                       
1456                       worstpsmax=0.;
1457                       for(l3=0;l3<360;l3++){
1458                           ctl=cos(l3/constant);
1459                           stl=sin(l3/constant);
1460                           for(l2=0;l2<6;l2++){
1461                               l=6*l1+l2;
1462                               stnl[l2]=ctl*stn[l]-stl*stn[l+6**nk];
1463                           }
1464                           
1465                           /* determining the eigenvalues */
1466                           
1467                           v1=stnl[0]+stnl[1]+stnl[2];
1468                           v2=stnl[1]*stnl[2]+stnl[0]*stnl[2]+stnl[0]*stnl[1]-
1469                               (stnl[5]*stnl[5]+stnl[4]*stnl[4]+stnl[3]*stnl[3]);
1470                           v3=stnl[0]*(stnl[1]*stnl[2]-stnl[5]*stnl[5])
1471                               -stnl[3]*(stnl[3]*stnl[2]-stnl[4]*stnl[5])
1472                               +stnl[4]*(stnl[3]*stnl[5]-stnl[4]*stnl[1]);
1473                           bb=v2-v1*v1/3.;
1474                           cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1475                           if(fabs(bb)<=1.e-10){
1476                               if(fabs(cc)>1.e-10){
1477                                   al[0]=-pow(cc,(1./3.));
1478                               }else{
1479                                   al[0]=0.;
1480                               }
1481                               al[1]=al[0];
1482                               al[2]=al[0];
1483                           }else{
1484                               cm=2.*sqrt(-bb/3.);
1485                               cn=3.*cc/(cm*bb);
1486                               if(fabs(cn)>1.){
1487                                   if(cn>1.){
1488                                       cn=1.;
1489                                   }else{
1490                                       cn=-1.;
1491                                   }
1492                               }
1493                               tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1494                               al[0]=cm*cos(tt);
1495                               al[1]=cm*cos(tt+2.*pi/3.);
1496                               al[2]=cm*cos(tt+4.*pi/3.);
1497                           }
1498                           for(l2=0;l2<3;l2++){
1499                               al[l2]+=v1/3.;
1500                           }
1501                           dd=fabs(al[0]);
1502                           if(fabs(al[1])>dd) dd=fabs(al[1]);
1503                           if(fabs(al[2])>dd) dd=fabs(al[2]);
1504                           if(dd>worstpsmax){
1505                               worstpsmax=dd;
1506                               stnmax[7*l1]=dd;
1507                               for(l2=1;l2<7;l2++){
1508                                   stnmax[7*l1+l2]=stnl[l2-1];
1509                               }
1510                           }
1511                       }
1512                       
1513                   }else{
1514                       l1=ialset[i1-2]-1;
1515                       do{
1516                           l1=l1-ialset[i1];
1517                           if(l1>=ialset[i1-1]-1) break;
1518                           
1519                           worstpsmax=0.;
1520                           for(l3=0;l3<360;l3++){
1521                               ctl=cos(l3/constant);
1522                               stl=sin(l3/constant);
1523                               for(l2=0;l2<6;l2++){
1524                                   l=6*l1+l2;
1525                                   stnl[l2]=ctl*stn[l]-stl*stn[l+6**nk];
1526                               }
1527                               
1528                               /* determining the eigenvalues */
1529                               
1530                               v1=stnl[0]+stnl[1]+stnl[2];
1531                               v2=stnl[1]*stnl[2]+stnl[0]*stnl[2]+stnl[0]*stnl[1]-
1532                                   (stnl[5]*stnl[5]+stnl[4]*stnl[4]+stnl[3]*stnl[3]);
1533                               v3=stnl[0]*(stnl[1]*stnl[2]-stnl[5]*stnl[5])
1534                                   -stnl[3]*(stnl[3]*stnl[2]-stnl[4]*stnl[5])
1535                                   +stnl[4]*(stnl[3]*stnl[5]-stnl[4]*stnl[1]);
1536                               bb=v2-v1*v1/3.;
1537                               cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1538                               if(fabs(bb)<=1.e-10){
1539                                   if(fabs(cc)>1.e-10){
1540                                       al[0]=-pow(cc,(1./3.));
1541                                   }else{
1542                                       al[0]=0.;
1543                                   }
1544                                   al[1]=al[0];
1545                                   al[2]=al[0];
1546                               }else{
1547                                   cm=2.*sqrt(-bb/3.);
1548                                   cn=3.*cc/(cm*bb);
1549                                   if(fabs(cn)>1.){
1550                                       if(cn>1.){
1551                                           cn=1.;
1552                                       }else{
1553                                           cn=-1.;
1554                                       }
1555                                   }
1556                                   tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1557                                   al[0]=cm*cos(tt);
1558                                   al[1]=cm*cos(tt+2.*pi/3.);
1559                                   al[2]=cm*cos(tt+4.*pi/3.);
1560                               }
1561                               for(l2=0;l2<3;l2++){
1562                                   al[l2]+=v1/3.;
1563                               }
1564                               dd=fabs(al[0]);
1565                               if(fabs(al[1])>dd) dd=fabs(al[1]);
1566                               if(fabs(al[2])>dd) dd=fabs(al[2]);
1567                               if(dd>worstpsmax){
1568                                   worstpsmax=dd;
1569                                   stnmax[7*l1]=dd;
1570                                   for(l2=1;l2<7;l2++){
1571                                       stnmax[7*l1+l2]=stnl[l2-1];
1572                                   }
1573                               }
1574                           }
1575                           
1576                       }while(1);
1577                   }
1578               }
1579           }
1580           
1581           /* determine the worst principal strain anywhere
1582              in the structure as a function of time; 
1583              the worst principal strain is the maximum
1584              of the absolute value of the principal strains,
1585              times its original sign */
1586           
1587     if(strcmp1(&filab[2523],"MAXE")==0){
1588 
1589       /* determining the set of nodes for the 
1590          worst principal strain calculation */
1591 
1592       ielset=0;
1593       for(i=0;i<*nset;i++){
1594         if(strcmp1(&set[81*i],"STRAINDOMAINN")==0){
1595           ielset=i+1;
1596           break;
1597         }
1598       }
1599       if(ielset==0){
1600         printf("\n*ERROR in arpackcs: no node set for MAXE\n");
1601         printf("       (must have the name STRAINDOMAIN)\n\n");
1602         FORTRAN(stop,());
1603       }
1604 
1605       for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
1606         if(ialset[i1]>0){
1607           l1=ialset[i1]-1;
1608 
1609           worstpsmax=0.;
1610           for(l3=0;l3<360;l3++){
1611             ctl=cos(l3/constant);
1612             stl=sin(l3/constant);
1613             for(l2=0;l2<6;l2++){
1614               l=6*l1+l2;
1615               eenl[l2]=ctl*een[l]-stl*een[l+6**nk];
1616             }
1617             
1618             /* determining the eigenvalues */
1619             
1620             v1=eenl[0]+eenl[1]+eenl[2];
1621             v2=eenl[1]*eenl[2]+eenl[0]*eenl[2]+eenl[0]*eenl[1]-
1622               (eenl[5]*eenl[5]+eenl[4]*eenl[4]+eenl[3]*eenl[3]);
1623             v3=eenl[0]*(eenl[1]*eenl[2]-eenl[5]*eenl[5])
1624               -eenl[3]*(eenl[3]*eenl[2]-eenl[4]*eenl[5])
1625               +eenl[4]*(eenl[3]*eenl[5]-eenl[4]*eenl[1]);
1626             bb=v2-v1*v1/3.;
1627             cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1628             if(fabs(bb)<=1.e-10){
1629               if(fabs(cc)>1.e-10){
1630                 al[0]=-pow(cc,(1./3.));
1631               }else{
1632                 al[0]=0.;
1633               }
1634               al[1]=al[0];
1635               al[2]=al[0];
1636             }else{
1637               cm=2.*sqrt(-bb/3.);
1638               cn=3.*cc/(cm*bb);
1639               if(fabs(cn)>1.){
1640                 if(cn>1.){
1641                   cn=1.;
1642                 }else{
1643                   cn=-1.;
1644                 }
1645               }
1646               tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1647               al[0]=cm*cos(tt);
1648               al[1]=cm*cos(tt+2.*pi/3.);
1649               al[2]=cm*cos(tt+4.*pi/3.);
1650             }
1651             for(l2=0;l2<3;l2++){
1652               al[l2]+=v1/3.;
1653             }
1654             dd=fabs(al[0]);
1655             if(fabs(al[1])>dd) dd=fabs(al[1]);
1656             if(fabs(al[2])>dd) dd=fabs(al[2]);
1657             if(dd>worstpsmax){
1658                 worstpsmax=dd;
1659                 eenmax[7*l1]=dd;
1660                 for(l2=1;l2<7;l2++){
1661                     eenmax[7*l1+l2]=eenl[l2-1];
1662                 }
1663             }
1664           }
1665           
1666         }else{
1667           l1=ialset[i1-2]-1;
1668           do{
1669             l1=l1-ialset[i1];
1670             if(l1>=ialset[i1-1]-1) break;
1671 
1672             worstpsmax=0.;
1673             for(l3=0;l3<360;l3++){
1674               ctl=cos(l3/constant);
1675               stl=sin(l3/constant);
1676               for(l2=0;l2<6;l2++){
1677                 l=6*l1+l2;
1678                 eenl[l2]=ctl*een[l]-stl*een[l+6**nk];
1679               }
1680               
1681               /* determining the eigenvalues */
1682               
1683               v1=eenl[0]+eenl[1]+eenl[2];
1684               v2=eenl[1]*eenl[2]+eenl[0]*eenl[2]+eenl[0]*eenl[1]-
1685                 (eenl[5]*eenl[5]+eenl[4]*eenl[4]+eenl[3]*eenl[3]);
1686               v3=eenl[0]*(eenl[1]*eenl[2]-eenl[5]*eenl[5])
1687                 -eenl[3]*(eenl[3]*eenl[2]-eenl[4]*eenl[5])
1688                 +eenl[4]*(eenl[3]*eenl[5]-eenl[4]*eenl[1]);
1689               bb=v2-v1*v1/3.;
1690               cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1691               if(fabs(bb)<=1.e-10){
1692                 if(fabs(cc)>1.e-10){
1693                   al[0]=-pow(cc,(1./3.));
1694                 }else{
1695                   al[0]=0.;
1696                 }
1697                 al[1]=al[0];
1698                 al[2]=al[0];
1699               }else{
1700                 cm=2.*sqrt(-bb/3.);
1701                 cn=3.*cc/(cm*bb);
1702                 if(fabs(cn)>1.){
1703                   if(cn>1.){
1704                     cn=1.;
1705                   }else{
1706                     cn=-1.;
1707                   }
1708                 }
1709                 tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1710                 al[0]=cm*cos(tt);
1711                 al[1]=cm*cos(tt+2.*pi/3.);
1712                 al[2]=cm*cos(tt+4.*pi/3.);
1713               }
1714               for(l2=0;l2<3;l2++){
1715                 al[l2]+=v1/3.;
1716               }
1717               dd=fabs(al[0]);
1718               if(fabs(al[1])>dd) dd=fabs(al[1]);
1719               if(fabs(al[2])>dd) dd=fabs(al[2]);
1720               if(dd>worstpsmax){
1721                   worstpsmax=dd;
1722                   eenmax[7*l1]=dd;
1723                   for(l2=1;l2<7;l2++){
1724                       eenmax[7*l1+l2]=eenl[l2-1];
1725                   }
1726               }
1727             }
1728             
1729           }while(1);
1730         }
1731       }
1732     }
1733     
1734     /* mapping the results to the other sectors */
1735 
1736     for(l=0;l<*nk;l++){inumt[l]=inum[l];}
1737 
1738     icntrl=2;imag=1;
1739 
1740     FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1741 
1742     if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1743       for(l=0;l<mt**nk;l++){vt[l]=v[l];}
1744       for(l=0;l<mt**nk;l++){vt[l+mt**nk*ngraph]=v[l+mt**nk];}}
1745     if(strcmp1(&filab[87],"NT  ")==0)
1746       for(l=0;l<*nk;l++){t1t[l]=t1[l];};
1747     if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1748         for(l=0;l<6**nk;l++){stnt[l]=stn[l];}
1749         for(l=0;l<6**nk;l++){stnt[l+6**nk*ngraph]=stn[l+6**nk];}}
1750     if(strcmp1(&filab[261],"E   ")==0){
1751         for(l=0;l<6**nk;l++){eent[l]=een[l];};
1752         for(l=0;l<6**nk;l++){eent[l+6**nk*ngraph]=een[l+6**nk];}}
1753     if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1754       for(l=0;l<mt**nk;l++){fnt[l]=fn[l];}
1755       for(l=0;l<mt**nk;l++){fnt[l+mt**nk*ngraph]=fn[l+mt**nk];}}
1756     if(strcmp1(&filab[522],"ENER")==0)
1757       for(l=0;l<*nk;l++){enernt[l]=enern[l];};
1758     if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
1759        ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0))){
1760       for(l=0;l<6*mi[0]**ne;l++){stxt[l]=stx[l];}
1761       for(l=0;l<6*mi[0]**ne;l++){stxt[l+6*mi[0]**ne*ngraph]=stx[l+6*mi[0]**ne];}}
1762     if(strcmp1(&filab[2697],"ME  ")==0){
1763         for(l=0;l<6**nk;l++){emnt[l]=emn[l];};
1764         for(l=0;l<6**nk;l++){emnt[l+6**nk*ngraph]=emn[l+6**nk];}}
1765     if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1766        &&(*mortar==1)){
1767         for(l=0;l<6**nk;l++){cdnt[l]=cdn[l];}
1768         for(l=0;l<6**nk;l++){cdnt[l+6**nk*ngraph]=cdn[l+6**nk];}}
1769 
1770     for(jj=0;jj<*mcs;jj++){
1771       ilength=cs[17*jj+3];
1772       is=cs[17*jj+4];
1773       lprev=cs[17*jj+13];
1774       for(i=1;i<is;i++){
1775         
1776         for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
1777         
1778         theta=i*nm*2.*pi/cs[17*jj];
1779         ctl=cos(theta);
1780         stl=sin(theta);
1781         
1782         if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1783           for(l1=0;l1<*nk;l1++){
1784             if(inocs[l1]==jj){
1785 
1786               /* check whether node lies on axis */
1787 
1788               ml1=-l1-1;
1789               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1790               if(id!=0){
1791                   if(ics[lprev+id-1]==ml1){
1792                       for(l2=0;l2<4;l2++){
1793                           l=mt*l1+l2;
1794                           vt[l+mt**nk*i]=v[l];
1795                       }
1796                       continue;
1797                   }
1798               }
1799               for(l2=0;l2<4;l2++){
1800                 l=mt*l1+l2;
1801                 vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
1802               }
1803             }
1804           }
1805         }
1806         
1807         /* imaginary part of the displacements in cylindrical
1808            coordinates */
1809 
1810         if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)){
1811           for(l1=0;l1<*nk;l1++){
1812             if(inocs[l1]==jj){
1813 
1814               /* check whether node lies on axis */
1815 
1816               ml1=-l1-1;
1817               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1818               if(id!=0){
1819                   if(ics[lprev+id-1]==ml1){
1820                       for(l2=0;l2<4;l2++){
1821                           l=mt*l1+l2;
1822                           vt[l+mt**nk*(i+ngraph)]=v[l+mt**nk];
1823                       }
1824                       continue;
1825                   }
1826               }
1827               for(l2=0;l2<4;l2++){
1828                 l=mt*l1+l2;
1829                 vt[l+mt**nk*(i+ngraph)]=stl*v[l]+ctl*v[l+mt**nk];
1830               }
1831             }
1832           }
1833         }
1834         
1835         if(strcmp1(&filab[87],"NT  ")==0){
1836           for(l=0;l<*nk;l++){
1837               if(inocs[l]==jj) t1t[l+*nk*i]=t1[l];
1838           }
1839         }
1840         
1841         if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1842           for(l1=0;l1<*nk;l1++){
1843             if(inocs[l1]==jj){
1844 
1845               /* check whether node lies on axis */
1846 
1847               ml1=-l1-1;
1848               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1849               if(id!=0){
1850                   if(ics[lprev+id-1]==ml1){
1851                       for(l2=0;l2<6;l2++){
1852                           l=6*l1+l2;
1853                           stnt[l+6**nk*i]=stn[l];
1854                       }
1855                       continue;
1856                   }
1857               }
1858               for(l2=0;l2<6;l2++){
1859                 l=6*l1+l2;
1860                 stnt[l+6**nk*i]=ctl*stn[l]-stl*stn[l+6**nk];
1861               }
1862             }
1863           }
1864         }
1865         
1866         /* imaginary part of the stresses in cylindrical
1867            coordinates */
1868         
1869         if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1870           for(l1=0;l1<*nk;l1++){
1871             if(inocs[l1]==jj){
1872 
1873               /* check whether node lies on axis */
1874 
1875               ml1=-l1-1;
1876               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1877               if(id!=0){
1878                   if(ics[lprev+id-1]==ml1){
1879                       for(l2=0;l2<6;l2++){
1880                           l=6*l1+l2;
1881                           stnt[l+6**nk*(i+ngraph)]=stn[l+6**nk];
1882                       }
1883                       continue;
1884                   }
1885               }
1886               for(l2=0;l2<6;l2++){
1887                 l=6*l1+l2;
1888                 stnt[l+6**nk*(i+ngraph)]=stl*stn[l]+ctl*stn[l+6**nk];
1889               }
1890             }
1891           }
1892         }
1893         
1894         if(strcmp1(&filab[261],"E   ")==0){
1895           for(l1=0;l1<*nk;l1++){
1896             if(inocs[l1]==jj){
1897 
1898               /* check whether node lies on axis */
1899 
1900               ml1=-l1-1;
1901               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1902               if(id!=0){
1903                   if(ics[lprev+id-1]==ml1){
1904                       for(l2=0;l2<6;l2++){
1905                           l=6*l1+l2;
1906                           eent[l+6**nk*i]=een[l];
1907                       }
1908                       continue;
1909                   }
1910               }
1911               for(l2=0;l2<6;l2++){
1912                 l=6*l1+l2;
1913                 eent[l+6**nk*i]=ctl*een[l]-stl*een[l+6**nk];
1914               }
1915             }
1916           }
1917         }
1918         
1919         /* imaginary part of the strains in cylindrical
1920            coordinates */
1921         
1922         if(strcmp1(&filab[261],"E   ")==0){
1923           for(l1=0;l1<*nk;l1++){
1924             if(inocs[l1]==jj){
1925 
1926               /* check whether node lies on axis */
1927 
1928               ml1=-l1-1;
1929               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1930               if(id!=0){
1931                   if(ics[lprev+id-1]==ml1){
1932                       for(l2=0;l2<6;l2++){
1933                           l=6*l1+l2;
1934                           eent[l+6**nk*(i+ngraph)]=een[l+6**nk];
1935                       }
1936                       continue;
1937                   }
1938               }
1939               for(l2=0;l2<6;l2++){
1940                 l=6*l1+l2;
1941                 eent[l+6**nk*(i+ngraph)]=stl*een[l]+ctl*een[l+6**nk];
1942               }
1943             }
1944           }
1945         }
1946         
1947         /* real part of the mechanical strains */
1948 
1949         if(strcmp1(&filab[2697],"ME  ")==0){
1950           for(l1=0;l1<*nk;l1++){
1951             if(inocs[l1]==jj){
1952 
1953               /* check whether node lies on axis */
1954 
1955               ml1=-l1-1;
1956               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1957               if(id!=0){
1958                   if(ics[lprev+id-1]==ml1){
1959                       for(l2=0;l2<6;l2++){
1960                           l=6*l1+l2;
1961                           emnt[l+6**nk*i]=emn[l];
1962                       }
1963                       continue;
1964                   }
1965               }
1966               for(l2=0;l2<6;l2++){
1967                 l=6*l1+l2;
1968                 emnt[l+6**nk*i]=ctl*emn[l]-stl*emn[l+6**nk];
1969               }
1970             }
1971           }
1972         }
1973         
1974         /* imaginary part of the mechanical strains in cylindrical
1975            coordinates */
1976         
1977         if(strcmp1(&filab[2697],"ME  ")==0){
1978           for(l1=0;l1<*nk;l1++){
1979             if(inocs[l1]==jj){
1980 
1981               /* check whether node lies on axis */
1982 
1983               ml1=-l1-1;
1984               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1985               if(id!=0){
1986                   if(ics[lprev+id-1]==ml1){
1987                       for(l2=0;l2<6;l2++){
1988                           l=6*l1+l2;
1989                           emnt[l+6**nk*(i+ngraph)]=emn[l+6**nk];
1990                       }
1991                       continue;
1992                   }
1993               }
1994               for(l2=0;l2<6;l2++){
1995                 l=6*l1+l2;
1996                 emnt[l+6**nk*(i+ngraph)]=stl*emn[l]+ctl*emn[l+6**nk];
1997               }
1998             }
1999           }
2000         }
2001         
2002         /* real part of the contact stresses */
2003 
2004         if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2005            &&(*mortar==1)){
2006           for(l1=0;l1<*nk;l1++){
2007             if(inocs[l1]==jj){
2008 
2009               /* check whether node lies on axis */
2010 
2011               ml1=-l1-1;
2012               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2013               if(id!=0){
2014                   if(ics[lprev+id-1]==ml1){
2015                       for(l2=0;l2<6;l2++){
2016                           l=6*l1+l2;
2017                           cdnt[l+6**nk*i]=cdn[l];
2018                       }
2019                       continue;
2020                   }
2021               }
2022               for(l2=0;l2<6;l2++){
2023                 l=6*l1+l2;
2024                 cdnt[l+6**nk*i]=ctl*cdn[l]-stl*cdn[l+6**nk];
2025               }
2026             }
2027           }
2028         }
2029         
2030         /* imaginary part of the contact stresses */
2031         
2032         if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2033            &&(*mortar==1)){
2034           for(l1=0;l1<*nk;l1++){
2035             if(inocs[l1]==jj){
2036 
2037               /* check whether node lies on axis */
2038 
2039               ml1=-l1-1;
2040               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2041               if(id!=0){
2042                   if(ics[lprev+id-1]==ml1){
2043                       for(l2=0;l2<6;l2++){
2044                           l=6*l1+l2;
2045                           cdnt[l+6**nk*(i+ngraph)]=cdn[l+6**nk];
2046                       }
2047                       continue;
2048                   }
2049               }
2050               for(l2=0;l2<6;l2++){
2051                 l=6*l1+l2;
2052                 cdnt[l+6**nk*(i+ngraph)]=stl*cdn[l]+ctl*cdn[l+6**nk];
2053               }
2054             }
2055           }
2056         }
2057         
2058         if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
2059           for(l1=0;l1<*nk;l1++){
2060             if(inocs[l1]==jj){
2061 
2062               /* check whether node lies on axis */
2063 
2064               ml1=-l1-1;
2065               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2066               if(id!=0){
2067                   if(ics[lprev+id-1]==ml1){
2068                       for(l2=0;l2<4;l2++){
2069                           l=mt*l1+l2;
2070                           fnt[l+mt**nk*i]=fn[l];
2071                       }
2072                       continue;
2073                   }
2074               }
2075               for(l2=0;l2<4;l2++){
2076                 l=mt*l1+l2;
2077                 fnt[l+mt**nk*i]=ctl*fn[l]-stl*fn[l+mt**nk];
2078               }
2079             }
2080           }
2081         }
2082         
2083         /* imaginary part of the forces in cylindrical
2084            coordinates */
2085 
2086         if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
2087           for(l1=0;l1<*nk;l1++){
2088             if(inocs[l1]==jj){
2089               
2090               /* check whether node lies on axis */
2091               
2092               ml1=-l1-1;
2093               FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2094               if(id!=0){
2095                 if(ics[lprev+id-1]==ml1){
2096                   for(l2=0;l2<4;l2++){
2097                     l=mt*l1+l2;
2098                     fnt[l+mt**nk*(i+ngraph)]=fn[l+mt**nk];
2099                   }
2100                   continue;
2101                 }
2102               }
2103               for(l2=0;l2<4;l2++){
2104                 l=mt*l1+l2;
2105                 fnt[l+mt**nk*(i+ngraph)]=stl*fn[l]+ctl*fn[l+mt**nk];
2106               }
2107             }
2108           }
2109         }
2110         
2111         if(strcmp1(&filab[522],"ENER")==0){
2112           for(l=0;l<*nk;l++){
2113             if(inocs[l]==jj) enernt[l+*nk*i]=0.;
2114           }
2115         }
2116       }
2117     }
2118         
2119     icntrl=-2;imag=0;
2120 
2121     FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
2122                      &imag,mi,emnt));
2123 
2124     FORTRAN(rectcylvi,(cot,&vt[mt**nk*ngraph],&fnt[mt**nk*ngraph],
2125           &stnt[6**nk*ngraph],qfnt,&eent[6**nk*ngraph],cs,&nkt,&icntrl,
2126           t,filab,&imag,mi,&emnt[6**nk*ngraph]));
2127 
2128     /* determining magnitude and phase angle for the displacements */
2129 
2130     if(strcmp1(&filab[870],"PU")==0){
2131       for(l1=0;l1<nkt;l1++){
2132         for(l2=0;l2<4;l2++){
2133           l=mt*l1+l2;
2134           vreal=vt[l];
2135           vimag=vt[l+mt**nk*ngraph];
2136           vr[l]=sqrt(vreal*vreal+vimag*vimag);
2137           if(fabs(vreal)<1.e-10){
2138             if(vimag>0){vi[l]=90.;}
2139             else{vi[l]=-90.;}
2140           }
2141           else{
2142             vi[l]=atan(vimag/vreal)*constant;
2143             if(vreal<0) vi[l]+=180.;
2144           }
2145         }
2146       }
2147     }
2148 
2149     /* determining magnitude and phase for the stress */
2150 
2151     if(strcmp1(&filab[1479],"PHS")==0){
2152       for(l1=0;l1<nkt;l1++){
2153         for(l2=0;l2<6;l2++){
2154           l=6*l1+l2;
2155           stnreal=stnt[l];
2156           stnimag=stnt[l+6**nk*ngraph];
2157           stnr[l]=sqrt(stnreal*stnreal+stnimag*stnimag);
2158           if(fabs(stnreal)<1.e-10){
2159             if(stnimag>0){stni[l]=90.;}
2160             else{stni[l]=-90.;}
2161           }
2162           else{
2163             stni[l]=atan(stnimag/stnreal)*constant;
2164             if(stnreal<0) stni[l]+=180.;
2165           }
2166         }
2167       }
2168     }
2169 
2170     /* determining magnitude and phase angle for the forces */
2171 
2172     if(strcmp1(&filab[2610],"PRF")==0){
2173       for(l1=0;l1<nkt;l1++){
2174         for(l2=0;l2<4;l2++){
2175           l=mt*l1+l2;
2176           fnreal=fnt[l];
2177           fnimag=fnt[l+mt**nk*ngraph];
2178           fnr[l]=sqrt(fnreal*fnreal+fnimag*fnimag);
2179           if(fabs(fnreal)<1.e-10){
2180             if(fnimag>0){fni[l]=90.;}
2181             else{fni[l]=-90.;}
2182           }
2183           else{
2184             fni[l]=atan(fnimag/fnreal)*constant;
2185             if(fnreal<0) fni[l]+=180.;
2186           }
2187         }
2188       }
2189     }
2190 
2191     /* determining magnitude and phase for the contact stress */
2192 
2193     if((strcmp1(&filab[3915],"PCON")==0)&&(*mortar==1)){
2194       for(l1=0;l1<nkt;l1++){
2195         for(l2=0;l2<6;l2++){
2196           l=6*l1+l2;
2197           cdnreal=cdnt[l];
2198           cdnimag=cdnt[l+6**nk*ngraph];
2199           cdnr[l]=sqrt(cdnreal*cdnreal+cdnimag*cdnimag);
2200           if(fabs(cdnreal)<1.e-10){
2201             if(cdnimag>0){cdni[l]=90.;}
2202             else{cdni[l]=-90.;}
2203           }
2204           else{
2205             cdni[l]=atan(cdnimag/cdnreal)*constant;
2206             if(cdnreal<0) cdni[l]+=180.;
2207           }
2208         }
2209       }
2210     }
2211 
2212     ++*kode;
2213     if(d[j]>=0.){
2214         freq=sqrt(d[j])/6.283185308;
2215     }else{
2216         freq=0.;
2217     }
2218 
2219     if(strcmp1(&filab[1044],"ZZS")==0){
2220         NNEW(neigh,ITG,40*net);
2221         NNEW(ipneigh,ITG,nkt);
2222     }
2223 
2224     frd(cot,&nkt,kont,ipkont,lakont,&net,vt,stnt,inumt,nmethod,
2225             kode,filab,eent,t1t,fnt,&freq,epn,ielmatt,matname,enernt,xstaten,
2226             nstate_,istep,&iinc,ithermal,qfn,&j,&nm,trab,inotrt,
2227             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2228             mi,stxt,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&net,
2229             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emnt,
2230             thicke,jobnamec,output,qfx,cdnt,mortar,cdnr,cdni,nmat);
2231 
2232     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2233 
2234     if(*iperturb!=0){
2235         for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
2236             xstate[k]=xstateini[k];
2237         }         
2238     }
2239   }
2240 
2241   if((fmax>-0.5)&&(fmax*fmax>d[nev-1])){
2242     printf("\n*WARNING: not all frequencies in the requested interval might be found;\nincrease the number of requested frequencies\n");
2243   }
2244 
2245   SFREE(adb);SFREE(aub);SFREE(temp_array);SFREE(coefmpcnew);
2246 
2247   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1653],"MAXS")==0)|| 
2248      (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
2249      (strcmp1(&filab[1044],"ERR ")==0)) 
2250      SFREE(stn);
2251 
2252   SFREE(v);SFREE(fn);SFREE(inum);SFREE(stx);SFREE(z);SFREE(d);
2253 
2254   if((strcmp1(&filab[261],"E   ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) SFREE(een);
2255   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2256   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
2257 
2258   if((strcmp1(&filab[0],"U  ")==0)||(strcmp1(&filab[870],"PU  ")==0)) SFREE(vt);
2259   if(strcmp1(&filab[87],"NT  ")==0) SFREE(t1t);
2260   if((strcmp1(&filab[174],"S   ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
2261      (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stnt);
2262   if(strcmp1(&filab[261],"E   ")==0) SFREE(eent);
2263   if((strcmp1(&filab[348],"RF  ")==0)||(strcmp1(&filab[2610],"PRF ")==0)) SFREE(fnt);
2264   if(strcmp1(&filab[522],"ENER")==0) SFREE(enernt);
2265   if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
2266      ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0))) SFREE(stxt);
2267   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emnt);
2268   if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2269      &&(*mortar==1)) SFREE(cdnt);
2270 
2271   SFREE(cot);SFREE(kont);SFREE(ipkont);SFREE(lakont);SFREE(inumt);SFREE(ielmatt);
2272   if(*ntrans>0){SFREE(inotrt);}
2273 
2274   if(mei[3]==1){
2275       (*nevtot)+=nev;
2276       fclose(f1);
2277   }
2278 
2279   /* end loop over the nodal diameters */
2280 
2281   }
2282 
2283   if(*iperturb!=0){
2284       if(ncont!=0){
2285           *ne=ne0;*nkon=nkon0;
2286           if(*nener==1){
2287               RENEW(ener,double,mi[0]**ne*2);
2288           }
2289           RENEW(ipkon,ITG,*ne);
2290           RENEW(lakon,char,8**ne);
2291           RENEW(kon,ITG,*nkon);
2292           if(*norien>0){
2293               RENEW(ielorien,ITG,mi[2]**ne);
2294           }
2295           RENEW(ielmat,ITG,mi[2]**ne);
2296           SFREE(cg);
2297           SFREE(straight);
2298           SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
2299           SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
2300           SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
2301           SFREE(springarea);SFREE(xmastnor);
2302 
2303           if(*mortar==0){
2304               SFREE(areaslav);
2305           }else if(*mortar==1){
2306               SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
2307               SFREE(islavact);
2308           }
2309       }
2310   }
2311 
2312   SFREE(inocs);SFREE(ielcs);SFREE(xstiff);
2313   SFREE(ipobody);
2314 
2315   if(*nstate_!=0) SFREE(xstateini);
2316 
2317   if(strcmp1(&filab[870],"PU")==0){SFREE(vr);SFREE(vi);}
2318   if(strcmp1(&filab[1479],"PHS")==0){SFREE(stnr);SFREE(stni);}
2319   if(strcmp1(&filab[3915],"PCON")==0){SFREE(cdnr);SFREE(cdni);}
2320   if(strcmp1(&filab[1566],"MAXU")==0){SFREE(vmax);}
2321   if(strcmp1(&filab[1653],"MAXS")==0){SFREE(stnmax);}
2322   if(strcmp1(&filab[2523],"MAXE")==0){SFREE(eenmax);}
2323   if(strcmp1(&filab[2610],"PRF")==0){SFREE(fnr);SFREE(fni);}
2324 
2325   for(i=0;i<6*mi[0]*ne0;i++){eme[i]=0.;}
2326 
2327   if(*iperturb!=0){
2328       mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
2329       mpcinfo[3]=maxlenmpc;
2330   }
2331 
2332   *irowp=irow;*enerp=ener;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
2333   *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;
2334 
2335   *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
2336 
2337   return;
2338 }
2339 
2340 #endif

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