CalculiX  2.8
A Free Software Three-Dimensional Structural Finite Element Program
 All Classes Files Functions Variables Macros
nonlingeo.c File Reference
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "CalculiX.h"
Include dependency graph for nonlingeo.c:

Go to the source code of this file.

Functions

void nonlingeo (double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp, ITG *nmpc, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, ITG *nactdof, ITG **icolp, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG *ikboun, ITG *ilboun, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp, ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_, double *t0, double *t1, double *t1old, ITG *ithermal, double *prestr, ITG *iprestr, double **voldp, ITG *iperturb, double *sti, ITG *nzs, ITG *kode, char *filab, ITG *idrct, ITG *jmax, ITG *jout, double *timepar, double *eme, double *xbounold, double *xforcold, double *xloadold, double *veold, double *accold, char *amname, double *amta, ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload, ITG *iamt1, double *alpha, ITG *iexpl, ITG *iamboun, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double **xstatep, ITG *npmat_, ITG *istep, double *ttime, char *matname, double *qaold, ITG *mi, ITG *isolver, ITG *ncmat_, ITG *nstate_, ITG *iumat, double *cs, ITG *mcs, ITG *nkon, double **enerp, ITG *mpcinfo, char *output, double *shcon, ITG *nshcon, double *cocon, ITG *ncocon, double *physcon, ITG *nflow, double *ctrl, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, ITG *nener, ITG *ikforc, ITG *ilforc, double *trab, ITG *inotr, ITG *ntrans, double **fmpcp, char *cbody, ITG *ibody, double *xbody, ITG *nbody, double *xbodyold, ITG *ielprop, double *prop, ITG *ntie, char *tieset, ITG *itpamp, ITG *iviewfile, char *jobnamec, double *tietol, ITG *nslavs, double *thicke, ITG *ics, ITG *nintpoint, ITG *mortar, ITG *ifacecount, char *typeboun, ITG **islavsurfp, double **pslavsurfp, double **clearinip, ITG *nmat, double *xmodal, ITG *iaxial, ITG *inext)
 

Function Documentation

void nonlingeo ( double **  cop,
ITG nk,
ITG **  konp,
ITG **  ipkonp,
char **  lakonp,
ITG ne,
ITG nodeboun,
ITG ndirboun,
double *  xboun,
ITG nboun,
ITG **  ipompcp,
ITG **  nodempcp,
double **  coefmpcp,
char **  labmpcp,
ITG nmpc,
ITG nodeforc,
ITG ndirforc,
double *  xforc,
ITG nforc,
ITG nelemload,
char *  sideload,
double *  xload,
ITG nload,
ITG nactdof,
ITG **  icolp,
ITG jq,
ITG **  irowp,
ITG neq,
ITG nzl,
ITG nmethod,
ITG **  ikmpcp,
ITG **  ilmpcp,
ITG ikboun,
ITG ilboun,
double *  elcon,
ITG nelcon,
double *  rhcon,
ITG nrhcon,
double *  alcon,
ITG nalcon,
double *  alzero,
ITG **  ielmatp,
ITG **  ielorienp,
ITG norien,
double *  orab,
ITG ntmat_,
double *  t0,
double *  t1,
double *  t1old,
ITG ithermal,
double *  prestr,
ITG iprestr,
double **  voldp,
ITG iperturb,
double *  sti,
ITG nzs,
ITG kode,
char *  filab,
ITG idrct,
ITG jmax,
ITG jout,
double *  timepar,
double *  eme,
double *  xbounold,
double *  xforcold,
double *  xloadold,
double *  veold,
double *  accold,
char *  amname,
double *  amta,
ITG namta,
ITG nam,
ITG iamforc,
ITG iamload,
ITG iamt1,
double *  alpha,
ITG iexpl,
ITG iamboun,
double *  plicon,
ITG nplicon,
double *  plkcon,
ITG nplkcon,
double **  xstatep,
ITG npmat_,
ITG istep,
double *  ttime,
char *  matname,
double *  qaold,
ITG mi,
ITG isolver,
ITG ncmat_,
ITG nstate_,
ITG iumat,
double *  cs,
ITG mcs,
ITG nkon,
double **  enerp,
ITG mpcinfo,
char *  output,
double *  shcon,
ITG nshcon,
double *  cocon,
ITG ncocon,
double *  physcon,
ITG nflow,
double *  ctrl,
char *  set,
ITG nset,
ITG istartset,
ITG iendset,
ITG ialset,
ITG nprint,
char *  prlab,
char *  prset,
ITG nener,
ITG ikforc,
ITG ilforc,
double *  trab,
ITG inotr,
ITG ntrans,
double **  fmpcp,
char *  cbody,
ITG ibody,
double *  xbody,
ITG nbody,
double *  xbodyold,
ITG ielprop,
double *  prop,
ITG ntie,
char *  tieset,
ITG itpamp,
ITG iviewfile,
char *  jobnamec,
double *  tietol,
ITG nslavs,
double *  thicke,
ITG ics,
ITG nintpoint,
ITG mortar,
ITG ifacecount,
char *  typeboun,
ITG **  islavsurfp,
double **  pslavsurfp,
double **  clearinip,
ITG nmat,
double *  xmodal,
ITG iaxial,
ITG inext 
)

Definition at line 36 of file nonlingeo.c.

