root/src/electromagnetics.c

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

DEFINITIONS

This source file includes following definitions.
  1. electromagnetics

   1 /*     CalculiX - A 3-dimensional finite element program                 */
   2 /*              Copyright (C) 1998-2015 Guido Dhondt                          */
   3 
   4 /*     This program is free software; you can redistribute it and/or     */
   5 /*     modify it under the terms of the GNU General Public License as    */
   6 /*     published by the Free Software Foundation(version 2);    */
   7 /*                    */
   8 
   9 /*     This program is distributed in the hope that it will be useful,   */
  10 /*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
  11 /*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
  12 /*     GNU General Public License for more details.                      */
  13 
  14 /*     You should have received a copy of the GNU General Public License */
  15 /*     along with this program; if not, write to the Free Software       */
  16 /*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */
  17 
  18 #include <stdio.h>
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include "CalculiX.h"
  22 #ifdef SPOOLES
  23    #include "spooles.h"
  24 #endif
  25 #ifdef SGI
  26    #include "sgi.h"
  27 #endif
  28 #ifdef TAUCS
  29    #include "tau.h"
  30 #endif
  31 #ifdef PARDISO
  32    #include "pardiso.h"
  33 #endif
  34 
  35 
  36 void electromagnetics(double **cop, ITG *nk, ITG **konp, ITG **ipkonp, 
  37              char **lakonp,ITG *ne, 
  38              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  39              ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp,
  40              ITG *nmpc, 
  41              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  42              ITG **nelemloadp, char **sideloadp, double *xload,ITG *nload, 
  43              ITG *nactdof, 
  44              ITG **icolp, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, 
  45              ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG *ikboun, 
  46              ITG *ilboun,
  47              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  48              double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
  49              ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
  50              double *t0, double *t1, double *t1old, 
  51              ITG *ithermal,double *prestr, ITG *iprestr, 
  52              double **voldp,ITG *iperturb, double *sti, ITG *nzs,  
  53              ITG *kode,char *filab, 
  54              ITG *idrct, ITG *jmax, ITG *jout, double *timepar,
  55              double *eme,
  56              double *xbounold, double *xforcold, double *xloadold,
  57              double *veold, double *accold,
  58              char *amname, double *amta, ITG *namta, ITG *nam,
  59              ITG *iamforc, ITG **iamloadp,
  60              ITG *iamt1, double *alpha, ITG *iexpl,
  61              ITG *iamboun, double *plicon, ITG *nplicon, double *plkcon,
  62              ITG *nplkcon,
  63              double **xstatep, ITG *npmat_, ITG *istep, double *ttime,
  64              char *matname, double *qaold, ITG *mi,
  65              ITG *isolver, ITG *ncmat_, ITG *nstate_, ITG *iumat,
  66              double *cs, ITG *mcs, ITG *nkon, double **enerp, ITG *mpcinfo,
  67              char *output,
  68              double *shcon, ITG *nshcon, double *cocon, ITG *ncocon,
  69              double *physcon, ITG *nflow, double *ctrl, 
  70              char **setp, ITG *nset, ITG **istartsetp,
  71              ITG **iendsetp, ITG **ialsetp, ITG *nprint, char *prlab,
  72              char *prset, ITG *nener,ITG *ikforc,ITG *ilforc, double *trab, 
  73              ITG *inotr, ITG *ntrans, double **fmpcp, char *cbody,
  74              ITG *ibody, double *xbody, ITG *nbody, double *xbodyold,
  75              ITG *ielprop, double *prop, ITG *ntie, char **tiesetp,
  76              ITG *itpamp, ITG *iviewfile, char *jobnamec, double **tietolp,
  77              ITG *nslavs, double *thicke, ITG *ics, ITG *nalset, ITG *nmpc_,
  78              ITG *nmat, char *typeboun,ITG *iaxial,ITG **idefloadp,ITG *nload_){
  79 
  80   char description[13]="            ",*lakon=NULL,jobnamef[396]="",
  81       *labmpc=NULL,kind1[2]="E",kind2[2]="E",*set=NULL,*tieset=NULL,
  82       cflag[1]=" ",*sideloadref=NULL,*sideload=NULL; 
  83  
  84   ITG *inum=NULL,k,iout=0,icntrl,iinc=0,jprint=0,iit=-1,jnz=0,
  85       icutb=0,istab=0,ifreebody,uncoupled=0,maxfaces,indexe,nope,
  86       iperturb_sav[2],*icol=NULL,*irow=NULL,ielas=0,icmd=0,j,
  87       memmpc_,mpcfree,icascade,maxlenmpc,*nodempc=NULL,*iaux=NULL,
  88       *itg=NULL,*ineighe=NULL,null=0,iactive[3],neqterms,ntflag,
  89       *ieg=NULL,ntg=0,ntr,*kontri=NULL,*nloadtr=NULL,index,
  90       *ipiv=NULL,ntri,newstep,mode=-1,noddiam=-1,nasym=0,
  91       ntrit,*inocs=NULL,inewton=0,*ipobody=NULL,*nacteq=NULL,
  92       *nactdog=NULL,nteq,network,nmastnode,imast,massact[2],
  93       *ipkon=NULL,*kon=NULL,*ielorien=NULL,nmethodact,ne2=0,
  94       *ielmat=NULL,inext,itp=0,symmetryflag=0,inputformat=0,
  95       iitterm=0,ngraph=1,ithermalact=2,*islavact=NULL,
  96       *ipompc=NULL,*ikmpc=NULL,*ilmpc=NULL,i0ref,irref,icref,
  97       *islavnode=NULL,*imastnode=NULL,*nslavnode=NULL,mortar=0,
  98       mt=mi[1]+1,*nactdofinv=NULL, inode,idir,*islavsurf=NULL,
  99       iemchange=0,nzsrad,*mast1rad=NULL,*irowrad=NULL,*icolrad=NULL,
 100       *jqrad=NULL,*ipointerrad=NULL,*integerglob=NULL,im,
 101       mass[2]={0,0}, stiffness=1, buckling=0, rhsi=1, intscheme=0,idiscon=0,
 102       coriolis=0,*ipneigh=NULL,*neigh=NULL,i,icfd=0,id,node,networknode,
 103       iflagact=0,*nodorig=NULL,*ipivr=NULL,*inomat=NULL,*nodface=NULL,
 104       *ipoface=NULL,*istartset=NULL,*iendset=NULL,*ialset=NULL,
 105       *nelemloadref=NULL,*iamloadref=NULL,*idefloadref=NULL,nloadref,
 106       *nelemload=NULL,*iamload=NULL,*idefload=NULL;
 107 
 108   double *stn=NULL,*v=NULL,*een=NULL,cam[5],*epn=NULL,*cdn=NULL,
 109          *f=NULL,*fn=NULL,qa[3]={0.,0.,-1.},qam[2]={0.,0.},dtheta,theta,
 110          err,ram[4]={0.,0.,0.,0.},*springarea=NULL,*h0=NULL,
 111          ram1[2]={0.,0.},ram2[2]={0.,0.},deltmx,*clearini=NULL,
 112          uam[2]={0.,0.},*vini=NULL,*ac=NULL,qa0,qau,ea,ptime,
 113          *t1act=NULL,qamold[2],*xbounact=NULL,*bc=NULL,
 114          *xforcact=NULL,*xloadact=NULL,*fext=NULL,h0scale=1.,
 115          reltime,time,bet=0.,gam=0.,*aux1=NULL,*aux2=NULL,dtime,*fini=NULL,
 116          *fextini=NULL,*veini=NULL,*xstateini=NULL,
 117          *ampli=NULL,*eei=NULL,*t1ini=NULL,*tinc,*tper,*tmin,*tmax,
 118          *xbounini=NULL,*xstiff=NULL,*stx=NULL,*cv=NULL,*cvini=NULL,
 119          *enern=NULL,*coefmpc=NULL,*aux=NULL,*xstaten=NULL,
 120          *enerini=NULL,*emn=NULL,*xmastnor=NULL,
 121          *tarea=NULL,*tenv=NULL,*erad=NULL,*fnr=NULL,*fni=NULL,
 122          *adview=NULL,*auview=NULL,*qfx=NULL,
 123          *qfn=NULL,*co=NULL,*vold=NULL,*fenv=NULL,sigma=0.,
 124          *xbodyact=NULL,*cgr=NULL,dthetaref,*vr=NULL,*vi=NULL,
 125          *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,*fmpc=NULL,*ener=NULL,
 126          *f_cm=NULL,*f_cs=NULL,*tietol=NULL,
 127          *xstate=NULL,*eenmax=NULL,*adrad=NULL,*aurad=NULL,*bcr=NULL,
 128          *emeini=NULL,*doubleglob=NULL,*au=NULL,
 129          *ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
 130          *cdnr=NULL,*cdni=NULL;
 131 
 132 #ifdef SGI
 133   ITG token;
 134 #endif
 135  
 136   /* next line is needed to avoid that elements with negative ipkon
 137      are taken into account in extrapolate.f */
 138 
 139   strcpy1(&filab[2],"C",1);
 140 
 141   if(*nmethod==8){
 142       *nmethod=1;
 143   }else if(*nmethod==9){
 144       *nmethod=4;
 145   }else if(*nmethod==10){
 146       *nmethod=2;
 147   }
 148  
 149   for(k=0;k<3;k++){
 150       strcpy1(&jobnamef[k*132],&jobnamec[k*132],132);
 151   }
 152   
 153   qa0=ctrl[20];qau=ctrl[21];ea=ctrl[23];deltmx=ctrl[26];
 154   i0ref=ctrl[0];irref=ctrl[1];icref=ctrl[3];
 155   
 156   memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
 157   maxlenmpc=mpcinfo[3];
 158   
 159   icol=*icolp;irow=*irowp;co=*cop;vold=*voldp;
 160   ipkon=*ipkonp;lakon=*lakonp;kon=*konp;ielorien=*ielorienp;
 161   ielmat=*ielmatp;ener=*enerp;xstate=*xstatep;
 162   
 163   ipompc=*ipompcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
 164   fmpc=*fmpcp;nodempc=*nodempcp;coefmpc=*coefmpcp;
 165 
 166   set=*setp;istartset=*istartsetp;iendset=*iendsetp;ialset=*ialsetp;
 167   tieset=*tiesetp;tietol=*tietolp;
 168 
 169   nelemload=*nelemloadp;iamload=*iamloadp;idefload=*idefloadp;
 170   sideload=*sideloadp;
 171 
 172   tinc=&timepar[0];
 173   tper=&timepar[1];
 174   tmin=&timepar[2];
 175   tmax=&timepar[3];
 176   
 177   /* invert nactdof */
 178   
 179   NNEW(nactdofinv,ITG,mt**nk);
 180   NNEW(nodorig,ITG,*nk);
 181   FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
 182                          ipkon,lakon,kon,ne));
 183   SFREE(nodorig);
 184   
 185   /* allocating a field for the stiffness matrix */
 186   
 187   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 188   
 189   /* allocating force fields */
 190   
 191   NNEW(f,double,neq[1]);
 192   NNEW(fext,double,neq[1]);
 193   
 194   NNEW(b,double,neq[1]);
 195   NNEW(vini,double,mt**nk);
 196   
 197   NNEW(aux,double,7*maxlenmpc);
 198   NNEW(iaux,ITG,maxlenmpc);
 199   
 200   /* allocating fields for the actual external loading */
 201   
 202   NNEW(xbounact,double,*nboun);
 203   NNEW(xbounini,double,*nboun);
 204   for(k=0;k<*nboun;++k){xbounact[k]=xbounold[k];}
 205   NNEW(xforcact,double,*nforc);
 206   NNEW(xloadact,double,2**nload);
 207   NNEW(xbodyact,double,7**nbody);
 208   /* copying the rotation axis and/or acceleration vector */
 209   for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
 210   
 211   /* assigning the body forces to the elements */ 
 212   
 213   if(*nbody>0){
 214       ifreebody=*ne+1;
 215       NNEW(ipobody,ITG,2*ifreebody**nbody);
 216       for(k=1;k<=*nbody;k++){
 217           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 218                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 219           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 220       }
 221       RENEW(ipobody,ITG,2*(ifreebody-1));
 222   }
 223   
 224   /* for thermal calculations: forced convection and cavity
 225      radiation*/
 226   
 227   if(*ithermal>1){
 228       NNEW(itg,ITG,*nload+3**nflow);
 229       NNEW(ieg,ITG,*nflow);
 230       /* max 6 triangles per face, 4 entries per triangle */
 231       NNEW(kontri,ITG,24**nload);
 232       NNEW(nloadtr,ITG,*nload);
 233       NNEW(nacteq,ITG,4**nk);
 234       NNEW(nactdog,ITG,4**nk);
 235       NNEW(v,double,mt**nk);
 236       FORTRAN(envtemp,(itg,ieg,&ntg,&ntr,sideload,nelemload,
 237                        ipkon,kon,lakon,ielmat,ne,nload,
 238                        kontri,&ntri,nloadtr,nflow,ndirboun,nactdog,
 239                        nodeboun,nacteq,nboun,ielprop,prop,&nteq,
 240                        v,&network,physcon,shcon,ntmat_,co,
 241                        vold,set,nshcon,rhcon,nrhcon,mi,nmpc,nodempc,
 242                        ipompc,labmpc,ikboun,&nasym,iaxial));
 243       SFREE(v);
 244       
 245       if((*mcs>0)&&(ntr>0)){
 246           NNEW(inocs,ITG,*nk);
 247           radcyc(nk,kon,ipkon,lakon,ne,cs,mcs,nkon,ialset,istartset,
 248                iendset,&kontri,&ntri,&co,&vold,&ntrit,inocs,mi);
 249       }
 250       else{ntrit=ntri;}
 251       
 252       nzsrad=100*ntr;
 253       NNEW(mast1rad,ITG,nzsrad);
 254       NNEW(irowrad,ITG,nzsrad);
 255       NNEW(icolrad,ITG,ntr);
 256       NNEW(jqrad,ITG,ntr+1);
 257       NNEW(ipointerrad,ITG,ntr);
 258       
 259       if(ntr>0){
 260           mastructrad(&ntr,nloadtr,sideload,ipointerrad,
 261                       &mast1rad,&irowrad,&nzsrad,
 262                       jqrad,icolrad);
 263       }
 264       
 265       SFREE(ipointerrad);SFREE(mast1rad);
 266       RENEW(irowrad,ITG,nzsrad);
 267       
 268       RENEW(itg,ITG,ntg);
 269       NNEW(ineighe,ITG,ntg);
 270       RENEW(kontri,ITG,4*ntrit);
 271       RENEW(nloadtr,ITG,ntr);
 272       
 273       NNEW(adview,double,ntr);
 274       NNEW(auview,double,2*nzsrad);
 275       NNEW(tarea,double,ntr);
 276       NNEW(tenv,double,ntr);
 277       NNEW(fenv,double,ntr);
 278       NNEW(erad,double,ntr);
 279       
 280       NNEW(ac,double,nteq*nteq);
 281       NNEW(bc,double,nteq);
 282       NNEW(ipiv,ITG,nteq);
 283       NNEW(adrad,double,ntr);
 284       NNEW(aurad,double,2*nzsrad);
 285       NNEW(bcr,double,ntr);
 286       NNEW(ipivr,ITG,ntr);
 287   }
 288   if(*ithermal>1){NNEW(qfx,double,3*mi[0]**ne);}
 289   
 290   /* allocating a field for the instantaneous amplitude */
 291   
 292   NNEW(ampli,double,*nam);
 293   
 294   NNEW(fini,double,neq[1]);
 295   
 296   /* allocating fields for nonlinear dynamics */
 297   
 298   if(*nmethod==4){
 299       mass[0]=1;
 300       mass[1]=1;
 301       NNEW(aux2,double,neq[1]);
 302       NNEW(fextini,double,neq[1]);
 303       NNEW(cv,double,neq[1]);
 304       NNEW(cvini,double,neq[1]);
 305       NNEW(veini,double,mt**nk);
 306       NNEW(adb,double,neq[1]);
 307       NNEW(aub,double,nzs[1]);
 308   }
 309   
 310   qa[0]=qaold[0];
 311   qa[1]=qaold[1];
 312   
 313   /* normalizing the time */
 314   
 315   FORTRAN(checktime,(itpamp,namta,tinc,ttime,amta,tmin,&inext,&itp,istep));
 316   dtheta=(*tinc)/(*tper);
 317   dthetaref=dtheta;
 318   if(dtheta<=1.e-6){
 319       printf("\n *ERROR in nonlingeo\n");
 320       printf(" increment size smaller than one millionth of step size\n");
 321       printf(" increase increment size\n\n");
 322   }
 323   *tmin=*tmin/(*tper);
 324   *tmax=*tmax/(*tper);
 325   theta=0.;
 326   
 327   /* calculating an initial flux norm */
 328   
 329   if(*ithermal!=2){
 330       if(qau>1.e-10){qam[0]=qau;}
 331       else if(qa0>1.e-10){qam[0]=qa0;}
 332       else if(qa[0]>1.e-10){qam[0]=qa[0];}
 333       else {qam[0]=1.e-2;}
 334   }
 335   if(*ithermal>1){
 336       if(qau>1.e-10){qam[1]=qau;}
 337       else if(qa0>1.e-10){qam[1]=qa0;}
 338       else if(qa[1]>1.e-10){qam[1]=qa[1];}
 339       else {qam[1]=1.e-2;}
 340   }
 341   
 342   
 343   /*********************************************************************/
 344   
 345   /* calculating the initial quasi-static magnetic intensity due to 
 346      the coil current */
 347   
 348   /*********************************************************************/
 349   
 350   /* calculate the current density in the coils
 351 
 352      in this section nload, nforc, nbody and nam are set to zero; the
 353      electrical potential is supposed to be given (in the form of a
 354      temperature), the current is calculated (in the form of heat
 355      flux) by thermal analogy  */
 356   
 357   reltime=1.;
 358   time=0.;
 359   dtime=0.;
 360   ithermalact=2;
 361 
 362   nmethodact=1;
 363   massact[0]=0;
 364   massact[1]=0;
 365 
 366   if(*ithermal<=1){
 367       NNEW(qfx,double,3*mi[0]**ne);
 368       NNEW(t0,double,*nk);
 369   }
 370   if(strcmp1(&filab[3567],"ECD ")==0){NNEW(qfn,double,3**nk);}
 371   
 372   /* the coil current is assumed to be applied at once, i.e. as 
 373      step loading; the calculation, however, is a quasi-static
 374      calculation */
 375 
 376   FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,&null,xloadold,xload,
 377               xloadact,iamload,&null,ibody,xbody,&null,xbodyold,xbodyact,
 378               t1old,t1,t1act,iamt1,nk,amta,
 379               namta,&null,ampli,&time,&reltime,ttime,&dtime,&ithermalact,nmethod,
 380               xbounold,xboun,xbounact,iamboun,nboun,
 381               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 382               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 383               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 384               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
 385   
 386   cam[0]=0.;cam[1]=0.;
 387   
 388   /* deactivating all elements except the shells */
 389   
 390   for(i=0;i<*ne;i++){
 391       if(strcmp1(&lakon[8*i+6],"L")!=0){
 392           ipkon[i]=-ipkon[i]-2;
 393       }
 394   }
 395   
 396   remastruct(ipompc,&coefmpc,&nodempc,nmpc,
 397              &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
 398              labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
 399              kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
 400              neq,nzs,&nmethodact,&f,&fext,&b,&aux2,&fini,&fextini,
 401              &adb,&aub,&ithermalact,iperturb,mass,mi,iexpl,&mortar,
 402              typeboun,&cv,&cvini,&iit);
 403   
 404   /* invert nactdof */
 405   
 406   SFREE(nactdofinv);
 407   NNEW(nactdofinv,ITG,mt**nk);
 408   NNEW(nodorig,ITG,*nk);
 409   FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
 410                          ipkon,lakon,kon,ne));
 411   SFREE(nodorig);
 412   
 413   iout=-1;
 414   
 415   NNEW(fn,double,mt**nk);
 416   NNEW(inum,ITG,*nk);
 417   NNEW(v,double,mt**nk);
 418 
 419   results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 420           elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 421           ielorien,norien,orab,ntmat_,t0,t1act,&ithermalact,
 422           prestr,iprestr,filab,eme,emn,een,iperturb,f,fn,nactdof,&iout,
 423           qa,vold,b,nodeboun,ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,
 424           labmpc,nmpc,&nmethodact,cam,&neq[1],veold,accold,&bet,
 425           &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 426           xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 427           ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,emeini,xstaten,
 428           eei,enerini,alcon,nalcon,set,nset,istartset,iendset,
 429           ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 430           nelemload,&null,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
 431           ne,xforc,&null,thicke,shcon,nshcon,
 432           sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 433           &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 434           islavsurf,ielprop,prop);
 435   
 436   SFREE(fn);SFREE(inum);SFREE(v);
 437   
 438   iout=1;
 439   
 440   NNEW(ad,double,neq[1]);
 441   NNEW(au,double,nzs[1]);
 442   
 443   FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounold,nboun,
 444                     ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
 445                     &null,nelemload,sideload,xloadact,&null,xbodyact,ipobody,
 446                     &null,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
 447                     &nmethodact,ikmpc,ilmpc,ikboun,ilboun,
 448                     elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
 449                     ielmat,ielorien,norien,orab,ntmat_,
 450                     t0,t1act,&ithermalact,prestr,iprestr,vold,iperturb,sti,
 451                     nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 452                     xstiff,npmat_,&dtime,matname,mi,
 453                     ncmat_,massact,&stiffness,&buckling,&rhsi,&intscheme,
 454                     physcon,shcon,nshcon,alcon,nalcon,ttime,&time,istep,&iinc,
 455                     &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
 456                     xstateini,xstate,thicke,integerglob,doubleglob,
 457                     tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 458                     pmastsurf,&mortar,clearini,ielprop,prop));
 459   
 460   if(nmethodact==0){
 461       
 462       /* error occurred in mafill: storing the geometry in frd format */
 463       
 464       ++*kode;
 465       
 466       ptime=*ttime+time;
 467       frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,&nmethodact,
 468           kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
 469           nstate_,istep,&iinc,&ithermalact,qfn,&mode,&noddiam,trab,inotr,
 470           ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 471           mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 472           cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 473           thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 474       
 475       FORTRAN(stop,());
 476       
 477   }
 478   
 479   for(k=0;k<neq[1];++k){
 480       b[k]=fext[k]-f[k];
 481   }
 482   
 483   if(*isolver==0){
 484 #ifdef SPOOLES
 485       spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
 486               &symmetryflag,&inputformat,&nzs[2]);
 487 #else
 488       printf("*ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
 489       FORTRAN(stop,());
 490 #endif
 491   }
 492   else if((*isolver==2)||(*isolver==3)){
 493       preiter(ad,&au,b,&icol,&irow,&neq[1],&nzs[1],isolver,iperturb);
 494   }
 495   else if(*isolver==4){
 496 #ifdef SGI
 497       token=1;
 498       sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],token);
 499 #else
 500       printf("*ERROR in nonlingeo: the SGI library is not linked\n\n");
 501       FORTRAN(stop,());
 502 #endif
 503   }
 504   else if(*isolver==5){
 505 #ifdef TAUCS
 506       tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[1],&nzs[1]);
 507 #else
 508       printf("*ERROR in nonlingeo: the TAUCS library is not linked\n\n");
 509       FORTRAN(stop,());
 510 #endif
 511   }
 512   else if(*isolver==7){
 513 #ifdef PARDISO
 514       pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
 515                    &symmetryflag,&inputformat,jq,&nzs[2]);
 516 #else
 517       printf("*ERROR in nonlingeo: the PARDISO library is not linked\n\n");
 518       FORTRAN(stop,());
 519 #endif
 520   }
 521 
 522   SFREE(au);SFREE(ad);
 523 
 524   NNEW(v,double,mt**nk);
 525   memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
 526   
 527   NNEW(fn,double,mt**nk);
 528   
 529   NNEW(inum,ITG,*nk);
 530   results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 531           elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 532           ielorien,norien,orab,ntmat_,t0,t1act,&ithermalact,
 533           prestr,iprestr,filab,eme,emn,een,iperturb,
 534           f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 535           ndirboun,xbounact,nboun,ipompc,
 536           nodempc,coefmpc,labmpc,nmpc,&nmethodact,cam,&neq[1],veold,accold,
 537           &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 538           xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 539           &icmd,ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,
 540           emeini,xstaten,eei,enerini,alcon,nalcon,set,nset,istartset,
 541           iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 542           fmpc,nelemload,&null,ikmpc,ilmpc,istep,&iinc,springarea,
 543           &reltime,ne,xforc,&null,thicke,shcon,nshcon,
 544           sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 545           &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 546           islavsurf,ielprop,prop);
 547   
 548   memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
 549   
 550   /* reactivating the non-shell elements (for mesh output purposes) 
 551      deactivating the initial temperature for the non-shell nodes */
 552   
 553   for(i=0;i<*ne;i++){
 554       if(strcmp1(&lakon[8*i+6],"L")!=0){
 555           ipkon[i]=-ipkon[i]-2;
 556           indexe=ipkon[i];
 557           if(strcmp1(&lakon[8*i+3],"6")==0){nope=6;}
 558           else if(strcmp1(&lakon[8*i+3],"8")==0){nope=8;}
 559           else if(strcmp1(&lakon[8*i+3],"1")==0){nope=15;}
 560           else{nope=20;}
 561           for(j=0;j<nope;j++){
 562               node=kon[indexe+j];
 563               v[mt*(node-1)]=0.;
 564           }
 565       }
 566   }
 567 
 568   /* deactivating the output of temperatures */
 569 
 570   if(strcmp1(&filab[87],"NT  ")==0){
 571       ntflag=1;
 572       strcpy1(&filab[87],"    ",4);
 573   }else{ntflag=0;}
 574 
 575   /* createinum is called in order to store the nodes and elements
 576      of the complete structure, not only of the coil */
 577 
 578   FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
 579           nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat));
 580 
 581   ++*kode;
 582   if(*mcs!=0){
 583       ptime=*ttime+time;
 584       frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,&nmethodact,kode,filab,een,
 585              t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
 586              nstate_,istep,&iinc,iperturb,ener,mi,output,&ithermalact,qfn,
 587              ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
 588              norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,ne,
 589              cdn,&mortar,nmat);
 590   }else{
 591       
 592       ptime=*ttime+time;
 593       frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,&nmethodact,
 594           kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
 595           nstate_,istep,&iinc,&ithermalact,qfn,&mode,&noddiam,trab,inotr,
 596           ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 597           mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 598           cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 599           thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
 600       
 601   }
 602   SFREE(inum);SFREE(v);SFREE(fn);
 603 
 604   /* reactivating the temperature output, if previously deactivated */
 605 
 606   if(ntflag==1){
 607       strcpy1(&filab[87],"NT  ",4);
 608   }
 609 
 610   NNEW(inomat,ITG,*nk);
 611 
 612   /* calculating the magnetic intensity caused by the current */
 613 
 614   FORTRAN(assigndomtonodes,(ne,lakon,ipkon,kon,ielmat,inomat,
 615                             elcon,ncmat_,ntmat_,mi,&ne2));
 616 
 617   NNEW(h0,double,3**nk);
 618 
 619   biosav(ipkon,kon,lakon,ne,co,qfx,h0,mi,inomat,nk);
 620 
 621   if(*ithermal<=1){SFREE(qfx);SFREE(t0);}
 622   if(strcmp1(&filab[3567],"ECD ")==0)SFREE(qfn);
 623   
 624   /* deactivating the shell elements */
 625   
 626   for(i=0;i<*ne;i++){
 627       if(strcmp1(&lakon[8*i+6],"L")==0){
 628           ipkon[i]=-ipkon[i]-2;
 629       }
 630   }
 631   
 632 /**************************************************************/
 633 /* creating connecting MPC's between the domains              */
 634 /**************************************************************/
 635 
 636 /* creating contact ties between the domains */
 637 
 638   if(*istep==1){
 639       
 640       NNEW(nodface,ITG,5*6**ne);
 641       NNEW(ipoface,ITG,*nk);
 642       
 643       RENEW(set,char,81*(*nset+3));
 644       RENEW(istartset,ITG,*nset+3);
 645       RENEW(iendset,ITG,*nset+3);
 646       RENEW(ialset,ITG,*nalset+6**ne);
 647       RENEW(tieset,char,243*(*ntie+5));
 648       RENEW(tietol,double,3*(*ntie+5));
 649       
 650       FORTRAN(createtiedsurfs,(nodface,ipoface,set,istartset,
 651               iendset,ialset,tieset,inomat,ne,ipkon,lakon,kon,ntie,
 652               tietol,nalset,nk,nset,iactive));
 653       
 654       SFREE(nodface);SFREE(ipoface);
 655       RENEW(set,char,81**nset);
 656       RENEW(istartset,ITG,*nset);
 657       RENEW(iendset,ITG,*nset);
 658       RENEW(ialset,ITG,*nalset);
 659       RENEW(tieset,char,243**ntie);
 660       RENEW(tietol,double,3**ntie);
 661       
 662       /* tied contact constraints: generate appropriate MPC's */
 663       
 664       tiedcontact(ntie,tieset,nset,set,istartset,iendset,ialset,
 665        lakon,ipkon,kon,tietol,nmpc, &mpcfree,&memmpc_,
 666        &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 667        ithermal,co,vold,&icfd,nmpc_,mi,nk,istep,ikboun,nboun,
 668        kind1,kind2);
 669 
 670       /* mapping h0 from the phi domain onto the border of
 671          the A and A-V domains */
 672 
 673       FORTRAN(calch0interface,(nmpc,ipompc,nodempc,coefmpc,h0));
 674 
 675   }
 676   
 677 /**************************************************************/
 678 /* creating the A.n MPC                                       */
 679 /**************************************************************/
 680 
 681   /* identifying the interfaces between the A and A-V domains
 682      and the phi-domain */
 683 
 684   FORTRAN(generateeminterfaces,(istartset,iendset,
 685           ialset,iactive,ipkon,lakon,kon,ikmpc,nmpc,&maxfaces));
 686   
 687   for(i=1;i<3;i++){
 688       imast=iactive[i];
 689       if(imast==0) continue;
 690 
 691       /* determining the normals on the face */
 692 
 693       NNEW(imastnode,ITG,8*maxfaces);
 694       NNEW(xmastnor,double,3*8*maxfaces);
 695 
 696       FORTRAN(normalsoninterface,(istartset,iendset,
 697        ialset,&imast,ipkon,kon,lakon,imastnode,&nmastnode,
 698        xmastnor,co));
 699 
 700       /* enlarging the fields for MPC's */
 701 
 702       *nmpc_=*nmpc_+nmastnode;
 703       RENEW(ipompc,ITG,*nmpc_);
 704       RENEW(labmpc,char,20**nmpc_+1);
 705       RENEW(ikmpc,ITG,*nmpc_);
 706       RENEW(ilmpc,ITG,*nmpc_);
 707       RENEW(fmpc,double,*nmpc_);
 708       
 709       /* determining the maximum number of terms;
 710          expanding nodempc and coefmpc to accommodate
 711          those terms */
 712       
 713       neqterms=3*nmastnode;
 714       index=memmpc_;
 715       (memmpc_)+=neqterms;
 716       RENEW(nodempc,ITG,3*memmpc_);
 717       RENEW(coefmpc,double,memmpc_);
 718       for(k=index;k<memmpc_;k++){
 719           nodempc[3*k-1]=k+1;
 720       }
 721       nodempc[3*memmpc_-1]=0;
 722 
 723       /* creating the A.n MPC's */
 724 
 725       FORTRAN(createinterfacempcs,(imastnode,xmastnor,&nmastnode,
 726               ikmpc,ilmpc,nmpc,ipompc,nodempc,coefmpc,labmpc,&mpcfree,
 727               ikboun,nboun));
 728 
 729       SFREE(imastnode);SFREE(xmastnor);
 730   }
 731   
 732   /* determining the new matrix structure */
 733       
 734   remastructem(ipompc,&coefmpc,&nodempc,nmpc,
 735              &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
 736              labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
 737              kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
 738              neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
 739              &adb,&aub,ithermal,iperturb,mass,mi,ielmat,elcon,
 740              ncmat_,ntmat_,inomat);
 741 
 742 /**************************************************************/
 743 /* starting the loop over the increments                      */
 744 /**************************************************************/
 745   
 746   /* saving the distributed loads (volume heating will be
 747      added because of Joule heating) */
 748 
 749   if(*ithermal==3){
 750       nloadref=*nload;
 751       NNEW(nelemloadref,ITG,2**nload);
 752       if(*nam>0) NNEW(iamloadref,ITG,2**nload);
 753       NNEW(idefloadref,ITG,*nload);
 754       NNEW(sideloadref,char,20**nload);
 755       
 756       memcpy(&nelemloadref[0],&nelemload[0],sizeof(ITG)*2**nload);
 757       if(*nam>0) memcpy(&iamloadref[0],&iamload[0],sizeof(ITG)*2**nload);
 758       memcpy(&idefloadref[0],&idefload[0],sizeof(ITG)**nload);
 759       memcpy(&sideloadref[0],&sideload[0],sizeof(char)*20**nload);
 760       
 761       /* generating new fields; ne2 is the number of elements
 762          in domain 2 = A,V-domain (the only domain with heating) */
 763       
 764       (*nload_)+=ne2;
 765       RENEW(nelemload,ITG,2**nload_);
 766       if(*nam>0) RENEW(iamload,ITG,2**nload_);
 767       RENEW(idefload,ITG,*nload_);
 768       RENEW(xloadact,double,2**nload_);
 769       RENEW(sideload,char,20**nload_);
 770   }
 771 
 772   if((*ithermal==1)||(*ithermal>=3)){
 773       NNEW(t1ini,double,*nk);
 774       NNEW(t1act,double,*nk);
 775       for(k=0;k<*nk;++k){t1act[k]=t1old[k];}
 776   }
 777   
 778   newstep=1;
 779   
 780   while(1.-theta>1.e-6){
 781       
 782       if(icutb==0){
 783           
 784           /* previous increment converged: update the initial values */
 785           
 786           iinc++;
 787           jprint++;
 788           
 789           /* vold is copied into vini */
 790           
 791           memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
 792           
 793           for(k=0;k<*nboun;++k){xbounini[k]=xbounact[k];}
 794           if((*ithermal==1)||(*ithermal>=3)){
 795               for(k=0;k<*nk;++k){t1ini[k]=t1act[k];}
 796           }
 797           for(k=0;k<neq[1];++k){
 798               fini[k]=f[k];
 799           }
 800           if(*nmethod==4){
 801               for(k=0;k<mt**nk;++k){
 802                   veini[k]=veold[k];
 803               }
 804               for(k=0;k<neq[1];++k){
 805                   fextini[k]=fext[k];
 806               }
 807           }
 808       }
 809       
 810       /* check for max. # of increments */
 811       
 812       if(iinc>*jmax){
 813           printf(" *ERROR: max. # of increments reached\n\n");
 814           FORTRAN(stop,());
 815       }
 816       printf(" increment %" ITGFORMAT " attempt %" ITGFORMAT " \n",iinc,icutb+1);
 817       printf(" increment size= %e\n",dtheta**tper);
 818       printf(" sum of previous increments=%e\n",theta**tper);
 819       printf(" actual step time=%e\n",(theta+dtheta)**tper);
 820       printf(" actual total time=%e\n\n",*ttime+(theta+dtheta)**tper);
 821       
 822       printf(" iteration 1\n\n");
 823       
 824       qamold[0]=qam[0];
 825       qamold[1]=qam[1];
 826       
 827       /* determining the actual loads at the end of the new increment*/
 828       
 829       reltime=theta+dtheta;
 830       time=reltime**tper;
 831       dtime=dtheta**tper;
 832 
 833       /* restoring the distributed loading before adding the
 834          Joule heating */
 835 
 836       if(*ithermal==3){
 837           *nload=nloadref;
 838           DMEMSET(nelemload,0,2**nload_,0);
 839           memcpy(&nelemload[0],&nelemloadref[0],sizeof(ITG)*2**nload);
 840           if(*nam>0){
 841               DMEMSET(iamload,0,2**nload_,0);
 842               memcpy(&iamload[0],&iamloadref[0],sizeof(ITG)*2**nload);
 843           }
 844           DMEMSET(idefload,0,*nload_,0);memcpy(&idefload[0],&idefloadref[0],sizeof(ITG)**nload);
 845           DMEMSET(xloadact,0,2**nload_,0.);
 846           DMEMSET(sideload,0,'\0',0.);memcpy(&sideload[0],&sideloadref[0],sizeof(char)*20**nload);
 847       }
 848       
 849       /* determining the actual loading */
 850 
 851       for(i=0;i<3**nk;i++){h0[i]/=h0scale;}
 852       FORTRAN(tempload_em,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
 853               xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
 854               t1old,t1,t1act,iamt1,nk,amta,
 855               namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
 856               xbounold,xboun,xbounact,iamboun,nboun,
 857               nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 858               co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 859               ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 860               iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc,
 861               &h0scale,inomat));
 862       for(i=0;i<3**nk;i++){h0[i]*=h0scale;}
 863 
 864       for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
 865       if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,bcr,ipivr,
 866        ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,ntmat_,vold,
 867        shcon,nshcon,ipkon,kon,co,
 868        kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
 869        nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
 870        cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
 871        ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
 872        nactdog,nacteq,nodeboun,ndirboun,&network,
 873        rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
 874        ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,nset,
 875        ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,iamload,
 876        jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
 877       }
 878       
 879       /* prediction of the next solution (only for temperature) */
 880       
 881       NNEW(v,double,mt**nk);
 882       
 883 //      if(*ithermal>2){
 884           prediction_em(uam,nmethod,&bet,&gam,&dtime,ithermal,nk,veold,v,
 885                  &iinc,&idiscon,vold,nactdof,mi);
 886 //      }
 887       
 888       NNEW(fn,double,mt**nk);
 889       
 890       iout=-1;
 891       iperturb_sav[0]=iperturb[0];
 892       iperturb_sav[1]=iperturb[1];
 893       
 894       /* first iteration in first increment: heat tangent */
 895       
 896       NNEW(inum,ITG,*nk);
 897       resultsinduction(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 898               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 899               ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
 900               prestr,iprestr,filab,eme,emn,een,iperturb,
 901               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 902               ndirboun,xbounact,nboun,ipompc,
 903               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 904               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 905               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 906               &icmd,ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,
 907               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 908               iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 909               fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 910               &reltime,ne,xforc,nforc,thicke,shcon,nshcon,
 911               sideload,xloadact,xloadold,&icfd,inomat,h0,islavnode,
 912               nslavnode,ntie,ielprop,prop,iactive);
 913       SFREE(inum);
 914      
 915       /* the calculation of the electromagnetic fields is (quasi)linear,
 916          i.e. the solution of the equations is the fields;
 917          only the temperature calculation is nonlinear,
 918          i.e. the solution of the equations is a differential temperature */
 919  
 920       memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
 921       
 922       iout=0;
 923       
 924       SFREE(fn);SFREE(v);
 925       
 926       /***************************************************************/
 927       /* iteration counter and start of the loop over the iterations */
 928       /***************************************************************/
 929       
 930       iit=1;
 931       icntrl=0;
 932       ctrl[0]=i0ref;ctrl[1]=irref;ctrl[3]=icref;
 933       
 934       while(icntrl==0){
 935           
 936           if(iit!=1){
 937               
 938               printf(" iteration %" ITGFORMAT "\n\n",iit);
 939 
 940               /* restoring the distributed loading before adding the
 941                  Joule heating */
 942 
 943               if(*ithermal==3){
 944                   *nload=nloadref;
 945                   DMEMSET(nelemload,0,2**nload_,0);
 946                   memcpy(&nelemload[0],&nelemloadref[0],sizeof(ITG)*2**nload);
 947                   if(*nam>0){
 948                       DMEMSET(iamload,0,2**nload_,0);
 949                       memcpy(&iamload[0],&iamloadref[0],sizeof(ITG)*2**nload);
 950                   }
 951                   DMEMSET(idefload,0,*nload_,0);memcpy(&idefload[0],&idefloadref[0],sizeof(ITG)**nload);
 952                   DMEMSET(xloadact,0,2**nload_,0.);
 953                   DMEMSET(sideload,0,20**nload_,'\0');memcpy(&sideload[0],&sideloadref[0],sizeof(char)*20**nload);
 954               }
 955       
 956               /* determining the actual loading */
 957               
 958               for(i=0;i<3**nk;i++){h0[i]/=h0scale;}
 959               FORTRAN(tempload_em,(xforcold,xforc,xforcact,iamforc,nforc,
 960                xloadold,xload,
 961                xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
 962                t1old,t1,t1act,iamt1,nk,amta,
 963                namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
 964                xbounold,xboun,xbounact,iamboun,nboun,
 965                nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
 966                co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
 967                ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
 968                iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc,
 969                &h0scale,inomat));
 970               for(i=0;i<3**nk;i++){h0[i]*=h0scale;}
 971 
 972               for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
 973               if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,
 974                 bcr,ipivr,ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,
 975                 ntmat_,vold,shcon,nshcon,ipkon,kon,co,
 976                 kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
 977                 nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
 978                 cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
 979                 ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
 980                 nactdog,nacteq,nodeboun,ndirboun,&network,
 981                 rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
 982                 ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,
 983                 nset,ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,
 984                 iamload,jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
 985               }
 986               
 987           }
 988           
 989           /* add Joule heating  */
 990 
 991           if(*ithermal==3){
 992               FORTRAN(jouleheating,(ipkon,lakon,kon,co,elcon,nelcon,
 993                   mi,ne,sti,ielmat,nelemload,sideload,xloadact,nload,nload_,
 994                   iamload,nam,idefload,ncmat_,ntmat_,
 995                   alcon,nalcon,ithermal,vold,t1));
 996           }
 997 
 998           if(*ithermal==3){
 999               for(k=0;k<*nk;++k){t1act[k]=vold[mt*k];}
1000           }
1001               
1002           /* calculating the local stiffness matrix and external loading */
1003           
1004           NNEW(ad,double,neq[1]);
1005           NNEW(au,double,nzs[1]);
1006           
1007           FORTRAN(mafillem,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
1008                   xbounact,nboun,
1009                   ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1010                   nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1011                   nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
1012                   nmethod,ikmpc,ilmpc,ikboun,ilboun,
1013                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1014                   ielmat,ielorien,norien,orab,ntmat_,
1015                   t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1016                   nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1017                   xstiff,npmat_,&dtime,matname,mi,
1018                   ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
1019                   physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
1020                   &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
1021                   xstateini,xstate,thicke,integerglob,doubleglob,
1022                   tieset,istartset,iendset,ialset,ntie,&nasym,iactive,h0,
1023                   pslavsurf,pmastsurf,&mortar,clearini,ielprop,prop));
1024               
1025               iperturb[0]=iperturb_sav[0];
1026               iperturb[1]=iperturb_sav[1];
1027 
1028           /* calculating the residual (f is only for the temperature
1029              nonzero) */
1030 
1031           calcresidual_em(nmethod,neq,b,fext,f,iexpl,nactdof,aux1,aux2,vold,
1032             vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1033             islavnode,nslavnode,&mortar,ntie,f_cm,f_cs,mi,
1034             nzs,&nasym,ithermal);
1035           
1036           newstep=0;
1037           
1038           if(*nmethod==0){
1039               
1040               /* error occurred in mafill: storing the geometry in frd format */
1041               
1042               *nmethod=0;
1043               ++*kode;
1044               NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
1045               
1046               ptime=*ttime+time;
1047               frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
1048                   kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1049                   nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1050                   ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1051                   mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1052                   cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1053                   thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1054               
1055           }
1056           
1057           /* implicit step (static or dynamic */
1058           
1059           if(*nmethod==4){
1060               
1061               /* mechanical part */
1062               
1063               if(*ithermal!=2){
1064                   for(k=0;k<neq[0];++k){
1065                       ad[k]=adb[k]/dtime+ad[k];
1066                   }
1067                   for(k=0;k<nzs[0];++k){
1068                       au[k]=aub[k]/dtime+au[k];
1069                   }
1070                   
1071                   /* upper triangle of asymmetric matrix */
1072                   
1073                   if(nasym>0){
1074                       for(k=nzs[2];k<nzs[2]+nzs[0];++k){
1075                           au[k]=aub[k]/dtime+au[k];
1076                       }
1077                   }
1078               }
1079               
1080               /* thermal part */
1081               
1082               if(*ithermal>1){
1083                   for(k=neq[0];k<neq[1];++k){
1084                       ad[k]=adb[k]/dtime+ad[k];
1085                   }
1086                   for(k=nzs[0];k<nzs[1];++k){
1087                       au[k]=aub[k]/dtime+au[k];
1088                   }
1089                   
1090                   /* upper triangle of asymmetric matrix */
1091                   
1092                   if(nasym>0){
1093                       for(k=nzs[2]+nzs[0];k<nzs[2]+nzs[1];++k){
1094                           au[k]=aub[k]/dtime+au[k];
1095                       }
1096                   }
1097               }
1098           }
1099           
1100           if(*isolver==0){
1101 #ifdef SPOOLES
1102               if(*ithermal<2){
1103                   spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1104                           &symmetryflag,&inputformat,&nzs[2]);
1105               }else{
1106                   spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1107                           &symmetryflag,&inputformat,&nzs[2]);
1108               }
1109 #else
1110               printf(" *ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
1111               FORTRAN(stop,());
1112 #endif
1113           }
1114           else if((*isolver==2)||(*isolver==3)){
1115               if(nasym>0){
1116                   printf(" *ERROR in nonlingeo: the iterative solver cannot be used for asymmetric matrices\n\n");
1117                   FORTRAN(stop,());
1118               }
1119               preiter(ad,&au,b,&icol,&irow,&neq[1],&nzs[1],isolver,iperturb);
1120           }
1121           else if(*isolver==4){
1122 #ifdef SGI
1123               if(nasym>0){
1124                   printf(" *ERROR in nonlingeo: the SGI solver cannot be used for asymmetric matrices\n\n");
1125                   FORTRAN(stop,());
1126               }
1127               token=1;
1128               if(*ithermal<2){
1129                   sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
1130               }else{
1131                   sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],token);
1132               }
1133 #else
1134               printf(" *ERROR in nonlingeo: the SGI library is not linked\n\n");
1135               FORTRAN(stop,());
1136 #endif
1137           }
1138           else if(*isolver==5){
1139               if(nasym>0){
1140                   printf(" *ERROR in nonlingeo: the TAUCS solver cannot be used for asymmetric matrices\n\n");
1141                   FORTRAN(stop,());
1142               }
1143 #ifdef TAUCS
1144               tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[1],&nzs[1]);
1145 #else
1146               printf(" *ERROR in nonlingeo: the TAUCS library is not linked\n\n");
1147               FORTRAN(stop,());
1148 #endif
1149           }
1150           else if(*isolver==7){
1151 #ifdef PARDISO
1152               if(*ithermal<2){
1153                   pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1154                                &symmetryflag,&inputformat,jq,&nzs[2]);
1155               }else{
1156                   pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1157                                &symmetryflag,&inputformat,jq,&nzs[2]);
1158               }
1159 #else
1160               printf(" *ERROR in nonlingeo: the PARDISO library is not linked\n\n");
1161               FORTRAN(stop,());
1162 #endif
1163           }
1164           
1165           /* calculating the electromagnetic fields and temperatures 
1166              only the temperature calculation is differential */
1167           
1168           NNEW(v,double,mt**nk);
1169           memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1170           
1171           NNEW(fn,double,mt**nk);
1172           
1173           NNEW(inum,ITG,*nk);
1174           resultsinduction(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1175                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1176                   ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1177                   prestr,iprestr,filab,eme,emn,een,iperturb,
1178                   f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1179                   ndirboun,xbounact,nboun,ipompc,
1180                   nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1181                   &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1182                   xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1183                   &icmd,ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,
1184                   emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1185                   iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1186                   fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1187                   &reltime,ne,xforc,nforc,thicke,shcon,nshcon,
1188                   sideload,xloadact,xloadold,&icfd,inomat,h0,islavnode,
1189                   nslavnode,ntie,ielprop,prop,iactive);
1190           SFREE(inum);
1191           
1192           SFREE(ad);SFREE(au);
1193           
1194           if(*ithermal!=2){
1195               if(cam[0]>uam[0]){uam[0]=cam[0];}      
1196               if(qau<1.e-10){
1197                   if(qa[0]>ea*qam[0]){qam[0]=(qamold[0]*jnz+qa[0])/(jnz+1);}
1198                   else {qam[0]=qamold[0];}
1199               }
1200           }
1201           if(*ithermal>1){
1202               if(cam[1]>uam[1]){uam[1]=cam[1];}      
1203               if(qau<1.e-10){
1204                   if(qa[1]>ea*qam[1]){qam[1]=(qamold[1]*jnz+qa[1])/(jnz+1);}
1205                   else {qam[1]=qamold[1];}
1206               }
1207           }
1208           
1209           memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1210           
1211           SFREE(v);SFREE(fn);
1212           
1213           /* calculating the residual */
1214           
1215           calcresidual_em(nmethod,neq,b,fext,f,iexpl,nactdof,aux1,aux2,vold,
1216               vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1217               islavnode,nslavnode,&mortar,ntie,f_cm,f_cs,mi,
1218               nzs,&nasym,ithermal);
1219           
1220           /* calculating the maximum residual (only thermal part)*/
1221           
1222           for(k=0;k<2;++k){
1223               ram2[k]=ram1[k];
1224               ram1[k]=ram[k];
1225               ram[k]=0.;
1226           }
1227 
1228           if(*ithermal!=2) ram[0]=0.;
1229 
1230           if(*ithermal>1){
1231               for(k=neq[0];k<neq[1];++k){
1232                   err=fabs(b[k]);
1233                   if(err>ram[1]){ram[1]=err;ram[3]=k+0.5;}
1234               }
1235           }
1236           
1237           /* printing residuals */
1238           
1239           if(*ithermal>1){
1240               if(ram[1]<1.e-6) ram[1]=0.;      
1241               printf(" average flux= %f\n",qa[1]);
1242               printf(" time avg. flux= %f\n",qam[1]);
1243               if((ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1==0){
1244                   printf(" largest residual flux= %f\n",
1245                          ram[1]);
1246               }else{
1247                   inode=(ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1;
1248                   idir=nactdofinv[(ITG)ram[3]]-mt*(inode-1);
1249                   printf(" largest residual flux= %f in node %" ITGFORMAT " and dof %" ITGFORMAT "\n",ram[1],inode,idir);
1250               }
1251               printf(" largest increment of temp= %e\n",uam[1]);
1252               if((ITG)cam[4]==0){
1253                   printf(" largest correction to temp= %e\n\n",
1254                          cam[1]);
1255               }else{
1256                   inode=(ITG)((double)nactdofinv[(ITG)cam[4]]/mt)+1;
1257                   idir=nactdofinv[(ITG)cam[4]]-mt*(inode-1);
1258                   printf(" largest correction to temp= %e in node %" ITGFORMAT " and dof %" ITGFORMAT "\n\n",cam[1],inode,idir);
1259               }
1260           }
1261           
1262           /* athermal electromagnetic calculations are linear:
1263              set iit=2 to force convergence */
1264 
1265           if(*ithermal<=1) iit=2;
1266 
1267           checkconvergence(co,nk,kon,ipkon,lakon,ne,stn,nmethod, 
1268             kode,filab,een,t1act,&time,epn,ielmat,matname,enern, 
1269             xstaten,nstate_,istep,&iinc,iperturb,ener,mi,output,
1270             ithermal,qfn,&mode,&noddiam,trab,inotr,ntrans,orab,
1271             ielorien,norien,description,sti,&icutb,&iit,&dtime,qa,
1272             vold,qam,ram1,ram2,ram,cam,uam,&ntg,ttime,&icntrl,
1273             &theta,&dtheta,veold,vini,idrct,tper,&istab,tmax, 
1274             nactdof,b,tmin,ctrl,amta,namta,itpamp,&inext,&dthetaref,
1275             &itp,&jprint,jout,&uncoupled,t1,&iitterm,nelemload,
1276             nload,nodeboun,nboun,itg,ndirboun,&deltmx,&iflagact,
1277             set,nset,istartset,iendset,ialset,emn,thicke,jobnamec,
1278             &mortar,nmat,ielprop,prop);
1279       }
1280       
1281       /*********************************************************/
1282       /*   end of the iteration loop                          */
1283       /*********************************************************/
1284       
1285       /* icutb=0 means that the iterations in the increment converged,
1286          icutb!=0 indicates that the increment has to be reiterated with
1287          another increment size (dtheta) */
1288       
1289       if(((qa[0]>ea*qam[0])||(qa[1]>ea*qam[1]))&&(icutb==0)){jnz++;}
1290       iit=0;
1291       
1292       if(icutb!=0){
1293           memcpy(&vold[0],&vini[0],sizeof(double)*mt**nk);
1294           
1295           for(k=0;k<*nboun;++k){xbounact[k]=xbounini[k];}
1296           if((*ithermal==1)||(*ithermal>=3)){
1297               for(k=0;k<*nk;++k){t1act[k]=t1ini[k];}
1298           }
1299           for(k=0;k<neq[1];++k){
1300               f[k]=fini[k];
1301           }
1302           if(*nmethod==4){
1303               for(k=0;k<mt**nk;++k){
1304                   veold[k]=veini[k];
1305               }
1306               for(k=0;k<neq[1];++k){
1307                   fext[k]=fextini[k];
1308               }
1309           }
1310           
1311           qam[0]=qamold[0];
1312           qam[1]=qamold[1];
1313       }
1314       
1315       if((jout[0]==jprint)&&(icutb==0)){
1316           
1317           jprint=0;
1318           
1319           /* calculating the displacements and the stresses and storing */
1320           /* the results in frd format  */
1321           
1322           NNEW(v,double,mt**nk);
1323           NNEW(fn,double,mt**nk);
1324           if(*ithermal>1) NNEW(qfn,double,3**nk);
1325           if((strcmp1(&filab[3741],"EMFE")==0)||
1326              (strcmp1(&filab[3828],"EMFB")==0)) NNEW(stn,double,6**nk);
1327           NNEW(inum,ITG,*nk);
1328           
1329           memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1330           
1331           iout=2;
1332           icmd=3;
1333           
1334           resultsinduction(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1335                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1336                   ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1337                   prestr,iprestr,filab,eme,emn,een,iperturb,
1338                   f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1339                   ndirboun,xbounact,nboun,ipompc,
1340                   nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1341                   &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1342                   xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1343                   ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,emeini,
1344                   xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1345                   ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1346                   nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1347                   &reltime,ne,xforc,nforc,thicke,shcon,nshcon,
1348                   sideload,xloadact,xloadold,&icfd,inomat,h0,islavnode,
1349                   nslavnode,ntie,ielprop,prop,iactive);
1350           
1351           memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1352           
1353           iout=0;
1354           icmd=0;
1355           FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
1356           
1357           ++*kode;
1358           if(*mcs!=0){
1359               ptime=*ttime+time;
1360               frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
1361                  t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
1362                  nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
1363                  ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
1364                   norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,ne,
1365                  cdn,&mortar,nmat);
1366           }else{
1367               
1368               ptime=*ttime+time;
1369               frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
1370                   kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,
1371                   enern,xstaten,
1372                   nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1373                   ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1374                   mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1375                   cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1376                   thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1377               
1378           }
1379           
1380           SFREE(v);SFREE(fn);SFREE(inum);
1381           if(*ithermal>1){SFREE(qfn);}
1382           if(strcmp1(&filab[3741],"EMF ")==0) SFREE(stn);
1383           
1384       }
1385       
1386   }
1387   
1388   /*********************************************************/
1389   /*   end of the increment loop                          */
1390   /*********************************************************/
1391   
1392   if(jprint!=0){
1393       
1394       /* calculating the displacements and the stresses and storing  
1395          the results in frd format */
1396       
1397       NNEW(v,double,mt**nk);
1398       NNEW(fn,double,mt**nk);
1399       if(*ithermal>1) NNEW(qfn,double,3**nk);
1400       if(strcmp1(&filab[3741],"EMF ")==0) NNEW(stn,double,6**nk);
1401       NNEW(inum,ITG,*nk);
1402       
1403       memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1404       iout=2;
1405       icmd=3;
1406       
1407       resultsinduction(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1408               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1409               ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1410               prestr,iprestr,filab,eme,emn,een,iperturb,
1411               f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1412               ndirboun,xbounact,nboun,ipompc,
1413               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1414               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1415               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1416               ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,emeini,
1417               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1418               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1419               nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1420               &reltime,ne,xforc,nforc,thicke,shcon,nshcon,
1421               sideload,xloadact,xloadold,&icfd,inomat,h0,islavnode,
1422               nslavnode,ntie,ielprop,prop,iactive);
1423       
1424       memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1425       
1426       iout=0;
1427       icmd=0;
1428       FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
1429       
1430       ++*kode;
1431       if(*mcs>0){
1432           ptime=*ttime+time;
1433           frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
1434                  t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
1435                  nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
1436                  ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
1437                  norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,ne,
1438                  cdn,&mortar,nmat);
1439       }else{
1440           
1441           ptime=*ttime+time;
1442           frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
1443               kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1444               nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1445               ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1446               mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1447               cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1448               thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1449           
1450       }
1451       
1452       SFREE(v);SFREE(fn);SFREE(inum);
1453       if(*ithermal>1){SFREE(qfn);}
1454       if(strcmp1(&filab[3741],"EMF ")==0) SFREE(stn);
1455       
1456   }
1457 
1458       /* restoring the distributed loading  */
1459 
1460   if(*ithermal==3){
1461       *nload=nloadref;
1462       (*nload_)-=ne2;
1463       RENEW(nelemload,ITG,2**nload);memcpy(&nelemload[0],&nelemloadref[0],sizeof(ITG)*2**nload);
1464       if(*nam>0){
1465           RENEW(iamload,ITG,2**nload);
1466           memcpy(&iamload[0],&iamloadref[0],sizeof(ITG)*2**nload);
1467       }
1468       RENEW(idefload,ITG,*nload);memcpy(&idefload[0],&idefloadref[0],sizeof(ITG)**nload);
1469       RENEW(sideload,char,20**nload);memcpy(&sideload[0],&sideloadref[0],sizeof(char)*20**nload);
1470       
1471       /* freeing the temporary fields */
1472       
1473       SFREE(nelemloadref);if(*nam>0){SFREE(iamloadref);};SFREE(idefloadref);
1474       SFREE(sideloadref);
1475   }
1476   
1477   /* setting the velocity to zero at the end of a quasistatic or stationary
1478      step */
1479   
1480   if(abs(*nmethod)==1){
1481       for(k=0;k<mt**nk;++k){veold[k]=0.;}
1482   }
1483   
1484   /* updating the loading at the end of the step; 
1485      important in case the amplitude at the end of the step
1486      is not equal to one */
1487   
1488   for(k=0;k<*nboun;++k){
1489       
1490       /* thermal boundary conditions are updated only if the
1491          step was thermal or thermomechanical */
1492       
1493       if(ndirboun[k]==0){
1494           if(*ithermal<2) continue;
1495           
1496           /* mechanical boundary conditions are updated only
1497              if the step was not thermal or the node is a
1498              network node */
1499           
1500       }else if((ndirboun[k]>0)&&(ndirboun[k]<4)){
1501           node=nodeboun[k];
1502           FORTRAN(nident,(itg,&node,&ntg,&id));
1503           networknode=0;
1504           if(id>0){
1505               if(itg[id-1]==node) networknode=1;
1506           }
1507           if((*ithermal==2)&&(networknode==0)) continue;
1508       }
1509       xbounold[k]=xbounact[k];
1510   }
1511   for(k=0;k<*nforc;++k){xforcold[k]=xforcact[k];}
1512   for(k=0;k<2**nload;++k){xloadold[k]=xloadact[k];}
1513   for(k=0;k<7**nbody;k=k+7){xbodyold[k]=xbodyact[k];}
1514   if(*ithermal==1){
1515       for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
1516       for(k=0;k<*nk;++k){vold[mt*k]=t1act[k];}
1517   }
1518   else if(*ithermal>1){
1519       for(k=0;k<*nk;++k){t1[k]=vold[mt*k];}
1520       if(*ithermal>=3){
1521           for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
1522       }
1523   }
1524   
1525   qaold[0]=qa[0];
1526   qaold[1]=qa[1];
1527   
1528   SFREE(f);
1529   SFREE(b);
1530   SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(xbodyact);
1531   if(*nbody>0) SFREE(ipobody);
1532   SFREE(fext);SFREE(ampli);SFREE(xbounini);SFREE(xstiff);
1533   if((*ithermal==1)||(*ithermal>=3)){SFREE(t1act);SFREE(t1ini);}
1534   
1535   if(*ithermal>1){
1536       SFREE(itg);SFREE(ieg);SFREE(kontri);SFREE(nloadtr);
1537       SFREE(nactdog);SFREE(nacteq);SFREE(ineighe);
1538       SFREE(tarea);SFREE(tenv);SFREE(fenv);SFREE(qfx);
1539       SFREE(erad);SFREE(ac);SFREE(bc);SFREE(ipiv);
1540       SFREE(bcr);SFREE(ipivr);SFREE(adview);SFREE(auview);SFREE(adrad);
1541       SFREE(aurad);SFREE(irowrad);SFREE(jqrad);SFREE(icolrad);
1542       if((*mcs>0)&&(ntr>0)){SFREE(inocs);}
1543   }
1544   
1545   SFREE(fini);
1546   if(*nmethod==4){
1547       SFREE(aux2);SFREE(fextini);SFREE(veini);
1548       SFREE(adb);SFREE(aub);SFREE(cv);SFREE(cvini);
1549   }
1550   
1551   SFREE(aux);SFREE(iaux);SFREE(vini);SFREE(h0);SFREE(inomat);
1552   
1553   /* reset icascade */
1554   
1555   if(icascade==1){icascade=0;}
1556   
1557   mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
1558   mpcinfo[3]=maxlenmpc;
1559   
1560   *icolp=icol;*irowp=irow;*cop=co;*voldp=vold;
1561   
1562   *ipompcp=ipompc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
1563   *fmpcp=fmpc;*nodempcp=nodempc;*coefmpcp=coefmpc;
1564   
1565   *ipkonp=ipkon;*lakonp=lakon;*konp=kon;*ielorienp=ielorien;
1566   *ielmatp=ielmat;*enerp=ener;*xstatep=xstate;
1567 
1568   *setp=set;*istartsetp=istartset;*iendsetp=iendset;*ialsetp=ialset;
1569   *tiesetp=tieset;*tietolp=tietol;
1570   
1571   *nelemloadp=nelemload;*iamloadp=iamload;*idefloadp=idefload;
1572   *sideloadp=sideload;
1573   
1574   (*tmin)*=(*tper);
1575   (*tmax)*=(*tper);
1576   
1577   SFREE(nactdofinv);
1578  
1579   if(*nmethod==1){
1580       *nmethod=8;
1581   }else if(*nmethod==4){
1582       *nmethod=9;
1583   }else if(*nmethod==2){
1584       *nmethod=10;
1585   }
1586 
1587   (*ttime)+=(*tper);
1588   
1589   return;
1590 }

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