80  {
81 
82  char description[13]=" ",*lakon=NULL,jobnamef[396]="",
83  *sideface=NULL,*labmpc=NULL,fnffrd[132]="",*lakonf=NULL;
84 
85  ITG *inum=NULL,k,iout=0,icntrl,iinc=0,jprint=0,iit=-1,jnz=0,
86  icutb=0,istab=0,ifreebody,uncoupled,n1,n2,itruecontact,
87  iperturb_sav[2],ilin,*icol=NULL,*irow=NULL,ielas=0,icmd=0,
88  memmpc_,mpcfree,icascade,maxlenmpc,*nodempc=NULL,*iaux=NULL,
89  *nodempcref=NULL,memmpcref_,mpcfreeref,*itg=NULL,*ineighe=NULL,
90  *ieg=NULL,ntg=0,ntr,*kontri=NULL,*nloadtr=NULL,idamping=0,
91  *ipiv=NULL,ntri,newstep,mode=-1,noddiam=-1,nasym=0,im,
92  ntrit,*inocs=NULL,inewton=0,*ipobody=NULL,*nacteq=NULL,
93  *nactdog=NULL,nteq,network,*itietri=NULL,*koncont=NULL,
94  ncont,ne0,nkon0,*ipkon=NULL,*kon=NULL,*ielorien=NULL,
95  *ielmat=NULL,itp=0,symmetryflag=0,inputformat=0,
96  *iruc=NULL,iitterm=0,iturbulent,ngraph=1,ismallsliding=0,
97  *ipompc=NULL,*ikmpc=NULL,*ilmpc=NULL,i0ref,irref,icref,
98  *itiefac=NULL,*islavsurf=NULL,*islavnode=NULL,*imastnode=NULL,
99  *nslavnode=NULL,*nmastnode=NULL,*imastop=NULL,imat,
100  *iponoels=NULL,*inoels=NULL,*islavsurfold=NULL,maxlenmpcref,
101  *islavact=NULL,mt=mi[1]+1,*nactdofinv=NULL,*ipe=NULL,
102  *ime=NULL,*ikactmech=NULL,nactmech,inode,idir,neold,
103  iemchange=0,nzsrad,*mast1rad=NULL,*irowrad=NULL,*icolrad=NULL,
104  *jqrad=NULL,*ipointerrad=NULL,*integerglob=NULL,negpres=0,
105  mass[2]={0,0}, stiffness=1, buckling=0, rhsi=1, intscheme=0,idiscon=0,
106  coriolis=0,*ipneigh=NULL,*neigh=NULL,maxprevcontel,nslavs_prev_step,
107  *nelemface=NULL,*ipoface=NULL,*nodface=NULL,*ifreestream=NULL,iex,
108  *isolidsurf=NULL,*neighsolidsurf=NULL,*iponoel=NULL,*inoel=NULL,
109  nef=0,nface,nfreestream,nsolidsurf,i,indexe,icfd=0,id,*neij=NULL,
110  node,networknode,iflagact=0,*nodorig=NULL,*ipivr=NULL,iglob=0,
111  *inomat=NULL,*ipnei=NULL,ntrimax,*nx=NULL,*ny=NULL,*nz=NULL,
112  *neifa=NULL,*neiel=NULL,*ielfa=NULL,*ifaext=NULL,nflnei,nfaext,
113  idampingwithoutcontact=0,*nactdoh=NULL,*nactdohinv=NULL,*ipkonf=NULL,
114  *ielmatf=NULL,*ielorienf=NULL;
115 
116  double *stn=NULL,*v=NULL,*een=NULL,cam[5],*epn=NULL,*cg=NULL,
117  *cdn=NULL,*vel=NULL,*vfa=NULL,*pslavsurfold=NULL,
118  *f=NULL,*fn=NULL,qa[3]={0.,0.,-1.},qam[2]={0.,0.},dtheta,theta,
119  err,ram[8]={0.,0.,0.,0.,0.,0.,0.,0.},*areaslav=NULL,
120  *springarea=NULL,ram1[8]={0.,0.,0.,0.,0.,0.,0.,0.},
121  ram2[8]={0.,0.,0.,0.,0.,0.,0.,0.},deltmx,ptime,
122  uam[2]={0.,0.},*vini=NULL,*ac=NULL,qa0,qau,ea,*straight=NULL,
123  *t1act=NULL,qamold[2],*xbounact=NULL,*bc=NULL,
124  *xforcact=NULL,*xloadact=NULL,*fext=NULL,*clearini=NULL,
125  reltime,time,bet=0.,gam=0.,*aux2=NULL,dtime,*fini=NULL,
126  *fextini=NULL,*veini=NULL,*accini=NULL,*xstateini=NULL,
127  *ampli=NULL,scal1,*eei=NULL,*t1ini=NULL,pressureratio,
128  *xbounini=NULL,dev,*xstiff=NULL,*stx=NULL,*stiini=NULL,
129  *enern=NULL,*coefmpc=NULL,*aux=NULL,*xstaten=NULL,
130  *coefmpcref=NULL,*enerini=NULL,*emn=NULL,alpham,betam,
131  *tarea=NULL,*tenv=NULL,*erad=NULL,*fnr=NULL,*fni=NULL,
132  *adview=NULL,*auview=NULL,*qfx=NULL,*cvini=NULL,*cv=NULL,
133  *qfn=NULL,*co=NULL,*vold=NULL,*fenv=NULL,sigma=0.,
134  *xbodyact=NULL,*cgr=NULL,dthetaref, *vcontu=NULL,*vr=NULL,*vi=NULL,
135  *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,*fmpc=NULL,*ener=NULL,
136  *f_cm=NULL, *f_cs=NULL,*adc=NULL,*auc=NULL,
137  *xstate=NULL,*eenmax=NULL,*adrad=NULL,*aurad=NULL,*bcr=NULL,
138  *xmastnor=NULL,*emeini=NULL,*tinc,*tper,*tmin,*tmax,*tincf,
139  *doubleglob=NULL,*xnoels=NULL,*au=NULL,
140  *ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
141  *x=NULL,*y=NULL,*z=NULL,*xo=NULL,
142  *yo=NULL,*zo=NULL,*cdnr=NULL,*cdni=NULL;
143 
144 #ifdef SGI
145  ITG token;
146 #endif
147 
148  icol=*icolp;irow=*irowp;co=*cop;vold=*voldp;
149  ipkon=*ipkonp;lakon=*lakonp;kon=*konp;ielorien=*ielorienp;
150  ielmat=*ielmatp;ener=*enerp;xstate=*xstatep;
151 
152  ipompc=*ipompcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
153  fmpc=*fmpcp;nodempc=*nodempcp;coefmpc=*coefmpcp;
154 
155  islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
156 
157  tinc=&timepar[0];
158  tper=&timepar[1];
159  tmin=&timepar[2];
160  tmax=&timepar[3];
161  tincf=&timepar[4];
162 
163  if(*ithermal==4){
164  uncoupled=1;
165  *ithermal=3;
166  }else{
167  uncoupled=0;
168  }
169 
170  if(*mortar!=1){
171  maxprevcontel=*nslavs;
172  }else if(*mortar==1){
173  maxprevcontel=*nintpoint;
174  if(*nstate_!=0){
175  if(maxprevcontel!=0){
176  NNEW(islavsurfold,ITG,2**ifacecount+2);
177  NNEW(pslavsurfold,double,3**nintpoint);
178  memcpy(&islavsurfold[0],&islavsurf[0],
179  sizeof(ITG)*(2**ifacecount+2));
180  memcpy(&pslavsurfold[0],&pslavsurf[0],
181  sizeof(double)*(3**nintpoint));
182  }
183  }
184  }
185  nslavs_prev_step=*nslavs;
186 
187  /* turbulence model
188  iturbulent==0: laminar
189  iturbulent==1: k-epsilon
190  iturbulent==2: q-omega
191  iturbulent==3: SST */
192 
193  iturbulent=(ITG)physcon[8];
194 
195  for(k=0;k<3;k++){
196  strcpy1(&jobnamef[k*132],&jobnamec[k*132],132);
197  }
198 
199  qa0=ctrl[20];qau=ctrl[21];ea=ctrl[23];deltmx=ctrl[26];
200  i0ref=ctrl[0];irref=ctrl[1];icref=ctrl[3];
201 
202  memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
203  maxlenmpc=mpcinfo[3];
204 
205  alpham=xmodal[0];
206  betam=xmodal[1];
207 
208  /* check whether, for a dynamic calculation, damping is involved */
209 
210  if(*nmethod==4){
211  if(*iexpl<=1){
212  if((fabs(alpham)>1.e-30)||(fabs(betam)>1.e-30)){
213  idamping=1;
214  idampingwithoutcontact=1;
215  }else{
216  for(i=0;i<*ne;i++){
217  if(ipkon[i]<0) continue;
218  if(strcmp1(&lakon[i*8],"ED")==0){
219  idamping=1;idampingwithoutcontact=1;break;
220  }
221  }
222  }
223  }
224  }
225 
226  if((icascade==2)&&(*iexpl>=2)){
227  printf(" *ERROR in nonlingeo: linear and nonlinear MPC's depend on each other\n");
228  printf(" This is not allowed in a explicit dynamic calculation\n");
229  FORTRAN(stop,());
230  }
231 
232  /* check whether the submodel is meant for a fluid-structure
233  interaction */
234 
235  strcpy(fnffrd,jobnamec);
236  strcat(fnffrd,"f.frd");
237  if((jobnamec[396]!='\0')&&(strcmp1(&jobnamec[396],fnffrd)==0)){
238 
239  /* fluid-structure interaction: wait till after the compfluid
240  call */
241 
242  NNEW(integerglob,ITG,1);
243  NNEW(doubleglob,double,1);
244  }else{
245 
246  /* determining the global values to be used as boundary conditions
247  for a submodel */
248 
249  getglobalresults(jobnamec,&integerglob,&doubleglob,nboun,iamboun,xboun,
250  nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
251  ithermal,nk,t1,iamt1);
252  }
253 
254  /* invert nactdof */
255 
256  NNEW(nactdofinv,ITG,mt**nk);
257  NNEW(nodorig,ITG,*nk);
258  FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
259  ipkon,lakon,kon,ne));
260  SFREE(nodorig);
261 
262  /* allocating a field for the stiffness matrix */
263 
264  NNEW(xstiff,double,(long long)27*mi[0]**ne);
265 
266  /* allocating force fields */
267 
268  NNEW(f,double,neq[1]);
269  NNEW(fext,double,neq[1]);
270 
271  NNEW(b,double,neq[1]);
272  NNEW(vini,double,mt**nk);
273 
274  NNEW(aux,double,7*maxlenmpc);
275  NNEW(iaux,ITG,2*maxlenmpc);
276 
277  /* allocating fields for the actual external loading */
278 
279  NNEW(xbounact,double,*nboun);
280  NNEW(xbounini,double,*nboun);
281  for(k=0;k<*nboun;++k){xbounact[k]=xbounold[k];}
282  NNEW(xforcact,double,*nforc);
283  NNEW(xloadact,double,2**nload);
284  NNEW(xbodyact,double,7**nbody);
285  /* copying the rotation axis and/or acceleration vector */
286  for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
287 
288  /* assigning the body forces to the elements */
289 
290  if(*nbody>0){
291  ifreebody=*ne+1;
292  NNEW(ipobody,ITG,2*ifreebody**nbody);
293  for(k=1;k<=*nbody;k++){
294  FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
295  iendset,ialset,&inewton,nset,&ifreebody,&k));
296  RENEW(ipobody,ITG,2*(*ne+ifreebody));
297  }
298  RENEW(ipobody,ITG,2*(ifreebody-1));
299  if(inewton==1){NNEW(cgr,double,4**ne);}
300  }
301 
302  /* for mechanical calculations: updating boundary conditions
303  calculated in a previous thermal step */
304 
305  if(*ithermal<2) FORTRAN(gasmechbc,(vold,nload,sideload,
306  nelemload,xload,mi));
307 
308  /* for thermal calculations: forced convection and cavity
309  radiation*/
310 
311  if(*ithermal>1){
312  NNEW(itg,ITG,*nload+3**nflow);
313  NNEW(ieg,ITG,*nflow);
314  /* max 6 triangles per face, 4 entries per triangle */
315  NNEW(kontri,ITG,24**nload);
316  NNEW(nloadtr,ITG,*nload);
317  NNEW(nacteq,ITG,4**nk);
318  NNEW(nactdog,ITG,4**nk);
319  NNEW(v,double,mt**nk);
320  FORTRAN(envtemp,(itg,ieg,&ntg,&ntr,sideload,nelemload,
321  ipkon,kon,lakon,ielmat,ne,nload,
322  kontri,&ntri,nloadtr,nflow,ndirboun,nactdog,
323  nodeboun,nacteq,nboun,ielprop,prop,&nteq,
324  v,&network,physcon,shcon,ntmat_,co,
325  vold,set,nshcon,rhcon,nrhcon,mi,nmpc,nodempc,
326  ipompc,labmpc,ikboun,&nasym,iaxial));
327  SFREE(v);
328 
329  if((*mcs>0)&&(ntr>0)){
330  NNEW(inocs,ITG,*nk);
331  radcyc(nk,kon,ipkon,lakon,ne,cs,mcs,nkon,ialset,istartset,
332  iendset,&kontri,&ntri,&co,&vold,&ntrit,inocs,mi);
333  }
334  else{ntrit=ntri;}
335 
336  nzsrad=100*ntr;
337  NNEW(mast1rad,ITG,nzsrad);
338  NNEW(irowrad,ITG,nzsrad);
339  NNEW(icolrad,ITG,ntr);
340  NNEW(jqrad,ITG,ntr+1);
341  NNEW(ipointerrad,ITG,ntr);
342 
343  if(ntr>0){
344  mastructrad(&ntr,nloadtr,sideload,ipointerrad,
345  &mast1rad,&irowrad,&nzsrad,
346  jqrad,icolrad);
347  }
348 
349  SFREE(ipointerrad);SFREE(mast1rad);
350  RENEW(irowrad,ITG,nzsrad);
351 
352  RENEW(itg,ITG,ntg);
353  NNEW(ineighe,ITG,ntg);
354  RENEW(kontri,ITG,4*ntrit);
355  RENEW(nloadtr,ITG,ntr);
356 
357  NNEW(adview,double,ntr);
358  NNEW(auview,double,2*nzsrad);
359  NNEW(tarea,double,ntr);
360  NNEW(tenv,double,ntr);
361  NNEW(fenv,double,ntr);
362  NNEW(erad,double,ntr);
363 
364  NNEW(ac,double,nteq*nteq);
365  NNEW(bc,double,nteq);
366  NNEW(ipiv,ITG,nteq);
367  NNEW(adrad,double,ntr);
368  NNEW(aurad,double,2*nzsrad);
369  NNEW(bcr,double,ntr);
370  NNEW(ipivr,ITG,ntr);
371  }
372 
373  /* check for fluid elements */
374 
375  NNEW(nactdoh,ITG,*ne);
376  NNEW(nactdohinv,ITG,*ne);
377  for(i=0;i<*ne;++i){
378  if(ipkon[i]<0) continue;
379  indexe=ipkon[i];
380  if(strcmp1(&lakon[8*i],"F")==0){
381  icfd=1;nactdohinv[nef]=i+1;nef++;nactdoh[i]=nef;}
382  }
383  if(icfd==1){
384  RENEW(nactdohinv,ITG,nef);
385  NNEW(ipkonf,ITG,nef);
386  NNEW(lakonf,char,8*nef);
387  NNEW(ielmatf,ITG,mi[2]*nef);
388  if(*norien>0) NNEW(ielorienf,ITG,mi[2]*nef);
389  NNEW(ipoface,ITG,*nk);
390  NNEW(nodface,ITG,5*6*nef);
391  NNEW(ipnei,ITG,*ne);
392  NNEW(neifa,ITG,6**ne);
393  NNEW(neiel,ITG,6**ne);
394  NNEW(neij,ITG,6**ne);
395  NNEW(ielfa,ITG,24**ne);
396  NNEW(ifaext,ITG,6**ne);
397  NNEW(isolidsurf,ITG,6**ne);
398  NNEW(vel,double,6**ne);
399  NNEW(vfa,double,6*6**ne);
400 
401  FORTRAN(precfd,(ne,ipkon,kon,lakon,ipnei,neifa,neiel,ipoface,
402  nodface,ielfa,&nflnei,&nface,ifaext,&nfaext,
403  isolidsurf,&nsolidsurf,set,nset,istartset,iendset,ialset,
404  vel,vold,mi,neij,&nef,nactdoh,ipkonf,lakonf,ielmatf,ielmat,
405  ielorienf,ielorien,norien));
406 
407  SFREE(ipoface);
408  SFREE(nodface);
409  RENEW(neifa,ITG,nflnei);
410  RENEW(neiel,ITG,nflnei);
411  RENEW(neij,ITG,nflnei);
412  RENEW(ielfa,ITG,4*nface);
413  RENEW(ifaext,ITG,nfaext);
414  RENEW(isolidsurf,ITG,nsolidsurf);
415  RENEW(vfa,double,6*nface);
416  }else{
417  SFREE(nactdoh);SFREE(nactdohinv);
418  }
419  if(*ithermal>1){NNEW(qfx,double,3*mi[0]**ne);}
420 
421  /* contact conditions */
422 
423  inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
424  lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
425  &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
426  mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
427  nmpc,&mpcfree,&memmpc_,
428  &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
429  iperturb,ikboun,nboun,co,istep,&xnoels);
430 
431  if(ncont!=0){
432 
433  NNEW(cg,double,3*ncont);
434  NNEW(straight,double,16*ncont);
435 
436  /* 11 instead of 10: last position is reserved for the
437  local contact spring element number; needed as
438  pointer into springarea */
439 
440  if(*mortar==0){
441  RENEW(kon,ITG,*nkon+11**nslavs);
442  NNEW(springarea,double,2**nslavs);
443  if(*nener==1){
444  RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
445  }
446  RENEW(ipkon,ITG,*ne+*nslavs);
447  RENEW(lakon,char,8*(*ne+*nslavs));
448 
449  if(*norien>0){
450  RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
451  for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
452  }
453 
454  RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
455  for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
456 
457  if((maxprevcontel==0)&&(*nslavs!=0)){
458  RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
459  for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
460  xstate[k]=0.;
461  }
462  }
463  maxprevcontel=*nslavs;
464 
465  NNEW(areaslav,double,*ifacecount);
466  }else if(*mortar==1){
467  NNEW(islavact,ITG,nslavnode[*ntie]);
468  if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
469 
470  /* check whether at least one contact definition involves true contact
471  and not just tied contact */
472 
473  FORTRAN(checktruecontact,(ntie,tieset,tietol,elcon,&itruecontact,
474  ncmat_,ntmat_));
475  }
476 
477  NNEW(xmastnor,double,3*nmastnode[*ntie]);
478  }
479 
480  if((icascade==2)||(ncont!=0)){
481  memmpcref_=memmpc_;mpcfreeref=mpcfree;maxlenmpcref=maxlenmpc;
482  NNEW(nodempcref,ITG,3*memmpc_);
483  for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
484  NNEW(coefmpcref,double,memmpc_);
485  for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
486  }
487 
488  if((*ithermal==1)||(*ithermal>=3)){
489  NNEW(t1ini,double,*nk);
490  NNEW(t1act,double,*nk);
491  for(k=0;k<*nk;++k){t1act[k]=t1old[k];}
492  }
493 
494  /* allocating a field for the instantaneous amplitude */
495 
496  NNEW(ampli,double,*nam);
497 
498  /* fini is also needed in static calculations if ilin=1
499  to get correct values of f after a divergent increment */
500 
501  NNEW(fini,double,neq[1]);
502 
503  /* allocating fields for nonlinear dynamics */
504 
505  if(*nmethod==4){
506  mass[0]=1;
507  mass[1]=1;
508  NNEW(aux2,double,neq[1]);
509 // NNEW(fini,double,neq[1]);
510  NNEW(fextini,double,neq[1]);
511  NNEW(veini,double,mt**nk);
512  NNEW(accini,double,mt**nk);
513  NNEW(adb,double,neq[1]);
514  NNEW(aub,double,nzs[1]);
515  NNEW(cvini,double,neq[1]);
516  NNEW(cv,double,neq[1]);
517  }
518 
519  if((*nstate_!=0)&&((*mortar!=1)||(ncont==0))){
520  NNEW(xstateini,double,*nstate_*mi[0]*(*ne+*nslavs));
521  for(k=0;k<*nstate_*mi[0]*(*ne+*nslavs);++k){
522  xstateini[k]=xstate[k];
523  }
524  }
525  if((*nstate_!=0)&&(*mortar==1)) NNEW(xstateini,double,1);
526  NNEW(eei,double,6*mi[0]**ne);
527  NNEW(stiini,double,6*mi[0]**ne);
528  NNEW(emeini,double,6*mi[0]**ne);
529  if(*nener==1)
530  NNEW(enerini,double,mi[0]**ne);
531 
532  qa[0]=qaold[0];
533  qa[1]=qaold[1];
534 
535  /* normalizing the time */
536 
537  FORTRAN(checktime,(itpamp,namta,tinc,ttime,amta,tmin,inext,&itp,istep));
538  dtheta=(*tinc)/(*tper);
539 
540  /* taking care of a small increment at the end of the step
541  for face-to-face penalty contact */
542 
543  dthetaref=dtheta;
544  if((dtheta<=1.e-6)&&(*iexpl<=1)){
545  printf("\n *ERROR in nonlingeo\n");
546  printf(" increment size smaller than one millionth of step size\n");
547  printf(" increase increment size\n\n");
548  }
549  *tmin=*tmin/(*tper);
550  *tmax=*tmax/(*tper);
551  theta=0.;
552 
553  /* calculating an initial flux norm */
554 
555  if(*ithermal!=2){
556  if(qau>1.e-10){qam[0]=qau;}
557  else if(qa0>1.e-10){qam[0]=qa0;}
558  else if(qa[0]>1.e-10){qam[0]=qa[0];}
559  else {qam[0]=1.e-2;}
560  }
561  if(*ithermal>1){
562  if(qau>1.e-10){qam[1]=qau;}
563  else if(qa0>1.e-10){qam[1]=qa0;}
564  else if(qa[1]>1.e-10){qam[1]=qa[1];}
565  else {qam[1]=1.e-2;}
566  }
567 
568  /* storing the element and topology information before introducing
569  contact elements */
570 
571  ne0=*ne;nkon0=*nkon;neold=*ne;
572 
573  /*********************************************************************/
574 
575  /* calculating of the acceleration due to force discontinuities
576  (external - internal force) at the start of a step */
577 
578  /*********************************************************************/
579 
580  if((*nmethod==4)&&(*ithermal!=2)){
581  bet=(1.-*alpha)*(1.-*alpha)/4.;
582  gam=0.5-*alpha;
583 
584  /* calculating the stiffness and mass matrix */
585 
586  reltime=0.;
587  time=0.;
588  dtime=0.;
589 
590  FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
591  xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
592  t1old,t1,t1act,iamt1,nk,amta,
593  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
594  xbounold,xboun,xbounact,iamboun,nboun,
595  nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
596  co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
597  ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
598  iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
599 
600  time=0.;
601  dtime=1.;
602 
603  /* updating the nonlinear mpc's (also affects the boundary
604  conditions through the nonhomogeneous part of the mpc's)
605  if contact arises the number of MPC's can also change */
606 
607  cam[0]=0.;cam[1]=0.;cam[2]=0.;
608 
609  if(icascade==2){
610  memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
611  RENEW(nodempc,ITG,3*memmpcref_);
612  for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
613  RENEW(coefmpc,double,memmpcref_);
614  for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
615  }
616 
617  newstep=0;
618  FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
619  nmpc,ikboun,ilboun,nboun,xbounold,aux,iaux,
620  &maxlenmpc,ikmpc,ilmpc,&icascade,
621  kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
622  &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
623  if(icascade==2){
624 // memmpcref_=memmpc_;mpcfreeref=mpcfree;
625 // RENEW(nodempcref,ITG,3*memmpc_);
626  for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
627 // RENEW(coefmpcref,double,memmpc_);
628  for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
629  }
630 
631  if(icascade>0) remastruct(ipompc,&coefmpc,&nodempc,nmpc,
632  &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
633  labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
634  kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
635  neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
636  &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
637  typeboun,&cv,&cvini,&iit);
638 
639  /* invert nactdof */
640 
641  SFREE(nactdofinv);
642  NNEW(nactdofinv,ITG,1);
643 
644  iout=-1;
645  ielas=1;
646 
647  NNEW(fn,double,mt**nk);
648  NNEW(stx,double,6*mi[0]**ne);
649 
650  NNEW(inum,ITG,*nk);
651  results(co,nk,kon,ipkon,lakon,ne,vold,stn,inum,stx,
652  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
653  ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
654  prestr,iprestr,filab,eme,emn,een,iperturb,
655  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
656  ndirboun,xbounold,nboun,ipompc,
657  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,&bet,
658  &gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
659  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
660  ncmat_,nstate_,sti,vini,ikboun,ilboun,ener,enern,emeini,xstaten,
661  eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
662  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
663  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
664  &ne0,xforc,nforc,thicke,shcon,nshcon,
665  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
666  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
667  islavsurf,ielprop,prop);
668 
669  SFREE(fn);SFREE(stx);SFREE(inum);
670 
671  iout=0;
672  ielas=0;
673 
674  reltime=0.;
675  time=0.;
676  dtime=0.;
677 
678  if(*iexpl<=1){intscheme=1;}
679 
680  /* in mafillsm the stiffness and mass matrix are computed;
681  The primary aim is to calculate the mass matrix (not
682  lumped for an implicit dynamic calculation, lumped for an
683  explicit dynamic calculation). However:
684  - for an implicit calculation the mass matrix is "doped" with
685  a small amount of stiffness matrix, therefore the calculation
686  of the stiffness matrix is needed.
687  - for an explicit calculation the stiffness matrix is not
688  needed at all. Since the calculation of the mass matrix alone
689  is not possible in mafillsm, the determination of the stiffness
690  matrix is taken as unavoidable "ballast". */
691 
692  NNEW(ad,double,neq[1]);
693  NNEW(au,double,nzs[1]);
694 
695  FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounold,nboun,
696  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
697  nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
698  nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
699  nmethod,ikmpc,ilmpc,ikboun,ilboun,
700  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
701  ielmat,ielorien,norien,orab,ntmat_,
702  t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
703  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
704  xstiff,npmat_,&dtime,matname,mi,
705  ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
706  physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
707  &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
708  xstateini,xstate,thicke,integerglob,doubleglob,
709  tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
710  pmastsurf,mortar,clearini,ielprop,prop));
711 
712  if(*nmethod==0){
713 
714  /* error occurred in mafill: storing the geometry in frd format */
715 
716  ++*kode;
717  if(strcmp1(&filab[1044],"ZZS")==0){
718  NNEW(neigh,ITG,40**ne);
719  NNEW(ipneigh,ITG,*nk);
720  }
721 
722  ptime=*ttime+time;
723  frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
724  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
725  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
726  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
727  mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
728  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
729  thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
730 
731  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
732 #ifdef COMPANY
733  FORTRAN(uout,(v,mi,ithermal));
734 #endif
735  FORTRAN(stop,());
736 
737  }
738 
739  /* mass x acceleration = f(external)-f(internal)
740  only for the mechanical loading*/
741 
742  for(k=0;k<neq[0];++k){
743  b[k]=fext[k]-f[k];
744  }
745 
746  if(*iexpl<=1){
747 
748  /* a small amount of stiffness is added to the mass matrix
749  otherwise the system leads to huge accelerations in
750  case of discontinuous load changes at the start of the step */
751 
752  dtime=*tinc/10.;
753  scal1=bet*dtime*dtime*(1.+*alpha);
754  for(k=0;k<neq[0];++k){
755  ad[k]=adb[k]+scal1*ad[k];
756  }
757  for(k=0;k<nzs[0];++k){
758  au[k]=aub[k]+scal1*au[k];
759  }
760  if(*isolver==0){
761 #ifdef SPOOLES
762  spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
763  &symmetryflag,&inputformat,&nzs[2]);
764 #else
765  printf(" *ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
766  FORTRAN(stop,());
767 #endif
768  }
769  else if((*isolver==2)||(*isolver==3)){
770  preiter(ad,&au,b,&icol,&irow,&neq[0],&nzs[0],isolver,iperturb);
771  }
772  else if(*isolver==4){
773 #ifdef SGI
774  token=1;
775  sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
776 #else
777  printf(" *ERROR in nonlingeo: the SGI library is not linked\n\n");
778  FORTRAN(stop,());
779 #endif
780  }
781  else if(*isolver==5){
782 #ifdef TAUCS
783  tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[0],&nzs[0]);
784 #else
785  printf(" *ERROR in nonlingeo: the TAUCS library is not linked\n\n");
786  FORTRAN(stop,());
787 #endif
788  }
789  else if(*isolver==7){
790 #ifdef PARDISO
791  pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
792  &symmetryflag,&inputformat,jq,&nzs[2]);
793 #else
794  printf(" *ERROR in nonlingeo: the PARDISO library is not linked\n\n");
795  FORTRAN(stop,());
796 #endif
797  }
798  }
799 
800  else{
801  for(k=0;k<neq[0];++k){
802  b[k]=(fext[k]-f[k])/adb[k];
803  }
804  }
805 
806  /* for thermal loading the acceleration is set to zero */
807 
808  for(k=neq[0];k<neq[1];++k){
809  b[k]=0.;
810  }
811 
812  /* calculating the displacements, stresses and forces */
813 
814  NNEW(v,double,mt**nk);
815  memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
816 
817  NNEW(stx,double,6*mi[0]**ne);
818  NNEW(fn,double,mt**nk);
819 
820  NNEW(inum,ITG,*nk);
821  dtime=1.e-20;
822  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
823  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
824  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
825  prestr,iprestr,filab,eme,emn,een,iperturb,
826  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
827  ndirboun,xbounact,nboun,ipompc,
828  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
829  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
830  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
831  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
832  emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
833  iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
834  fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
835  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
836  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
837  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
838  islavsurf,ielprop,prop);
839  SFREE(inum);
840  dtime=0.;
841 
842  memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
843  if(*ithermal!=2){
844  for(k=0;k<6*mi[0]*ne0;++k){
845  sti[k]=stx[k];
846  }
847  }
848 
849  SFREE(v);SFREE(stx);SFREE(fn);
850  SFREE(ad);SFREE(au);
851 
852  /* the mass matrix is kept for subsequent calculations, therefore,
853  no new mass calculation is necessary for the remaining iterations
854  in the present step */
855 
856  mass[0]=0;intscheme=0;
857 /* DMEMSET(f,0,neq[1],0.);
858  DMEMSET(fext,0,neq[1],0.);*/
859 
860  }
861 
862  if(*iexpl>1) icmd=3;
863 
864  /**************************************************************/
865  /* starting the loop over the increments */
866  /**************************************************************/
867 
868  newstep=1;
869 
870  while((1.-theta>1.e-6)||(negpres==1)){
871 
872  if(icutb==0){
873 
874  /* previous increment converged: update the initial values */
875 
876  iinc++;
877  jprint++;
878 
879  /* vold is copied into vini */
880 
881  memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
882 
883  for(k=0;k<*nboun;++k){xbounini[k]=xbounact[k];}
884  if((*ithermal==1)||(*ithermal>=3)){
885  for(k=0;k<*nk;++k){t1ini[k]=t1act[k];}
886  }
887  for(k=0;k<neq[1];++k){
888  fini[k]=f[k];
889  }
890  if(*nmethod==4){
891  for(k=0;k<mt**nk;++k){
892  veini[k]=veold[k];
893  accini[k]=accold[k];
894  }
895  for(k=0;k<neq[1];++k){
896 // fini[k]=f[k];
897  fextini[k]=fext[k];
898  cvini[k]=cv[k];
899  }
900  }
901  if(*ithermal!=2){
902  for(k=0;k<6*mi[0]*ne0;++k){
903  stiini[k]=sti[k];
904  emeini[k]=eme[k];
905  }
906  }
907  if(*nener==1)
908  for(k=0;k<mi[0]*ne0;++k){enerini[k]=ener[k];}
909 
910  if(*mortar!=1){
911  if(*nstate_!=0){
912  for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
913  xstateini[k]=xstate[k];
914  }
915  }
916  }
917  }
918 
919  /* check for max. # of increments */
920 
921  if(iinc>*jmax){
922  printf(" *ERROR: max. # of increments reached\n\n");
923  FORTRAN(stop,());
924  }
925  printf(" increment %" ITGFORMAT " attempt %" ITGFORMAT " \n",iinc,icutb+1);
926  printf(" increment size= %e\n",dtheta**tper);
927  printf(" sum of previous increments=%e\n",theta**tper);
928  printf(" actual step time=%e\n",(theta+dtheta)**tper);
929  printf(" actual total time=%e\n\n",*ttime+(theta+dtheta)**tper);
930 
931  printf(" iteration 1\n\n");
932 
933  qamold[0]=qam[0];
934  qamold[1]=qam[1];
935 
936  /* determining the actual loads at the end of the new increment*/
937 
938  reltime=theta+dtheta;
939  time=reltime**tper;
940  dtime=dtheta**tper;
941 
942  FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
943  xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
944  t1old,t1,t1act,iamt1,nk,amta,
945  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
946  xbounold,xboun,xbounact,iamboun,nboun,
947  nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
948  co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
949  ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
950  iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
951 
952  for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
953  if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,bcr,ipivr,
954  ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,ntmat_,vold,
955  shcon,nshcon,ipkon,kon,co,
956  kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
957  nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
958  cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
959  ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
960  nactdog,nacteq,nodeboun,ndirboun,&network,
961  rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
962  ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,nset,
963  ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,iamload,
964  jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
965  }
966 
967  if(icfd==1){
968  compfluid(&co,nk,&ipkonf,&kon,&lakonf,&sideface,
969  ifreestream,&nfreestream,isolidsurf,neighsolidsurf,&nsolidsurf,
970  &iponoel,&inoel,nshcon,shcon,nrhcon,rhcon,&vold,ntmat_,nodeboun,
971  ndirboun,nboun,&ipompc,&nodempc,nmpc,&ikmpc,&ilmpc,ithermal,
972  ikboun,ilboun,&iturbulent,isolver,iexpl,vcontu,ttime,
973  &time,&dtime,nodeforc,ndirforc,xforc,nforc,nelemload,sideload,
974  xload,nload,xbody,ipobody,nbody,ielmatf,matname,mi,ncmat_,
975  physcon,istep,&iinc,ibody,xloadold,xboun,&coefmpc,
976  nmethod,xforcold,xforcact,iamforc,iamload,xbodyold,xbodyact,
977  t1old,t1,t1act,iamt1,amta,namta,nam,ampli,xbounold,xbounact,
978  iamboun,itg,&ntg,amname,t0,&nelemface,&nface,cocon,ncocon,xloadact,
979  tper,jmax,jout,set,nset,istartset,iendset,ialset,prset,prlab,
980  nprint,trab,inotr,ntrans,filab,&labmpc,sti,norien,orab,jobnamef,
981  tieset,ntie,mcs,ics,cs,nkon,&mpcfree,&memmpc_,&fmpc,&nef,&inomat,
982  qfx,neifa,neiel,ielfa,ifaext,vfa,vel,ipnei,&nflnei,&nfaext,
983  typeboun,neij,tincf,nactdoh,nactdohinv,ielorienf);
984 
985  /* determining the global values to be used as boundary conditions
986  for a submodel */
987 
988  SFREE(integerglob);SFREE(doubleglob);
989  getglobalresults(jobnamec,&integerglob,&doubleglob,nboun,iamboun,
990  xboun,nload,sideload,iamload,&iglob,nforc,iamforc,
991  xforc,ithermal,nk,t1,iamt1);
992  }
993 
994  if((icascade==2)||(ncont!=0)){
995  memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
996  RENEW(nodempc,ITG,3*memmpcref_);
997  for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
998  RENEW(coefmpc,double,memmpcref_);
999  for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
1000  }
1001 
1002  /* generating contact elements */
1003 
1004  if((ncont!=0)&&(*mortar<=1)){
1005  *ne=ne0;*nkon=nkon0;
1006 
1007  /* at start of new increment:
1008  - copy state variables (node-to-face)
1009  - determine slave integration points (face-to-face)
1010  - interpolate state variables (face-to-face) */
1011 
1012  if(icutb==0){
1013  if(*mortar==1){
1014 
1015  if(*nstate_!=0){
1016  if(maxprevcontel!=0){
1017  if(iit!=-1){
1018  NNEW(islavsurfold,ITG,2**ifacecount+2);
1019  NNEW(pslavsurfold,double,3**nintpoint);
1020  memcpy(&islavsurfold[0],&islavsurf[0],
1021  sizeof(ITG)*(2**ifacecount+2));
1022  memcpy(&pslavsurfold[0],&pslavsurf[0],
1023  sizeof(double)*(3**nintpoint));
1024  }
1025  }
1026  }
1027 
1028  *nintpoint=0;
1029 
1030  precontact(&ncont,ntie,tieset,nset,set,istartset,
1031  iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
1032  cg,straight,co,vold,istep,&iinc,&iit,itiefac,
1033  islavsurf,islavnode,imastnode,nslavnode,nmastnode,
1034  imastop,mi,ipe,ime,tietol,&iflagact,
1035  nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
1036  nslavs);
1037 
1038  /* changing the dimension of element-related fields */
1039 
1040  RENEW(kon,ITG,*nkon+22**nintpoint);
1041  RENEW(springarea,double,2**nintpoint);
1042  RENEW(pmastsurf,double,6**nintpoint);
1043 
1044  if(*nener==1){
1045  RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
1046  }
1047  RENEW(ipkon,ITG,*ne+*nintpoint);
1048  RENEW(lakon,char,8*(*ne+*nintpoint));
1049 
1050  if(*norien>0){
1051  RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
1052  for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
1053  }
1054  RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
1055  for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
1056 
1057  }
1058  }
1059 
1060  contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1061  ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
1062  co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
1063  &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
1064  &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
1065  iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,
1066  islavsurf,
1067  itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
1068  imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
1069  xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta);
1070 
1071  /* check whether, for a dynamic calculation, damping is involved */
1072 
1073  if(*nmethod==4){
1074  if(*iexpl<=1){
1075  if(idampingwithoutcontact==0){
1076  for(i=0;i<*ne;i++){
1077  if(ipkon[i]<0) continue;
1078  if(*ncmat_>=5){
1079  if(strcmp1(&lakon[i*8],"ES")==0){
1080  if(strcmp1(&lakon[i*8+6],"C")==0){
1081  imat=ielmat[i*mi[2]];
1082  if(elcon[(*ncmat_+1)**ntmat_*(imat-1)+4]>0.){
1083  idamping=1;break;
1084  }
1085  }
1086  }
1087  }
1088  }
1089  }
1090  }
1091  }
1092 
1093  printf(" Number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
1094 
1095  /* interpolating the state variables */
1096 
1097  if(icutb==0){
1098  if(*mortar==1){
1099  if(*nstate_!=0){
1100  if(maxprevcontel!=0){
1101  RENEW(xstateini,double,
1102  *nstate_*mi[0]*(ne0+maxprevcontel));
1103  for(k=*nstate_*mi[0]*ne0;
1104  k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
1105  xstateini[k]=xstate[k];
1106  }
1107  }
1108 
1109  RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
1110  for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
1111  xstate[k]=0.;
1112  }
1113 
1114  if((*nintpoint>0)&&(maxprevcontel>0)){
1115  iex=2;
1116 
1117  /* interpolation of xstate */
1118 
1119  FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
1120  &ne0,mi,xstate,pslavsurf,nstate_,
1121  xstateini,islavsurf,islavsurfold,
1122  pslavsurfold));
1123 
1124  }
1125 
1126  if(maxprevcontel!=0){
1127  SFREE(islavsurfold);SFREE(pslavsurfold);
1128  }
1129 
1130  maxprevcontel=*nintpoint;
1131 
1132  RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
1133  for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
1134  xstateini[k]=xstate[k];
1135  }
1136  }
1137  }
1138  }
1139 
1140  if(icascade<1)icascade=1;
1141  }
1142 
1143  /* updating the nonlinear mpc's (also affects the boundary
1144  conditions through the nonhomogeneous part of the mpc's) */
1145 
1146  FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1147  nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1148  &maxlenmpc,ikmpc,ilmpc,&icascade,
1149  kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1150  &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1151 
1152  if((icascade==2)||(ncont!=0)){
1153 // memmpcref_=memmpc_;mpcfreeref=mpcfree;
1154 // RENEW(nodempcref,ITG,3*memmpc_);
1155  for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
1156 // RENEW(coefmpcref,double,memmpc_);
1157  for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
1158  }
1159 
1160  if(icascade>0) remastruct(ipompc,&coefmpc,&nodempc,nmpc,
1161  &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
1162  labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
1163  kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
1164  neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
1165  &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
1166  typeboun,&cv,&cvini,&iit);
1167 
1168  /* invert nactdof */
1169 
1170  SFREE(nactdofinv);
1171  NNEW(nactdofinv,ITG,mt**nk);
1172  NNEW(nodorig,ITG,*nk);
1173  FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
1174  ipkon,lakon,kon,ne));
1175  SFREE(nodorig);
1176 
1177  /* check whether the forced displacements changed; if so, and
1178  if the procedure is static, the first iteration has to be
1179  purely linear elastic, in order to get an equilibrium
1180  displacement field; otherwise huge (maybe nonelastic)
1181  stresses may occur, jeopardizing convergence */
1182 
1183  ilin=0;
1184 
1185  /* only for iinc=1 a linearized calculation is performed, since
1186  for iinc>1 a reasonable displacement field is predicted by using the
1187  initial velocity field at the end of the last increment */
1188 
1189  if((iinc==1)&&(*ithermal<2)){
1190  dev=0.;
1191  for(k=0;k<*nboun;++k){
1192  err=fabs(xbounact[k]-xbounini[k]);
1193  if(err>dev){dev=err;}
1194  }
1195  if(dev>1.e-5) ilin=1;
1196  }
1197 
1198  /* prediction of the kinematic vectors */
1199 
1200  NNEW(v,double,mt**nk);
1201 
1202  if(ncont>idiscon)idiscon=ncont;
1203  prediction(uam,nmethod,&bet,&gam,&dtime,ithermal,nk,veold,accold,v,
1204  &iinc,&idiscon,vold,nactdof,mi);
1205 
1206  NNEW(fn,double,mt**nk);
1207  NNEW(stx,double,6*mi[0]**ne);
1208 
1209  /* determining the internal forces at the start of the increment
1210 
1211  for a static calculation with increased forced displacements
1212  the linear strains are calculated corresponding to
1213 
1214  the displacements at the end of the previous increment, extrapolated
1215  if appropriate (for nondispersive media) +
1216  the forced displacements at the end of the present increment +
1217  the temperatures at the end of the present increment (this sum is
1218  v) -
1219  the displacements at the end of the previous increment (this is vold)
1220 
1221  these linear strains are converted in stresses by multiplication
1222  with the tangent element stiffness matrix and converted into nodal
1223  forces.
1224 
1225  this boils down to the fact that the effect of forced displacements
1226  should be handled in a purely linear way at the
1227  start of a new increment, in order to speed up the convergence and
1228  (for dissipative media) guarantee smooth loading within the increment.
1229 
1230  for all other cases the nodal force calculation is based on
1231  the true stresses derived from the appropriate strain tensor taking
1232  into account the extrapolated displacements at the end of the
1233  previous increment + the forced displacements and the temperatures
1234  at the end of the present increment */
1235 
1236  iout=-1;
1237  iperturb_sav[0]=iperturb[0];
1238  iperturb_sav[1]=iperturb[1];
1239 
1240  /* first iteration in first increment: elastic tangent */
1241 
1242  if((*nmethod!=4)&&(ilin==1)){
1243 
1244  ielas=1;
1245 
1246  iperturb[0]=-1;
1247  iperturb[1]=0;
1248 
1249  for(k=0;k<neq[1];++k){b[k]=f[k];}
1250  NNEW(inum,ITG,*nk);
1251  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1252  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1253  ielorien,norien,orab,ntmat_,t1ini,t1act,ithermal,
1254  prestr,iprestr,filab,eme,emn,een,iperturb,
1255  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1256  ndirboun,xbounact,nboun,ipompc,
1257  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1258  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1259  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1260  &icmd, ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1261  emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1262  iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1263  fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1264  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1265  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1266  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1267  islavsurf,ielprop,prop);
1268  iperturb[0]=0;SFREE(inum);
1269 
1270  /* check whether any displacements or temperatures are changed
1271  in the new increment */
1272 
1273  for(k=0;k<neq[1];++k){f[k]=f[k]+b[k];}
1274 
1275  }
1276  else{
1277 
1278  NNEW(inum,ITG,*nk);
1279  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1280  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1281  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1282  prestr,iprestr,filab,eme,emn,een,iperturb,
1283  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1284  ndirboun,xbounact,nboun,ipompc,
1285  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1286  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1287  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1288  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1289  emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1290  iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1291  fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1292  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1293  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1294  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1295  islavsurf,ielprop,prop);
1296  SFREE(inum);
1297 
1298  memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1299 
1300  if(*ithermal!=2){
1301  for(k=0;k<6*mi[0]*ne0;++k){
1302  sti[k]=stx[k];
1303  }
1304  }
1305 
1306  }
1307 
1308  ielas=0;
1309  iout=0;
1310 
1311  SFREE(fn);SFREE(stx);SFREE(v);
1312 
1313  /***************************************************************/
1314  /* iteration counter and start of the loop over the iterations */
1315  /***************************************************************/
1316 
1317  iit=1;
1318  icntrl=0;
1319  ctrl[0]=i0ref;ctrl[1]=irref;ctrl[3]=icref;
1320  if(uncoupled){
1321  *ithermal=2;
1322  NNEW(iruc,ITG,nzs[1]-nzs[0]);
1323  for(k=0;k<nzs[1]-nzs[0];k++) iruc[k]=irow[k+nzs[0]]-neq[0];
1324  }
1325 
1326  while(icntrl==0){
1327 
1328  /* updating the nonlinear mpc's (also affects the boundary
1329  conditions through the nonhomogeneous part of the mpc's) */
1330 
1331  if((iit!=1)||((uncoupled)&&(*ithermal==1))){
1332 
1333  printf(" iteration %" ITGFORMAT "\n\n",iit);
1334 
1335  FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
1336  xloadold,xload,
1337  xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
1338  t1old,t1,t1act,iamt1,nk,amta,
1339  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
1340  xbounold,xboun,xbounact,iamboun,nboun,
1341  nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
1342  co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
1343  ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
1344  iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
1345 
1346  for(i=0;i<3;i++){cam[i]=0.;}for(i=3;i<5;i++){cam[i]=0.5;}
1347  if(*ithermal>1){radflowload(itg,ieg,&ntg,&ntr,adrad,aurad,bcr,ipivr,
1348  ac,bc,nload,sideload,nelemload,xloadact,lakon,ipiv,
1349  ntmat_,vold,shcon,nshcon,ipkon,kon,co,
1350  kontri,&ntri,nloadtr,tarea,tenv,physcon,erad,&adview,&auview,
1351  nflow,ikboun,xbounact,nboun,ithermal,&iinc,&iit,
1352  cs,mcs,inocs,&ntrit,nk,fenv,istep,&dtime,ttime,&time,ilboun,
1353  ikforc,ilforc,xforcact,nforc,cam,ielmat,&nteq,prop,ielprop,
1354  nactdog,nacteq,nodeboun,ndirboun,&network,
1355  rhcon,nrhcon,ipobody,ibody,xbodyact,nbody,iviewfile,jobnamef,
1356  ctrl,xloadold,&reltime,nmethod,set,mi,istartset,iendset,ialset,
1357  nset,ineighe,nmpc,nodempc,ipompc,coefmpc,labmpc,&iemchange,nam,
1358  iamload,jqrad,irowrad,&nzsrad,icolrad,ne,iaxial);
1359  }
1360 
1361  if((icascade==2)||
1362  ((ncont!=0)&&(ismallsliding==0))){
1363  memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
1364  RENEW(nodempc,ITG,3*memmpcref_);
1365  for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
1366  RENEW(coefmpc,double,memmpcref_);
1367  for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
1368  }
1369 
1370  if((ncont!=0)&&(*mortar<=1)&&(ismallsliding==0)&&((iit<=8)||(*mortar==1))){
1371  neold=*ne;
1372  *ne=ne0;*nkon=nkon0;
1373  contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1374  ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,
1375  straight,nkon,co,vold,ielmat,cs,elcon,istep,
1376  &iinc,&iit,ncmat_,ntmat_,&ne0,
1377  vini,nmethod,nmpc,&mpcfree,&memmpc_,&ipompc,&labmpc,
1378  &ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,iperturb,
1379  ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
1380  itiefac,areaslav,iponoels,inoels,springarea,tietol,
1381  &reltime,imastnode,nmastnode,xmastnor,
1382  filab,mcs,ics,&nasym,xnoels,mortar,pslavsurf,pmastsurf,
1383  clearini,&theta);
1384 
1385  /* check whether, for a dynamic calculation, damping is involved */
1386 
1387  if(*nmethod==4){
1388  if(*iexpl<=1){
1389  if(idampingwithoutcontact==0){
1390  for(i=0;i<*ne;i++){
1391  if(ipkon[i]<0) continue;
1392  if(*ncmat_>=5){
1393  if(strcmp1(&lakon[i*8],"ES")==0){
1394  if(strcmp1(&lakon[i*8+6],"C")==0){
1395  imat=ielmat[i*mi[2]];
1396  if(elcon[(*ncmat_+1)**ntmat_*(imat-1)+4]>0.){
1397  idamping=1;break;
1398  }
1399  }
1400  }
1401  }
1402  }
1403  }
1404  }
1405  }
1406 
1407  if(*mortar==0){
1408  if(*ne!=neold){iflagact=1;}
1409  }else if(*mortar==1){
1410  if(((*ne-ne0)<(neold-ne0)*0.999)||
1411  ((*ne-ne0)>(neold-ne0)*1.001)){iflagact=1;}
1412  }
1413 
1414  printf(" Number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
1415 
1416  if(icascade<1)icascade=1;
1417  }
1418 
1419  if(*ithermal==3){
1420  for(k=0;k<*nk;++k){t1act[k]=vold[mt*k];}
1421  }
1422 
1423  FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1424  nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1425  &maxlenmpc,ikmpc,ilmpc,&icascade,
1426  kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,&iit,
1427  &idiscon,&ncont,trab,ntrans,ithermal,mi));
1428 
1429  if((icascade==2)||
1430  ((ncont!=0)&&(ismallsliding==0))){
1431 // memmpcref_=memmpc_;mpcfreeref=mpcfree;
1432 // RENEW(nodempcref,ITG,3*memmpc_);
1433  for(k=0;k<3*memmpc_;k++){nodempcref[k]=nodempc[k];}
1434 // RENEW(coefmpcref,double,memmpc_);
1435  for(k=0;k<memmpc_;k++){coefmpcref[k]=coefmpc[k];}
1436  }
1437 
1438  if(icascade>0){
1439  remastruct(ipompc,&coefmpc,&nodempc,nmpc,
1440  &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
1441  labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
1442  kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
1443  neq,nzs,nmethod,&f,&fext,&b,&aux2,&fini,&fextini,
1444  &adb,&aub,ithermal,iperturb,mass,mi,iexpl,mortar,
1445  typeboun,&cv,&cvini,&iit);
1446 
1447  /* invert nactdof */
1448 
1449  SFREE(nactdofinv);
1450  NNEW(nactdofinv,ITG,mt**nk);
1451  NNEW(nodorig,ITG,*nk);
1452  FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
1453  ipkon,lakon,kon,ne));
1454  SFREE(nodorig);
1455 
1456  NNEW(v,double,mt**nk);
1457  NNEW(stx,double,6*mi[0]**ne);
1458  NNEW(fn,double,mt**nk);
1459 
1460  memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1461  iout=-1;
1462 
1463  NNEW(inum,ITG,*nk);
1464  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1465  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1466  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1467  prestr,iprestr,filab,eme,emn,een,iperturb,
1468  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1469  ndirboun,xbounact,nboun,ipompc,
1470  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1471  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1472  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1473  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1474  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1475  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1476  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1477  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1478  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1479  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1480  islavsurf,ielprop,prop);
1481 
1482  /*for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}*/
1483 
1484  SFREE(v);SFREE(stx);SFREE(fn);SFREE(inum);
1485  iout=0;
1486 
1487  }else{
1488 
1489  /*for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}*/
1490  }
1491  }
1492 
1493  if(*iexpl<=1){
1494 
1495  /* calculating the local stiffness matrix and external loading */
1496 
1497  NNEW(ad,double,neq[1]);
1498  NNEW(au,double,nzs[1]);
1499 
1500  FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounact,nboun,
1501  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1502  nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1503  nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
1504  nmethod,ikmpc,ilmpc,ikboun,ilboun,
1505  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1506  ielmat,ielorien,norien,orab,ntmat_,
1507  t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1508  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1509  xstiff,npmat_,&dtime,matname,mi,
1510  ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
1511  physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
1512  &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
1513  xstateini,xstate,thicke,integerglob,doubleglob,
1514  tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
1515  pmastsurf,mortar,clearini,ielprop,prop));
1516 
1517  if(nasym==1){
1518  RENEW(au,double,2*nzs[1]);
1519  if(*nmethod==4) RENEW(aub,double,2*nzs[1]);
1520  symmetryflag=2;
1521  inputformat=1;
1522 
1523  FORTRAN(mafillsmas,(co,nk,kon,ipkon,lakon,ne,nodeboun,
1524  ndirboun,xbounact,nboun,
1525  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1526  nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1527  nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
1528  nmethod,ikmpc,ilmpc,ikboun,ilboun,
1529  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1530  ielmat,ielorien,norien,orab,ntmat_,
1531  t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1532  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1533  xstiff,npmat_,&dtime,matname,mi,
1534  ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
1535  physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
1536  &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
1537  xstateini,xstate,thicke,
1538  integerglob,doubleglob,tieset,istartset,iendset,
1539  ialset,ntie,&nasym,pslavsurf,pmastsurf,mortar,clearini,
1540  ielprop,prop));
1541  }
1542 
1543 
1544  iperturb[0]=iperturb_sav[0];
1545  iperturb[1]=iperturb_sav[1];
1546 
1547  }else{
1548 
1549  /* calculating the external loading
1550 
1551  This is only done once per increment. In reality, the
1552  external loading is a function of vold (specifically,
1553  the body forces and surface loading). This effect is
1554  neglected, since the increment size in dynamic explicit
1555  calculations is usually small */
1556 
1557  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1558  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1559  nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,
1560  nbody,cgr,fext,nactdof,&neq[1],
1561  nmethod,ikmpc,ilmpc,
1562  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1563  ielmat,ielorien,norien,orab,ntmat_,
1564  t0,t1act,ithermal,iprestr,vold,iperturb,
1565  iexpl,plicon,nplicon,plkcon,nplkcon,
1566  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1567  xbodyold,&reltime,veold,matname,mi,ikactmech,
1568  &nactmech,ielprop,prop));
1569  }
1570 
1571 /* for(k=0;k<neq[1];++k){printf("f=%" ITGFORMAT ",%f\n",k,f[k]);}
1572  for(k=0;k<neq[1];++k){printf("fext=%" ITGFORMAT ",%f\n",k,fext[k]);}
1573  for(k=0;k<neq[1];++k){printf("ad=%" ITGFORMAT ",%f\n",k,ad[k]);}
1574  for(k=0;k<nzs[1];++k){printf("au=%" ITGFORMAT ",%f\n",k,au[k]);}*/
1575 
1576  /* calculating the damping matrix for implicit dynamic
1577  calculations */
1578 
1579  if(idamping==1){
1580 
1581  /* Rayleigh damping */
1582 
1583  NNEW(adc,double,neq[1]);
1584  for(k=0;k<neq[0];k++){adc[k]=alpham*adb[k]+betam*ad[k];}
1585  if(nasym==0){
1586  NNEW(auc,double,nzs[1]);
1587  for(k=0;k<nzs[0];k++){auc[k]=alpham*aub[k]+betam*au[k];}
1588  }else{
1589  NNEW(auc,double,2*nzs[1]);
1590  for(k=0;k<2*nzs[0];k++){auc[k]=alpham*aub[k]+betam*au[k];}
1591  }
1592 
1593  /* dashpots and contact damping */
1594 
1595  FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
1596  ndirboun,xbounact,nboun,
1597  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1598  nforc,nelemload,sideload,xloadact,nload,xbodyact,
1599  ipobody,nbody,cgr,
1600  adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
1601  ikmpc,ilmpc,ikboun,ilboun,
1602  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1603  ielorien,norien,orab,ntmat_,
1604  t0,t1act,ithermal,prestr,iprestr,vold,iperturb,sti,
1605  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
1606  xstiff,npmat_,&dtime,matname,mi,ncmat_,
1607  ttime,&time,istep,&iinc,ibody,clearini,mortar,springarea,
1608  pslavsurf,pmastsurf,&reltime,&nasym));
1609  }
1610 
1611  /* calculating the residual */
1612 
1613  calcresidual(nmethod,neq,b,fext,f,iexpl,nactdof,aux2,vold,
1614  vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1615  islavnode,nslavnode,mortar,ntie,f_cm,f_cs,mi,
1616  nzs,&nasym,&idamping,veold,adc,auc,cvini,cv);
1617 
1618  newstep=0;
1619 
1620  if(*nmethod==0){
1621 
1622  /* error occurred in mafill: storing the geometry in frd format */
1623 
1624  *nmethod=0;
1625  ++*kode;
1626  NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
1627  if(strcmp1(&filab[1044],"ZZS")==0){
1628  NNEW(neigh,ITG,40**ne);
1629  NNEW(ipneigh,ITG,*nk);
1630  }
1631 
1632  ptime=*ttime+time;
1633  frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
1634  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1635  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1636  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1637  mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1638  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1639  thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
1640 
1641  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1642  #ifdef COMPANY
1643  FORTRAN(uout,(v,mi,ithermal));
1644 #endif
1645  SFREE(inum);FORTRAN(stop,());
1646 
1647  }
1648 
1649  /* implicit step (static or dynamic */
1650 
1651  if(*iexpl<=1){
1652  if((*nmethod==4)&&(*mortar<2)){
1653 
1654  /* mechanical part */
1655 
1656  if(*ithermal!=2){
1657  scal1=bet*dtime*dtime*(1.+*alpha);
1658  for(k=0;k<neq[0];++k){
1659  ad[k]=adb[k]+scal1*ad[k];
1660  }
1661  for(k=0;k<nzs[0];++k){
1662  au[k]=aub[k]+scal1*au[k];
1663  }
1664 
1665  /* upper triangle of asymmetric matrix */
1666 
1667  if(nasym>0){
1668  for(k=nzs[2];k<nzs[2]+nzs[0];++k){
1669  au[k]=aub[k]+scal1*au[k];
1670  }
1671  }
1672 
1673  /* damping */
1674 
1675  if(idamping==1){
1676  scal1=gam*dtime*(1.+*alpha);
1677  for(k=0;k<neq[0];++k){
1678  ad[k]+=scal1*adc[k];
1679  }
1680  for(k=0;k<nzs[0];++k){
1681  au[k]+=scal1*auc[k];
1682  }
1683 
1684  /* upper triangle of asymmetric matrix */
1685 
1686  if(nasym>0){
1687  for(k=nzs[2];k<nzs[2]+nzs[0];++k){
1688  au[k]+=scal1*auc[k];
1689  }
1690  }
1691  }
1692 
1693  }
1694 
1695  /* thermal part */
1696 
1697  if(*ithermal>1){
1698  for(k=neq[0];k<neq[1];++k){
1699  ad[k]=adb[k]/dtime+ad[k];
1700  }
1701  for(k=nzs[0];k<nzs[1];++k){
1702  au[k]=aub[k]/dtime+au[k];
1703  }
1704 
1705  /* upper triangle of asymmetric matrix */
1706 
1707  if(nasym>0){
1708  for(k=nzs[2]+nzs[0];k<nzs[2]+nzs[1];++k){
1709  au[k]=aub[k]/dtime+au[k];
1710  }
1711  }
1712  }
1713  }
1714 
1715  if(*isolver==0){
1716 #ifdef SPOOLES
1717  if(*ithermal<2){
1718  spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1719  &symmetryflag,&inputformat,&nzs[2]);
1720  }else if((*ithermal==2)&&(uncoupled)){
1721  n1=neq[1]-neq[0];
1722  n2=nzs[1]-nzs[0];
1723  spooles(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1724  &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1725  &n1,&n2,&symmetryflag,&inputformat,&nzs[2]);
1726  }else{
1727  spooles(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1728  &symmetryflag,&inputformat,&nzs[2]);
1729  }
1730 #else
1731  printf(" *ERROR in nonlingeo: the SPOOLES library is not linked\n\n");
1732  FORTRAN(stop,());
1733 #endif
1734  }
1735  else if((*isolver==2)||(*isolver==3)){
1736  if(nasym>0){
1737  printf(" *ERROR in nonlingeo: the iterative solver cannot be used for asymmetric matrices\n\n");
1738  FORTRAN(stop,());
1739  }
1740  preiter(ad,&au,b,&icol,&irow,&neq[1],&nzs[1],isolver,iperturb);
1741  }
1742  else if(*isolver==4){
1743 #ifdef SGI
1744  if(nasym>0){
1745  printf(" *ERROR in nonlingeo: the SGI solver cannot be used for asymmetric matrices\n\n");
1746  FORTRAN(stop,());
1747  }
1748  token=1;
1749  if(*ithermal<2){
1750  sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],token);
1751  }else if((*ithermal==2)&&(uncoupled)){
1752  n1=neq[1]-neq[0];
1753  n2=nzs[1]-nzs[0];
1754  sgi_main(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1755  &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1756  &n1,&n2,token);
1757  }else{
1758  sgi_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],token);
1759  }
1760 #else
1761  printf(" *ERROR in nonlingeo: the SGI library is not linked\n\n");
1762  FORTRAN(stop,());
1763 #endif
1764  }
1765  else if(*isolver==5){
1766  if(nasym>0){
1767  printf(" *ERROR in nonlingeo: the TAUCS solver cannot be used for asymmetric matrices\n\n");
1768  FORTRAN(stop,());
1769  }
1770 #ifdef TAUCS
1771  tau(ad,&au,adb,aub,&sigma,b,icol,&irow,&neq[1],&nzs[1]);
1772 #else
1773  printf(" *ERROR in nonlingeo: the TAUCS library is not linked\n\n");
1774  FORTRAN(stop,());
1775 #endif
1776  }
1777  else if(*isolver==7){
1778 #ifdef PARDISO
1779  if(*ithermal<2){
1780  pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[0],&nzs[0],
1781  &symmetryflag,&inputformat,jq,&nzs[2]);
1782  }else if((*ithermal==2)&&(uncoupled)){
1783  n1=neq[1]-neq[0];
1784  n2=nzs[1]-nzs[0];
1785  pardiso_main(&ad[neq[0]],&au[nzs[0]],&adb[neq[0]],&aub[nzs[0]],
1786  &sigma,&b[neq[0]],&icol[neq[0]],iruc,
1787  &n1,&n2,&symmetryflag,&inputformat,jq,&nzs[2]);
1788  }else{
1789  pardiso_main(ad,au,adb,aub,&sigma,b,icol,irow,&neq[1],&nzs[1],
1790  &symmetryflag,&inputformat,jq,&nzs[2]);
1791  }
1792 #else
1793  printf(" *ERROR in nonlingeo: the PARDISO library is not linked\n\n");
1794  FORTRAN(stop,());
1795 #endif
1796  }
1797 
1798  if(*mortar<=1){SFREE(ad);SFREE(au);}
1799  }
1800 
1801  /* explicit dynamic step */
1802 
1803  else{
1804  if(*ithermal!=2){
1805  for(k=0;k<neq[0];++k){
1806  b[k]=b[k]/adb[k];
1807  }
1808  }
1809  if(*ithermal>1){
1810  for(k=neq[0];k<neq[1];++k){
1811  b[k]=b[k]*dtime/adb[k];
1812  }
1813  }
1814  }
1815 
1816  /* calculating the displacements, stresses and forces */
1817 
1818  NNEW(v,double,mt**nk);
1819  memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1820 
1821  NNEW(stx,double,6*mi[0]**ne);
1822  NNEW(fn,double,mt**nk);
1823 
1824  NNEW(inum,ITG,*nk);
1825  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1826  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1827  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
1828  prestr,iprestr,filab,eme,emn,een,iperturb,
1829  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
1830  ndirboun,xbounact,nboun,ipompc,
1831  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1832  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1833  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1834  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1835  emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1836  iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
1837  fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1838  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
1839  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1840  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1841  islavsurf,ielprop,prop);
1842  SFREE(inum);
1843 
1844  if(*ithermal!=2){
1845  if(cam[0]>uam[0]){uam[0]=cam[0];}
1846  if(qau<1.e-10){
1847  if(qa[0]>ea*qam[0]){qam[0]=(qamold[0]*jnz+qa[0])/(jnz+1);}
1848  else {qam[0]=qamold[0];}
1849  }
1850  }
1851  if(*ithermal>1){
1852  if(cam[1]>uam[1]){uam[1]=cam[1];}
1853  if(qau<1.e-10){
1854  if(qa[1]>ea*qam[1]){qam[1]=(qamold[1]*jnz+qa[1])/(jnz+1);}
1855  else {qam[1]=qamold[1];}
1856  }
1857  }
1858 
1859  memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1860  if(*ithermal!=2){
1861  for(k=0;k<6*mi[0]*ne0;++k){
1862  sti[k]=stx[k];
1863  }
1864  }
1865 
1866  /* calculating the ratio of the smallest to largest pressure
1867  for face-to-face contact
1868  only done at the end of a step */
1869 
1870  if((*mortar==1)&&(1.-theta-dtheta<=1.e-6)){
1871  FORTRAN(negativepressure,(&ne0,ne,mi,stx,&pressureratio));
1872  }else{pressureratio=0.;}
1873 
1874  SFREE(v);SFREE(stx);SFREE(fn);
1875 
1876  /* calculating the residual */
1877 
1878  calcresidual(nmethod,neq,b,fext,f,iexpl,nactdof,aux2,vold,
1879  vini,&dtime,accold,nk,adb,aub,jq,irow,nzl,alpha,fextini,fini,
1880  islavnode,nslavnode,mortar,ntie,f_cm,f_cs,mi,
1881  nzs,&nasym,&idamping,veold,adc,auc,cvini,cv);
1882 
1883  if(idamping==1){SFREE(adc);SFREE(auc);}
1884 
1885  /* calculating the maximum residual */
1886 
1887  for(k=0;k<2;++k){
1888  ram2[k]=ram1[k];
1889  ram1[k]=ram[k];
1890  ram[k]=0.;
1891  }
1892  if(*ithermal!=2){
1893  for(k=0;k<neq[0];++k){
1894  err=fabs(b[k]);
1895  if(err>ram[0]){ram[0]=err;ram[2]=k+0.5;}
1896  }
1897  }
1898  if(*ithermal>1){
1899  for(k=neq[0];k<neq[1];++k){
1900  err=fabs(b[k]);
1901  if(err>ram[1]){ram[1]=err;ram[3]=k+0.5;}
1902  }
1903  }
1904 
1905  /* Divergence criteria for face-to-face penalty is different */
1906 
1907  if(*mortar==1){
1908  for(k=4;k<8;++k){
1909  ram2[k]=ram1[k];
1910  ram1[k]=ram[k];
1911  }
1912  ram[4]=ram[0]+ram1[0];
1913  if((iflagact==0)&&(iit>1)){
1914  ram[5]=1.5;
1915  }else{ram[5]=0.5;}
1916  ram[6]=abs((neold-ne0)-(*ne-ne0))+0.5;
1917  if(iit>3){
1918  if((ram[6]>=ram1[6])&&(ram[6]>=ram2[6])){
1919  ram[7]=1.5;
1920  }else{ram[7]=0.5;}
1921  }
1922 
1923 /* for(k=2;k<4;++k){
1924  ram2[k]=ram1[k];
1925  ram1[k]=ram[k+2];
1926  }
1927  ram[4]=ram[0]+ram1[0];
1928 
1929  ram[5]=abs((neold-ne0)-(*ne-ne0));
1930 
1931  if(iit==2){ram[6]=ram[4];}
1932  if((iit>=3)&&(ram[6]>ram[4])){ram[6]=ram[4];}*/
1933  }
1934 
1935  /* next line is inserted to cope with stress-less
1936  temperature calculations */
1937 
1938  if(*ithermal!=2){
1939  if(ram[0]<1.e-6) ram[0]=0.;
1940  printf(" average force= %f\n",qa[0]);
1941  printf(" time avg. forc= %f\n",qam[0]);
1942  if((ITG)((double)nactdofinv[(ITG)ram[2]]/mt)+1==0){
1943  printf(" largest residual force= %f\n",
1944  ram[0]);
1945  }else{
1946  inode=(ITG)((double)nactdofinv[(ITG)ram[2]]/mt)+1;
1947  idir=nactdofinv[(ITG)ram[2]]-mt*(inode-1);
1948  printf(" largest residual force= %f in node %" ITGFORMAT " and dof %" ITGFORMAT "\n",
1949  ram[0],inode,idir);
1950  }
1951  printf(" largest increment of disp= %e\n",uam[0]);
1952  if((ITG)cam[3]==0){
1953  printf(" largest correction to disp= %e\n\n",
1954  cam[0]);
1955  }else{
1956  inode=(ITG)((double)nactdofinv[(ITG)cam[3]]/mt)+1;
1957  idir=nactdofinv[(ITG)cam[3]]-mt*(inode-1);
1958  printf(" largest correction to disp= %e in node %" ITGFORMAT " and dof %" ITGFORMAT "\n\n",cam[0],inode,idir);
1959  }
1960  }
1961  if(*ithermal>1){
1962  if(ram[1]<1.e-6) ram[1]=0.;
1963  printf(" average flux= %f\n",qa[1]);
1964  printf(" time avg. flux= %f\n",qam[1]);
1965  if((ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1==0){
1966  printf(" largest residual flux= %f\n",
1967  ram[1]);
1968  }else{
1969  inode=(ITG)((double)nactdofinv[(ITG)ram[3]]/mt)+1;
1970  idir=nactdofinv[(ITG)ram[3]]-mt*(inode-1);
1971  printf(" largest residual flux= %f in node %" ITGFORMAT " and dof %" ITGFORMAT "\n",ram[1],inode,idir);
1972  }
1973  printf(" largest increment of temp= %e\n",uam[1]);
1974  if((ITG)cam[4]==0){
1975  printf(" largest correction to temp= %e\n\n",
1976  cam[1]);
1977  }else{
1978  inode=(ITG)((double)nactdofinv[(ITG)cam[4]]/mt)+1;
1979  idir=nactdofinv[(ITG)cam[4]]-mt*(inode-1);
1980  printf(" largest correction to temp= %e in node %" ITGFORMAT " and dof %" ITGFORMAT "\n\n",cam[1],inode,idir);
1981  }
1982  }
1983  fflush(stdout);
1984 
1985  FORTRAN(writecvg,(istep,&iinc,&iit,ne,&ne0,ram,qam,cam,uam,
1986  ithermal));
1987 
1988  checkconvergence(co,nk,kon,ipkon,lakon,ne,stn,nmethod,
1989  kode,filab,een,t1act,&time,epn,ielmat,matname,enern,
1990  xstaten,nstate_,istep,&iinc,iperturb,ener,mi,output,
1991  ithermal,qfn,&mode,&noddiam,trab,inotr,ntrans,orab,
1992  ielorien,norien,description,sti,&icutb,&iit,&dtime,qa,
1993  vold,qam,ram1,ram2,ram,cam,uam,&ntg,ttime,&icntrl,
1994  &theta,&dtheta,veold,vini,idrct,tper,&istab,tmax,
1995  nactdof,b,tmin,ctrl,amta,namta,itpamp,inext,&dthetaref,
1996  &itp,&jprint,jout,&uncoupled,t1,&iitterm,nelemload,
1997  nload,nodeboun,nboun,itg,ndirboun,&deltmx,&iflagact,
1998  set,nset,istartset,iendset,ialset,emn,thicke,jobnamec,
1999  mortar,nmat,ielprop,prop);
2000 
2001  }
2002 
2003  /*********************************************************/
2004  /* end of the iteration loop */
2005  /*********************************************************/
2006 
2007  /* icutb=0 means that the iterations in the increment converged,
2008  icutb!=0 indicates that the increment has to be reiterated with
2009  another increment size (dtheta) */
2010 
2011  if(uncoupled){
2012  SFREE(iruc);
2013  }
2014 
2015  if(((qa[0]>ea*qam[0])||(qa[1]>ea*qam[1]))&&(icutb==0)){jnz++;}
2016  iit=0;
2017 
2018  if(icutb!=0){
2019  memcpy(&vold[0],&vini[0],sizeof(double)*mt**nk);
2020 
2021  for(k=0;k<*nboun;++k){xbounact[k]=xbounini[k];}
2022  if((*ithermal==1)||(*ithermal>=3)){
2023  for(k=0;k<*nk;++k){t1act[k]=t1ini[k];}
2024  }
2025  for(k=0;k<neq[1];++k){
2026  f[k]=fini[k];
2027  }
2028  if(*nmethod==4){
2029  for(k=0;k<mt**nk;++k){
2030  veold[k]=veini[k];
2031  accold[k]=accini[k];
2032  }
2033  for(k=0;k<neq[1];++k){
2034 // f[k]=fini[k];
2035  fext[k]=fextini[k];
2036  cv[k]=cvini[k];
2037  }
2038  }
2039  if(*ithermal!=2){
2040  for(k=0;k<6*mi[0]*ne0;++k){
2041  sti[k]=stiini[k];
2042  eme[k]=emeini[k];
2043  }
2044  }
2045  if(*nener==1)
2046  for(k=0;k<mi[0]*ne0;++k){ener[k]=enerini[k];}
2047 
2048  for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
2049  xstate[k]=xstateini[k];
2050  }
2051 
2052  qam[0]=qamold[0];
2053  qam[1]=qamold[1];
2054  }
2055 
2056  /* face-to-face penalty */
2057 
2058  if((*mortar==1)&&(icutb==0)){
2059 
2060  ntrimax=0;
2061  for(i=0;i<*ntie;i++){
2062  if(itietri[2*i+1]-itietri[2*i]+1>ntrimax)
2063  ntrimax=itietri[2*i+1]-itietri[2*i]+1;
2064  }
2065  NNEW(xo,double,ntrimax);
2066  NNEW(yo,double,ntrimax);
2067  NNEW(zo,double,ntrimax);
2068  NNEW(x,double,ntrimax);
2069  NNEW(y,double,ntrimax);
2070  NNEW(z,double,ntrimax);
2071  NNEW(nx,ITG,ntrimax);
2072  NNEW(ny,ITG,ntrimax);
2073  NNEW(nz,ITG,ntrimax);
2074 
2075  /* Determination of active nodes (islavact) */
2076 
2077  FORTRAN(islavactive,(tieset,ntie,itietri,cg,straight,
2078  co,vold,xo,yo,zo,x,y,z,nx,ny,nz,mi,
2079  imastop,nslavnode,islavnode,islavact));
2080 
2081  SFREE(xo);SFREE(yo);SFREE(zo);SFREE(x);SFREE(y);SFREE(z);SFREE(nx);
2082  SFREE(ny);SFREE(nz);
2083 
2084  if(negpres==0){
2085  if((*mortar==1)&&(1.-theta-dtheta<=1.e-6)&&(itruecontact==1)){
2086  printf(" pressure ratio (smallest/largest pressure over all contact areas) =%e\n\n",pressureratio);
2087  if(pressureratio<-0.05){
2088  printf(" zero-size increment is appended\n\n");
2089  negpres=1;dtheta=0.;
2090  }
2091  }
2092  }else{negpres=0;}
2093 
2094  }
2095 
2096  /* output */
2097 
2098  if((jout[0]==jprint)&&(icutb==0)){
2099 
2100  jprint=0;
2101 
2102  /* calculating the displacements and the stresses and storing */
2103  /* the results in frd format */
2104 
2105  NNEW(v,double,mt**nk);
2106  NNEW(fn,double,mt**nk);
2107  NNEW(stn,double,6**nk);
2108  if(*ithermal>1) NNEW(qfn,double,3**nk);
2109  NNEW(inum,ITG,*nk);
2110  NNEW(stx,double,6*mi[0]**ne);
2111 
2112  if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
2113  if(strcmp1(&filab[435],"PEEQ")==0) NNEW(epn,double,*nk);
2114  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2115  if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
2116  if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
2117  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
2118 
2119  memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
2120 
2121  iout=2;
2122  icmd=3;
2123 
2124 #ifdef COMPANY
2125  FORTRAN(uinit,());
2126 #endif
2127  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
2128  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2129  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
2130  prestr,iprestr,filab,eme,emn,een,iperturb,
2131  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
2132  ndirboun,xbounact,nboun,ipompc,
2133  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
2134  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
2135  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
2136  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
2137  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
2138  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
2139  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
2140  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
2141  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2142  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
2143  islavsurf,ielprop,prop);
2144 
2145  memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
2146 
2147  iout=0;
2148  if(*iexpl<=1) icmd=0;
2149 // FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
2150 // for(k=0;k<ntg;k++)if(inum[itg[k]-1]>0){inum[itg[k]-1]*=-1;}
2151 
2152  ++*kode;
2153  if(*mcs!=0){
2154  ptime=*ttime+time;
2155  frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
2156  t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
2157  nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
2158  ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
2159  norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,&ne0,
2160  cdn,mortar,nmat);
2161 #ifdef COMPANY
2162  FORTRAN(uout,(v,mi,ithermal));
2163 #endif
2164  }
2165  else{
2166  if(strcmp1(&filab[1044],"ZZS")==0){
2167  NNEW(neigh,ITG,40**ne);
2168  NNEW(ipneigh,ITG,*nk);
2169  }
2170 
2171  ptime=*ttime+time;
2172  frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
2173  kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2174  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2175  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2176  mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
2177  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2178  thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
2179 
2180  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2181 #ifdef COMPANY
2182  FORTRAN(uout,(v,mi,ithermal));
2183 #endif
2184  }
2185 
2186  SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
2187  if(*ithermal>1){SFREE(qfn);}
2188 
2189  if(strcmp1(&filab[261],"E ")==0) SFREE(een);
2190  if(strcmp1(&filab[435],"PEEQ")==0) SFREE(epn);
2191  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2192  if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
2193  if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
2194  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
2195  }
2196 
2197  }
2198 
2199  /*********************************************************/
2200  /* end of the increment loop */
2201  /*********************************************************/
2202 
2203  if(jprint!=0){
2204 
2205  /* calculating the displacements and the stresses and storing
2206  the results in frd format */
2207 
2208  NNEW(v,double,mt**nk);
2209  NNEW(fn,double,mt**nk);
2210  NNEW(stn,double,6**nk);
2211  if(*ithermal>1) NNEW(qfn,double,3**nk);
2212  NNEW(inum,ITG,*nk);
2213  NNEW(stx,double,6*mi[0]**ne);
2214 
2215  if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
2216  if(strcmp1(&filab[435],"PEEQ")==0) NNEW(epn,double,*nk);
2217  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2218  if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
2219  if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
2220  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
2221 
2222  memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
2223  iout=2;
2224  icmd=3;
2225 
2226 #ifdef COMPANY
2227  FORTRAN(uinit,());
2228 #endif
2229  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
2230  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2231  ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
2232  prestr,iprestr,filab,eme,emn,een,iperturb,
2233  f,fn,nactdof,&iout,qa,vold,b,nodeboun,
2234  ndirboun,xbounact,nboun,ipompc,
2235  nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
2236  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
2237  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
2238  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
2239  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
2240  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
2241  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
2242  &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
2243  sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2244  mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
2245  islavsurf,ielprop,prop);
2246 
2247  memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
2248 
2249  iout=0;
2250  if(*iexpl<=1) icmd=0;
2251 // FORTRAN(networkinum,(ipkon,inum,kon,lakon,ne,itg,&ntg));
2252 // for(k=0;k<ntg;k++)if(inum[itg[k]-1]>0){inum[itg[k]-1]*=-1;}
2253 
2254  ++*kode;
2255  if(*mcs>0){
2256  ptime=*ttime+time;
2257  frdcyc(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,kode,filab,een,
2258  t1act,fn,&ptime,epn,ielmat,matname,cs,mcs,nkon,enern,xstaten,
2259  nstate_,istep,&iinc,iperturb,ener,mi,output,ithermal,qfn,
2260  ialset,istartset,iendset,trab,inotr,ntrans,orab,ielorien,
2261  norien,stx,veold,&noddiam,set,nset,emn,thicke,jobnamec,&ne0,
2262  cdn,mortar,nmat);
2263 #ifdef COMPANY
2264  FORTRAN(uout,(v,mi,ithermal));
2265 #endif
2266 
2267  }
2268  else{
2269  if(strcmp1(&filab[1044],"ZZS")==0){
2270  NNEW(neigh,ITG,40**ne);
2271  NNEW(ipneigh,ITG,*nk);
2272  }
2273 
2274  ptime=*ttime+time;
2275  frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,nmethod,
2276  kode,filab,een,t1act,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2277  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2278  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2279  mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
2280  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2281  thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
2282 
2283  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2284 #ifdef COMPANY
2285  FORTRAN(uout,(v,mi,ithermal));
2286 #endif
2287  }
2288 
2289  SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
2290  if(*ithermal>1){SFREE(qfn);}
2291 
2292  if(strcmp1(&filab[261],"E ")==0) SFREE(een);
2293  if(strcmp1(&filab[435],"PEEQ")==0) SFREE(epn);
2294  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2295  if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
2296  if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
2297  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
2298 
2299  }
2300 
2301  /* setting the velocity to zero at the end of a quasistatic or stationary
2302  step */
2303 
2304  if(abs(*nmethod)==1){
2305  for(k=0;k<mt**nk;++k){veold[k]=0.;}
2306  }
2307 
2308  /* updating the loading at the end of the step;
2309  important in case the amplitude at the end of the step
2310  is not equal to one */
2311 
2312  for(k=0;k<*nboun;++k){
2313 
2314  /* thermal boundary conditions are updated only if the
2315  step was thermal or thermomechanical */
2316 
2317  if(ndirboun[k]==0){
2318  if(*ithermal<2) continue;
2319 
2320  /* mechanical boundary conditions are updated only
2321  if the step was not thermal or the node is a
2322  network node */
2323 
2324  }else if((ndirboun[k]>0)&&(ndirboun[k]<4)){
2325  node=nodeboun[k];
2326  FORTRAN(nident,(itg,&node,&ntg,&id));
2327  networknode=0;
2328  if(id>0){
2329  if(itg[id-1]==node) networknode=1;
2330  }
2331  if((*ithermal==2)&&(networknode==0)) continue;
2332  }
2333  xbounold[k]=xbounact[k];
2334  }
2335  for(k=0;k<*nforc;++k){xforcold[k]=xforcact[k];}
2336  for(k=0;k<2**nload;++k){xloadold[k]=xloadact[k];}
2337  for(k=0;k<7**nbody;k=k+7){xbodyold[k]=xbodyact[k];}
2338  if(*ithermal==1){
2339  for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
2340  for(k=0;k<*nk;++k){vold[mt*k]=t1act[k];}
2341  }
2342  else if(*ithermal>1){
2343  for(k=0;k<*nk;++k){t1[k]=vold[mt*k];}
2344  if(*ithermal>=3){
2345  for(k=0;k<*nk;++k){t1old[k]=t1act[k];}
2346  }
2347  }
2348 
2349  qaold[0]=qa[0];
2350  qaold[1]=qa[1];
2351 
2352  SFREE(f);SFREE(b);
2353  SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(xbodyact);
2354  if(*nbody>0) SFREE(ipobody);if(inewton==1){SFREE(cgr);}
2355  SFREE(fext);SFREE(ampli);SFREE(xbounini);SFREE(xstiff);
2356  if((*ithermal==1)||(*ithermal>=3)){SFREE(t1act);SFREE(t1ini);}
2357 
2358  if(*ithermal>1){
2359  SFREE(itg);SFREE(ieg);SFREE(kontri);SFREE(nloadtr);
2360  SFREE(nactdog);SFREE(nacteq);SFREE(ineighe);
2361  SFREE(tarea);SFREE(tenv);SFREE(fenv);SFREE(qfx);
2362  SFREE(erad);SFREE(ac);SFREE(bc);SFREE(ipiv);
2363  SFREE(bcr);SFREE(ipivr);SFREE(adview);SFREE(auview);SFREE(adrad);
2364  SFREE(aurad);SFREE(irowrad);SFREE(jqrad);SFREE(icolrad);
2365  if((*mcs>0)&&(ntr>0)){SFREE(inocs);}
2366  }
2367 
2368  if(icfd==1){
2369  SFREE(neifa);SFREE(neiel);SFREE(neij);SFREE(ielfa);SFREE(ifaext);
2370  SFREE(vel);SFREE(vfa);SFREE(nactdoh);SFREE(nactdohinv);
2371  SFREE(ipkonf);SFREE(lakonf);SFREE(ielmatf);
2372  if(*norien>0) SFREE(ielorienf);
2373  }
2374 
2375  SFREE(fini);
2376  if(*nmethod==4){
2377  SFREE(aux2);SFREE(fextini);SFREE(veini);SFREE(accini);
2378  SFREE(adb);SFREE(aub);SFREE(cvini);SFREE(cv);
2379  }
2380  SFREE(eei);SFREE(stiini);SFREE(emeini);
2381  if(*nener==1)SFREE(enerini);
2382  if(*nstate_!=0){SFREE(xstateini);}
2383 
2384  SFREE(aux);SFREE(iaux);SFREE(vini);
2385 
2386  if((icascade==2)||(ncont!=0)){
2387  memmpc_=memmpcref_;mpcfree=mpcfreeref;maxlenmpc=maxlenmpcref;
2388  RENEW(nodempc,ITG,3*memmpcref_);
2389  for(k=0;k<3*memmpcref_;k++){nodempc[k]=nodempcref[k];}
2390  RENEW(coefmpc,double,memmpcref_);
2391  for(k=0;k<memmpcref_;k++){coefmpc[k]=coefmpcref[k];}
2392  SFREE(nodempcref);SFREE(coefmpcref);
2393  }
2394 
2395  if(ncont!=0){
2396  *ne=ne0;*nkon=nkon0;
2397  if(*nener==1){
2398  RENEW(ener,double,mi[0]**ne*2);
2399  }
2400  RENEW(ipkon,ITG,*ne);
2401  RENEW(lakon,char,8**ne);
2402  RENEW(kon,ITG,*nkon);
2403  if(*norien>0){
2404  RENEW(ielorien,ITG,mi[2]**ne);
2405  }
2406  RENEW(ielmat,ITG,mi[2]**ne);
2407 
2408  SFREE(cg);SFREE(straight);
2409  SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
2410  SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
2411  SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
2412  SFREE(springarea);SFREE(xmastnor);
2413 
2414  if(*mortar==0){
2415  SFREE(areaslav);
2416  }else if(*mortar==1){
2417  SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
2418  SFREE(islavact);
2419  }
2420  }
2421 
2422  /* reset icascade */
2423 
2424  if(icascade==1){icascade=0;}
2425 
2426  mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
2427  mpcinfo[3]=maxlenmpc;
2428 
2429  if(iglob==1){SFREE(integerglob);SFREE(doubleglob);}
2430 
2431  *icolp=icol;*irowp=irow;*cop=co;*voldp=vold;
2432 
2433  *ipompcp=ipompc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2434  *fmpcp=fmpc;*nodempcp=nodempc;*coefmpcp=coefmpc;
2435 
2436  *ipkonp=ipkon;*lakonp=lakon;*konp=kon;*ielorienp=ielorien;
2437  *ielmatp=ielmat;*enerp=ener;*xstatep=xstate;
2438 
2439  *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
2440 
2441  (*tmin)*=(*tper);
2442  (*tmax)*=(*tper);
2443 
2444  SFREE(nactdofinv);
2445 
2446  (*ttime)+=(*tper);
2447 
2448  return;
2449 }
subroutine negativepressure(ne0, ne, mi, stx, pressureratio)
Definition: negativepressure.f:19
#define ITGFORMAT
Definition: CalculiX.h:52
void pardiso_main(double *ad, double *au, double *adb, double *aub, double *sigma, double *b, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG *symmetryflag, ITG *inputformat, ITG *jq, ITG *nzs3)
void frd(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne0, double *v, double *stn, ITG *inum, ITG *nmethod, ITG *kode, char *filab, double *een, double *t1, double *fn, double *time, double *epn, ITG *ielmat, char *matname, double *enern, double *xstaten, ITG *nstate_, ITG *istep, ITG *iinc, ITG *ithermal, double *qfn, ITG *mode, ITG *noddiam, double *trab, ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien, char *description, ITG *ipneigh, ITG *neigh, ITG *mi, double *stx, double *vr, double *vi, double *stnr, double *stni, double *vmax, double *stnmax, ITG *ngraph, double *veold, double *ener, ITG *ne, double *cs, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, double *eenmax, double *fnr, double *fni, double *emn, double *thicke, char *jobnamec, char *output, double *qfx, double *cdn, ITG *mortar, double *cdnr, double *cdni, ITG *nmat)
Definition: frd.c:27
void remastruct(ITG *ipompc, double **coefmpcp, ITG **nodempcp, ITG *nmpc, ITG *mpcfree, ITG *nodeboun, ITG *ndirboun, ITG *nboun, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun, char *labmpc, ITG *nk, ITG *memmpc_, ITG *icascade, ITG *maxlenmpc, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, ITG *nactdof, ITG *icol, ITG *jq, ITG **irowp, ITG *isolver, ITG *neq, ITG *nzs, ITG *nmethod, double **fp, double **fextp, double **bp, double **aux2p, double **finip, double **fextinip, double **adbp, double **aubp, ITG *ithermal, ITG *iperturb, ITG *mass, ITG *mi, ITG *iexpl, ITG *mortar, char *typeboun, double **cvp, double **cvinip, ITG *iit)
Definition: remastruct.c:24
subroutine checktruecontact(ntie, tieset, tietol, elcon, itruecontact, ncmat_, ntmat_)
Definition: checktruecontact.f:19
subroutine mafillsmas(co, nk, kon, ipkon, lakon, ne, nodeboun, ndirboun, xboun, nboun, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, ad, au, bb, nactdof, icol, jq, irow, neq, nzl, nmethod, ikmpc, ilmpc, ikboun, ilboun, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, vold, iperturb, sti, nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, mass, stiffness, buckling, rhsi, intscheme, physcon, shcon, nshcon, cocon, ncocon, ttime, time, istep, iinc, coriolis, ibody, xloadold, reltime, veold, springarea, nstate_, xstateini, xstate, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, pslavsurf, pmastsurf, mortar, clearini, ielprop, prop)
Definition: mafillsmas.f:19
void preiter(double *ad, double **aup, double *b, ITG **icolp, ITG **irowp, ITG *neq, ITG *nzs, ITG *isolver, ITG *iperturb)
Definition: preiter.c:23
subroutine mafilldm(co, nk, kon, ipkon, lakon, ne, nodeboun, ndirboun, xboun, nboun, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, ad, au, nactdof, icol, jq, irow, neq, nzl, nmethod, ikmpc, ilmpc, ikboun, ilboun, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, vold, iperturb, sti, nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, ttime, time, istep, iinc, ibody, clearini, mortar, springarea, pslavsurf, pmastsurf, reltime, nasym)
Definition: mafilldm.f:19
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
subroutine precfd(ne, ipkon, kon, lakon, ipnei, neifa, neiel, ipoface, nodface, ielfa, nkonnei, nface, ifaext, nfaext, isolidsurf, nsolidsurf, set, nset, istartset, iendset, ialset, vel, vold, mi, neij, nef, nactdoh, ipkonf, lakonf, ielmatf, ielmat, ielorienf, ielorien, norien)
Definition: precfd.f:19
subroutine interpolatestate(ne, ipkon, kon, lakon, ne0, mi, xstate, pslavsurf, nstate_, xstateini, islavsurf, islavsurfold, pslavsurfold)
Definition: interpolatestate.f:28
ITG strcpy1(char *s1, const char *s2, ITG length)
Definition: strcpy1.c:24
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
void frdcyc(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *v, double *stn, ITG *inum, ITG *nmethod, ITG *kode, char *filab, double *een, double *t1, double *fn, double *time, double *epn, ITG *ielmat, char *matname, double *cs, ITG *mcs, ITG *nkon, double *enern, double *xstaten, ITG *nstate_, ITG *istep, ITG *iinc, ITG *iperturb, double *ener, ITG *mi, char *output, ITG *ithermal, double *qfn, ITG *ialset, ITG *istartset, ITG *iendset, double *trab, ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien, double *sti, double *veold, ITG *noddiam, char *set, ITG *nset, double *emn, double *thicke, char *jobnamec, ITG *ne0, double *cdn, ITG *mortar, ITG *nmat)
Definition: frdcyc.c:24
void mastructrad(ITG *ntr, ITG *nloadtr, char *sideload, ITG *ipointerrad, ITG **mast1radp, ITG **irowradp, ITG *nzsrad, ITG *jqrad, ITG *icolrad)
Definition: mastructrad.c:24
void precontact(ITG *ncont, ITG *ntie, char *tieset, ITG *nset, char *set, ITG *istartset, ITG *iendset, ITG *ialset, ITG *itietri, char *lakon, ITG *ipkon, ITG *kon, ITG *koncont, ITG *ne, double *cg, double *straight, double *co, double *vold, ITG *istep, ITG *iinc, ITG *iit, ITG *itiefac, ITG *islavsurf, ITG *islavnode, ITG *imastnode, ITG *nslavnode, ITG *nmastnode, ITG *imastop, ITG *mi, ITG *ipe, ITG *ime, double *tietol, ITG *iflagact, ITG *nintpoint, double **pslavsurfp, double *xmastnor, double *cs, ITG *mcs, ITG *ics, double *clearini, ITG *nslavs)
Definition: precontact.c:24
void spooles(double *ad, double *au, double *adb, double *aub, double *sigma, double *b, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG *symmtryflag, ITG *inputformat, ITG *nzs3)
subroutine stop()
Definition: stop.f:19
void tau(double *ad, double **aup, double *adb, double *aubp, double *sigma, double *b, ITG *icol, ITG **irowp, ITG *neq, ITG *nzs)
void inicont(ITG *nk, ITG *ncont, ITG *ntie, char *tieset, ITG *nset, char *set, ITG *istartset, ITG *iendset, ITG *ialset, ITG **itietrip, char *lakon, ITG *ipkon, ITG *kon, ITG **koncontp, ITG *ncone, double *tietol, ITG *ismallsliding, ITG **itiefacp, ITG **islavsurfp, ITG **islavnodep, ITG **imastnodep, ITG **nslavnodep, ITG **nmastnodep, ITG *mortar, ITG **imastopp, ITG *nkon, ITG **iponoels, ITG **inoelsp, ITG **ipep, ITG **imep, ITG *ne, ITG *ifacecount, ITG *nmpc, ITG *mpcfree, ITG *memmpc_, ITG **ipompcp, char **labmpcp, ITG **ikmpcp, ITG **ilmpcp, double **fmpcp, ITG **nodempcp, double **coefmpcp, ITG *iperturb, ITG *ikboun, ITG *nboun, double *co, ITG *istep, double **xnoelsp)
Definition: inicont.c:24
void radcyc(ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *cs, ITG *mcs, ITG *nkon, ITG *ialset, ITG *istartset, ITG *iendset, ITG **kontrip, ITG *ntri, double **cop, double **voldp, ITG *ntrit, ITG *inocs, ITG *mi)
Definition: radcyc.c:24
void sgi_main(double *ad, double *au, double *adb, double *aub, double *sigma, double *b, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG token)
subroutine writecvg(istep, iinc, iit, ne, ne0, ram, qam, cam, uam, ithermal)
Definition: writecvg.f:19
#define RENEW(a, b, c)
Definition: CalculiX.h:40
#define SFREE(a)
Definition: CalculiX.h:41
subroutine gennactdofinv(nactdof, nactdofinv, nk, mi, nodorig, ipkon, lakon, kon, ne)
Definition: gennactdofinv.f:19
void results(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *v, double *stn, ITG *inum, double *stx, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *alcon, ITG *nalcon, double *alzero, ITG *ielmat, ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_, double *t0, double *t1, ITG *ithermal, double *prestr, ITG *iprestr, char *filab, double *eme, double *emn, double *een, ITG *iperturb, double *f, double *fn, ITG *nactdof, ITG *iout, double *qa, double *vold, double *b, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, ITG *nmpc, ITG *nmethod, double *vmax, ITG *neq, double *veold, double *accold, double *beta, double *gamma, double *dtime, double *time, double *ttime, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *xstateini, double *xstiff, double *xstate, ITG *npmat_, double *epl, char *matname, ITG *mi, ITG *ielas, ITG *icmd, ITG *ncmat_, ITG *nstate_, double *stiini, double *vini, ITG *ikboun, ITG *ilboun, double *ener, double *enern, double *emeini, double *xstaten, double *eei, double *enerini, double *cocon, ITG *ncocon, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, double *qfx, double *qfn, double *trab, ITG *inotr, ITG *ntrans, double *fmpc, ITG *nelemload, ITG *nload, ITG *ikmpc, ITG *ilmpc, ITG *istep, ITG *iinc, double *springarea, double *reltime, ITG *ne0, double *xforc, ITG *nforc, double *thicke, double *shcon, ITG *nshcon, char *sideload, double *xload, double *xloadold, ITG *icfd, ITG *inomat, double *pslavsurf, double *pmastsurf, ITG *mortar, ITG *islavact, double *cdn, ITG *islavnode, ITG *nslavnode, ITG *ntie, double *clearini, ITG *islavsurf, ITG *ielprop, double *prop)
Definition: results.c:42
subroutine islavactive(tieset, ntie, itietri, cg, straight, co, vold, xo, yo, zo, x, y, z, nx, ny, nz, mi, imastop, nslavnode, islavnode, islavact)
Definition: islavactive.f:22
subroutine tempload(xforcold, xforc, xforcact, iamforc, nforc, xloadold, xload, xloadact, iamload, nload, ibody, xbody, nbody, xbodyold, xbodyact, t1old, t1, t1act, iamt1, nk, amta, namta, nam, ampli, time, reltime, ttime, dtime, ithermal, nmethod, xbounold, xboun, xbounact, iamboun, nboun, nodeboun, ndirboun, nodeforc, ndirforc, istep, iinc, co, vold, itg, ntg, amname, ikboun, ilboun, nelemload, sideload, mi, ntrans, trab, inotr, veold, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nmpc, ipompc, ikmpc, ilmpc, nodempc, coefmpc)
Definition: tempload.f:19
subroutine checktime(itpamp, namta, tinc, ttime, amta, tmin, inext, itp, istep)
Definition: checktime.f:19
static double * adview
Definition: radflowload.c:42
subroutine nident(x, px, n, id)
Definition: nident.f:25
subroutine mafillsm(co, nk, kon, ipkon, lakon, ne, nodeboun, ndirboun, xboun, nboun, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, ad, au, fext, nactdof, icol, jq, irow, neq, nzl, nmethod, ikmpc, ilmpc, ikboun, ilboun, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, vold, iperturb, sti, nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, mass, stiffness, buckling, rhsi, intscheme, physcon, shcon, nshcon, cocon, ncocon, ttime, time, istep, iinc, coriolis, ibody, xloadold, reltime, veold, springarea, nstate_, xstateini, xstate, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, pslavsurf, pmastsurf, mortar, clearini, ielprop, prop)
Definition: mafillsm.f:19
void getglobalresults(char *jobnamec, ITG **integerglobp, double **doubleglobp, ITG *nboun, ITG *iamboun, double *xboun, ITG *nload, char *sideload, ITG *iamload, ITG *iglob, ITG *nforc, ITG *iamforc, double *xforc, ITG *ithermal, ITG *nk, double *t1, ITG *iamt1)
void prediction(double *uam, ITG *nmethod, double *bet, double *gam, double *dtime, ITG *ithermal, ITG *nk, double *veold, double *accold, double *v, ITG *iinc, ITG *idiscon, double *vold, ITG *nactdof, ITG *mi)
Definition: prediction.c:33
real *8 function f_cm(x, phi, lambda1, zk0, Pup, Tup, rurd, xflow, kup)
Definition: moehring.f:579
void contact(ITG *ncont, ITG *ntie, char *tieset, ITG *nset, char *set, ITG *istartset, ITG *iendset, ITG *ialset, ITG *itietri, char *lakon, ITG *ipkon, ITG *kon, ITG *koncont, ITG *ne, double *cg, double *straight, ITG *ifree, double *co, double *vold, ITG *ielmat, double *cs, double *elcon, ITG *istep, ITG *iinc, ITG *iit, ITG *ncmat_, ITG *ntmat_, ITG *ne0, double *vini, ITG *nmethod, ITG *nmpc, ITG *mpcfree, ITG *memmpc_, ITG **ipompcp, char **labmpcp, ITG **ikmpcp, ITG **ilmpcp, double **fmpcp, ITG **nodempcp, double **coefmpcp, ITG *iperturb, ITG *ikboun, ITG *nboun, ITG *mi, ITG *imastop, ITG *nslavnode, ITG *islavnode, ITG *islavsurf, ITG *itiefac, double *areaslav, ITG *iponoels, ITG *inoels, double *springarea, double *tietol, double *reltime, ITG *imastnode, ITG *nmastnode, double *xmastnor, char *filab, ITG *mcs, ITG *ics, ITG *nasym, double *xnoels, ITG *mortar, double *pslavsurf, double *pmastsurf, double *clearini, double *theta)
Definition: contact.c:23
subroutine nonlinmpc(co, vold, ipompc, nodempc, coefmpc, labmpc, nmpc, ikboun, ilboun, nboun, xbounact, aux, iaux, maxlenmpc, ikmpc, ilmpc, icascade, kon, ipkon, lakon, ne, reltime, newstep, xboun, fmpc, iit, idiscon, ncont, trab, ntrans, ithermal, mi)
Definition: nonlinmpc.f:19
void checkconvergence(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, double *stn, ITG *nmethod, ITG *kode, char *filab, double *een, double *t1act, double *time, double *epn, ITG *ielmat, char *matname, double *enern, double *xstaten, ITG *nstate_, ITG *istep, ITG *iinc, ITG *iperturb, double *ener, ITG *mi, char *output, ITG *ithermal, double *qfn, ITG *mode, ITG *noddiam, double *trab, ITG *inotr, ITG *ntrans, double *orab, ITG *ielorien, ITG *norien, char *description, double *sti, ITG *icutb, ITG *iit, double *dtime, double *qa, double *vold, double *qam, double *ram1, double *ram2, double *ram, double *cam, double *uam, ITG *ntg, double *ttime, ITG *icntrl, double *theta, double *dtheta, double *veold, double *vini, ITG *idrct, double *tper, ITG *istab, double *tmax, ITG *nactdof, double *b, double *tmin, double *ctrl, double *amta, ITG *namta, ITG *itpamp, ITG *inext, double *dthetaref, ITG *itp, ITG *jprint, ITG *jout, ITG *uncoupled, double *t1, ITG *iitterm, ITG *nelemload, ITG *nload, ITG *nodeboun, ITG *nboun, ITG *itg, ITG *ndirboun, double *deltmx, ITG *iflagact, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, double *emn, double *thicke, char *jobnamec, ITG *mortar, ITG *nmat, ITG *ielprop, double *prop)
void calcresidual(ITG *nmethod, ITG *neq, double *b, double *fext, double *f, ITG *iexpl, ITG *nactdof, double *aux2, double *vold, double *vini, double *dtime, double *accold, ITG *nk, double *adb, double *aub, ITG *jq, ITG *irow, ITG *nzl, double *alpha, double *fextini, double *fini, ITG *islavnode, ITG *nslavnode, ITG *mortar, ITG *ntie, double *f_cm, double *f_cs, ITG *mi, ITG *nzs, ITG *nasym, ITG *idamping, double *veold, double *adc, double *auc, double *cvini, double *cv)
Definition: calcresidual.c:33
subroutine gasmechbc(vold, nload, sideload, nelemload, xload, mi)
Definition: gasmechbc.f:19
subroutine envtemp(itg, ieg, ntg, ntr, sideload, nelemload, ipkon, kon, lakon, ielmat, ne, nload, kontri, ntri, nloadtr, nflow, ndirboun, nactdog, nodeboun, nacteq, nboun, ielprop, prop, nteq, v, network, physcon, shcon, ntmat_, co, vold, set, nshcon, rhcon, nrhcon, mi, nmpc, nodempc, ipompc, labmpc, ikboun, nasym, iaxial)
Definition: envtemp.f:19
void radflowload(ITG *itg, ITG *ieg, ITG *ntg, ITG *ntr, double *adrad, double *aurad, double *bcr, ITG *ipivr, double *ac, double *bc, ITG *nload, char *sideload, ITG *nelemload, double *xloadact, char *lakon, ITG *ipiv, ITG *ntmat_, double *vold, double *shcon, ITG *nshcon, ITG *ipkon, ITG *kon, double *co, ITG *kontri, ITG *ntri, ITG *nloadtr, double *tarea, double *tenv, double *physcon, double *erad, double **adviewp, double **auviewp, ITG *nflow, ITG *ikboun, double *xboun, ITG *nboun, ITG *ithermal, ITG *iinc, ITG *iit, double *cs, ITG *mcs, ITG *inocs, ITG *ntrit, ITG *nk, double *fenv, ITG *istep, double *dtime, double *ttime, double *time, ITG *ilboun, ITG *ikforc, ITG *ilforc, double *xforc, ITG *nforc, double *cam, ITG *ielmat, ITG *nteq, double *prop, ITG *ielprop, ITG *nactdog, ITG *nacteq, ITG *nodeboun, ITG *ndirboun, ITG *network, double *rhcon, ITG *nrhcon, ITG *ipobody, ITG *ibody, double *xbody, ITG *nbody, ITG *iviewfile, char *jobnamef, double *ctrl, double *xloadold, double *reltime, ITG *nmethod, char *set, ITG *mi, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nset, ITG *ineighe, ITG *nmpc, ITG *nodempc, ITG *ipompc, double *coefmpc, char *labmpc, ITG *iemchange, ITG *nam, ITG *iamload, ITG *jqrad, ITG *irowrad, ITG *nzsrad, ITG *icolrad, ITG *ne, ITG *iaxial)
Definition: radflowload.c:45
void compfluid(double **cop, ITG *nk, ITG **ipkonp, ITG **konp, char **lakonp, char **sideface, ITG *ifreestream, ITG *nfreestream, ITG *isolidsurf, ITG *neighsolidsurf, ITG *nsolidsurf, ITG **iponoel, ITG **inoel, ITG *nshcon, double *shcon, ITG *nrhcon, double *rhcon, double **voldp, ITG *ntmat_, ITG *nodeboun, ITG *ndirboun, ITG *nboun, ITG **ipompcp, ITG **nodempcp, ITG *nmpc, ITG **ikmpcp, ITG **ilmpcp, ITG *ithermal, ITG *ikboun, ITG *ilboun, ITG *turbulent, ITG *isolver, ITG *iexpl, double *vcontu, double *ttime, double *time, double *dtime, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, double *xbody, ITG *ipobody, ITG *nbody, ITG *ielmat, char *matname, ITG *mi, ITG *ncmat_, double *physcon, ITG *istep, ITG *iinc, ITG *ibody, double *xloadold, double *xboun, double **coefmpcp, ITG *nmethod, double *xforcold, double *xforcact, ITG *iamforc, ITG *iamload, double *xbodyold, double *xbodyact, double *t1old, double *t1, double *t1act, ITG *iamt1, double *amta, ITG *namta, ITG *nam, double *ampli, double *xbounold, double *xbounact, ITG *iamboun, ITG *itg, ITG *ntg, char *amname, double *t0, ITG **nelemface, ITG *nface, double *cocon, ITG *ncocon, double *xloadact, double *tper, ITG *jmax, ITG *jout, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, char *prset, char *prlab, ITG *nprint, double *trab, ITG *inotr, ITG *ntrans, char *filab, char **labmpcp, double *sti, ITG *norien, double *orab, char *jobnamef, char *tieset, ITG *ntie, ITG *mcs, ITG *ics, double *cs, ITG *nkon, ITG *mpcfree, ITG *memmpc_, double **fmpcp, ITG *nef, ITG **inomat, double *qfx, ITG *neifa, ITG *neiel, ITG *ielfa, ITG *ifaext, double *vfa, double *vel, ITG *ipnei, ITG *nflnei, ITG *nfaext, char *typeboun, ITG *neij, double *tincf, ITG *nactdoh, ITG *nactdohinv, ITG *ielorien)
Definition: compfluid.c:39
#define ITG
Definition: CalculiX.h:51
#define NNEW(a, b, c)
Definition: CalculiX.h:39
subroutine bodyforce(cbody, ibody, ipobody, nbody, set, istartset, iendset, ialset, inewton, nset, ifreebody, k)
Definition: bodyforce.f:19
subroutine rhs(co, nk, kon, ipkon, lakon, ne, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, fext, nactdof, neq, nmethod, ikmpc, ilmpc, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, iprestr, vold, iperturb, iexpl, plicon, nplicon, plkcon, nplkcon, npmat_, ttime, time, istep, iinc, dtime, physcon, ibody, xloadold, reltime, veold, matname, mi, ikactmech, nactmech, ielprop, prop)
Definition: rhs.f:19
static double * auview
Definition: radflowload.c:42
Hosted by OpenAircraft.com, (Michigan UAV, LLC)