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

Go to the source code of this file.

Functions

void dyna (double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, ITG **nodebounp, ITG **ndirbounp, double **xbounp, 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 **nactdofp, ITG *neq, ITG *nzl, ITG *icol, ITG *irow, ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG **ikbounp, ITG **ilbounp, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *cocon, ITG *ncocon, double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp, ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_, double **t0p, double **t1p, ITG *ithermal, double *prestr, ITG *iprestr, double **voldp, ITG *iperturb, double **stip, ITG *nzs, double *timepar, double *xmodal, double **veoldp, char *amname, double *amta, ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload, ITG **iamt1p, ITG *jout, ITG *kode, char *filab, double **emep, double *xforcold, double *xloadold, double **t1oldp, ITG **iambounp, double **xbounoldp, ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double **xstatep, ITG *npmat_, char *matname, ITG *mi, ITG *ncmat_, ITG *nstate_, double **enerp, char *jobnamec, double *ttime, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG **ialsetp, ITG *nprint, char *prlab, char *prset, ITG *nener, double *trab, ITG **inotrp, ITG *ntrans, double **fmpcp, char *cbody, ITG *ibody, double *xbody, ITG *nbody, double *xbodyold, ITG *istep, ITG *isolver, ITG *jq, char *output, ITG *mcs, ITG *nkon, ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset, ITG *idrct, ITG *jmax, double *ctrl, ITG *itpamp, double *tietol, ITG *nalset, ITG *ikforc, ITG *ilforc, double *thicke, ITG *nslavs, ITG *nmat, char *typeboun, ITG *ielprop, double *prop)
 

Function Documentation

void dyna ( double **  cop,
ITG nk,
ITG **  konp,
ITG **  ipkonp,
char **  lakonp,
ITG ne,
ITG **  nodebounp,
ITG **  ndirbounp,
double **  xbounp,
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 **  nactdofp,
ITG neq,
ITG nzl,
ITG icol,
ITG irow,
ITG nmethod,
ITG **  ikmpcp,
ITG **  ilmpcp,
ITG **  ikbounp,
ITG **  ilbounp,
double *  elcon,
ITG nelcon,
double *  rhcon,
ITG nrhcon,
double *  cocon,
ITG ncocon,
double *  alcon,
ITG nalcon,
double *  alzero,
ITG **  ielmatp,
ITG **  ielorienp,
ITG norien,
double *  orab,
ITG ntmat_,
double **  t0p,
double **  t1p,
ITG ithermal,
double *  prestr,
ITG iprestr,
double **  voldp,
ITG iperturb,
double **  stip,
ITG nzs,
double *  timepar,
double *  xmodal,
double **  veoldp,
char *  amname,
double *  amta,
ITG namta,
ITG nam,
ITG iamforc,
ITG iamload,
ITG **  iamt1p,
ITG jout,
ITG kode,
char *  filab,
double **  emep,
double *  xforcold,
double *  xloadold,
double **  t1oldp,
ITG **  iambounp,
double **  xbounoldp,
ITG iexpl,
double *  plicon,
ITG nplicon,
double *  plkcon,
ITG nplkcon,
double **  xstatep,
ITG npmat_,
char *  matname,
ITG mi,
ITG ncmat_,
ITG nstate_,
double **  enerp,
char *  jobnamec,
double *  ttime,
char *  set,
ITG nset,
ITG istartset,
ITG iendset,
ITG **  ialsetp,
ITG nprint,
char *  prlab,
char *  prset,
ITG nener,
double *  trab,
ITG **  inotrp,
ITG ntrans,
double **  fmpcp,
char *  cbody,
ITG ibody,
double *  xbody,
ITG nbody,
double *  xbodyold,
ITG istep,
ITG isolver,
ITG jq,
char *  output,
ITG mcs,
ITG nkon,
ITG mpcend,
ITG ics,
double *  cs,
ITG ntie,
char *  tieset,
ITG idrct,
ITG jmax,
double *  ctrl,
ITG itpamp,
double *  tietol,
ITG nalset,
ITG ikforc,
ITG ilforc,
double *  thicke,
ITG nslavs,
ITG nmat,
char *  typeboun,
ITG ielprop,
double *  prop 
)

Definition at line 36 of file dyna.c.

71  {
72 
73  char fneig[132]="",description[13]=" ",*lakon=NULL,*labmpc=NULL,
74  *labmpcold=NULL,lakonl[9]=" \0",*tchar1=NULL,*tchar2=NULL,
75  *tchar3=NULL,cflag[1]=" ";
76 
77  ITG nev,i,j,k,idof,*inum=NULL,*ipobody=NULL,inewton=0,id,
78  iinc=0,jprint=0,l,iout,ielas=0,icmd,iprescribedboundary,init,ifreebody,
79  mode=-1,noddiam=-1,*kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,
80  *inotr=NULL,*nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,
81  *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
82  *ilmpc=NULL,nsectors,nmpcold,mpcendold,*ipompcold=NULL,*nodempcold=NULL,
83  *ikmpcold=NULL,*ilmpcold=NULL,kflag=2,nmd,nevd,*nm=NULL,*iamt1=NULL,
84  *itg=NULL,ntg=0,symmetryflag=0,inputformat=0,dashpot,lrw,liw,iddebdf=0,
85  *iwork=NULL,ngraph=1,nkg,neg,ncont,ne0,nkon0, *itietri=NULL,
86  *koncont=NULL,konl[20],imat,nope,kodem,indexe,j1,jdof,
87  *ipneigh=NULL,*neigh=NULL,inext,itp=0,*islavact=NULL,
88  ismallsliding=0,isteadystate,mpcfree,im,cyclicsymmetry,
89  memmpc_,imax,iener=0,*icole=NULL,*irowe=NULL,*jqe=NULL,nzse[3],
90  nalset_=*nalset,*ialset=*ialsetp,*istartset_=NULL,*iendset_=NULL,
91  *itiefac=NULL,*islavsurf=NULL,*islavnode=NULL,mt=mi[1]+1,
92  *imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,mortar=0,*imastop=NULL,
93  *iponoels=NULL,*inoels=NULL,*imddof=NULL,nmddof,
94  *ikactcont=NULL,nactcont,nactcont_=100,*ikactmech=NULL,nactmech,
95  iabsload=0,*ipe=NULL,*ime=NULL,iprev=1,inonlinmpc=0,ielem,
96  *imdnode=NULL,nmdnode,*imdboun=NULL,nmdboun,*imdmpc=NULL,
97  nmdmpc,intpointvar,kmin,kmax,i1,ifacecount,*izdof=NULL,
98  nzdof,iload,iforc,*iponoel=NULL,*inoel=NULL,*imdelem=NULL,nmdelem,
99  irenewxstate,nasym=0,*nshcon=NULL,nherm,icfd=0,*inomat=NULL;
100 
101  long long i2;
102 
103  double *d=NULL, *z=NULL, *b=NULL, *zeta=NULL,*stiini=NULL,
104  *cd=NULL, *cv=NULL, *xforcact=NULL, *xloadact=NULL,*cc=NULL,
105  *t1act=NULL, *ampli=NULL, *aa=NULL, *bb=NULL, *aanew=NULL, *bj=NULL,
106  *v=NULL,*aamech=NULL,*emn=NULL,*cdn=NULL,ptime,
107  *stn=NULL, *stx=NULL, *een=NULL, *adb=NULL,*xstiff=NULL,*bjp=NULL,
108  *aub=NULL, *temp_array1=NULL, *temp_array2=NULL, *aux=NULL,
109  *f=NULL, *fn=NULL, *xbounact=NULL,*epn=NULL,*xstateini=NULL,
110  *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
111  *qfx=NULL, *xbodyact=NULL, *cgr=NULL, *au=NULL, *vbounact=NULL,
112  *abounact=NULL,dtime,reltime,*t0=NULL,*t1=NULL,*t1old=NULL,
113  physcon[1],zetaj,dj,ddj,h1,h2,h3,h4,h5,h6,sum,aai,bbi,tstart,tend,
114  qa[3],cam[5],accold[1],bet,gam,*ad=NULL,sigma=0.,alpham,betam,
115  *bact=NULL,*bmin=NULL,*co=NULL,*xboun=NULL,*xbounold=NULL,*vold=NULL,
116  *eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,*coefmpcold,*veold=NULL,
117  *xini=NULL,*rwork=NULL,*adc=NULL,*auc=NULL,*zc=NULL, *rpar=NULL,
118  *cg=NULL,*straight=NULL,xl[27],voldl[mt*9],elas[21],fnl[27],t0l,t1l,
119  elconloc[21],veoldl[mt*9],setnull,deltmx,bbmax,dd,dtheta,dthetaref,
120  theta,*vini=NULL,dthetaold,*bcont=NULL,*vr=NULL,*vi=NULL,*bcontini=NULL,
121  *stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,precision,resultmaxprev,
122  resultmax,func,funcp,fexp,fexm,fcos,fsin,sump,*bp=NULL,h14,senergy=0.0,
123  *bv=NULL,*cstr=NULL,*aube=NULL,*adbe=NULL,*sti=*stip,time0=0.0,
124  time=0.0,*xforcdiff=NULL,*xloaddiff=NULL,*xbodydiff=NULL,*t1diff=NULL,
125  *xboundiff=NULL,*bprev=NULL,*bdiff=NULL,*areaslav=NULL,
126  *springarea=NULL, *bold=NULL,*eenmax=NULL,*fnr=NULL,*fni=NULL,
127  *xmastnor=NULL,*emeini=NULL,*xstate=NULL,*clearini=NULL,
128  *shcon=NULL,*xmr=NULL,*xmi=NULL,*xnoels=NULL,*pslavsurf=NULL,
129  *pmastsurf=NULL,*cdnr=NULL,*cdni=NULL,*tinc,*tper,*tmin,*tmax;
130 
131  FILE *f1;
132 
133  /* dummy variables for nonlinmpc */
134 
135  ITG *iaux=NULL,maxlenmpc,icascade=0,newstep=0,iit=-1,idiscon;
136 
137 #ifdef SGI
138  ITG token;
139 #endif
140 
141  /* if iabsload=0: aamech is modified by the present incremental
142  contribution of b
143  iabsload=1: the last incremental contribution is
144  subtracted before the new one is added to b;
145  this latter incremental contribution is used
146  to update aamech
147  iabsload=2: aamech is determined by the absolute
148  contribution of b (no incremental procedure
149  for the load; this is necessary if
150  - nonlinear MPC's are applied or
151  - user dloads are applied */
152 
153  co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
154  ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
155  ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;xstate=*xstatep;
156  xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
157  vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
158  coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
159  fmpc=*fmpcp;veold=*veoldp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
160 
161  tinc=&timepar[0];
162  tper=&timepar[1];
163  tmin=&timepar[2];
164  tmax=&timepar[3];
165 
166  if(ithermal[0]<=1){
167  kmin=1;kmax=3;
168  }else if(ithermal[0]==2){
169  kmin=0;kmax=mi[1];if(kmax>2)kmax=2;
170  }else{
171  kmin=0;kmax=3;
172  }
173 
174 //xstiffNNEW(// xstiff,double,(long long)27*mi[0]**ne);
175 
176  dtime=*tinc;
177 
178  alpham=xmodal[0];
179  betam=xmodal[1];
180 
181  dd=ctrl[16];deltmx=ctrl[26];
182 
183  /* determining nzl */
184 
185  *nzl=0;
186  for(i=neq[1];i>0;i--){
187  if(icol[i-1]>0){
188  *nzl=i;
189  break;
190  }
191  }
192 
193  /* opening the eigenvalue file and checking for cyclic symmetry */
194 
195  strcpy(fneig,jobnamec);
196  strcat(fneig,".eig");
197 
198  if((f1=fopen(fneig,"rb"))==NULL){
199  printf("*ERROR in dyna: cannot open eigenvalue file for reading");
200  exit(0);
201  }
202 
203  printf(" *INFO in dyna: if there are problems reading the .eig file this may be due to:\n");
204  printf(" 1) the nonexistence of the .eig file\n");
205  printf(" 2) other boundary conditions than in the input deck\n");
206  printf(" which created the .eig file\n\n");
207 
208  if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
209  printf("*ERROR in dyna reading the cyclic symmetry flag in the eigenvalue file");
210  exit(0);
211  }
212 
213  if(fread(&nherm,sizeof(ITG),1,f1)!=1){
214  printf("*ERROR in dyna reading the Hermitian flag in the eigenvalue file");
215  exit(0);
216  }
217 
218  if(nherm!=1){
219  printf("*ERROR in dyna: the eigenvectors in the .eig-file result\n");
220  printf(" from a non-Hermitian eigenvalue problem. The modal\n");
221  printf(" dynamic procedure cannot handle that yet\n\n");
222  FORTRAN(stop,());
223  }
224 
225  /* creating imddof containing the degrees of freedom
226  retained by the user and imdnode containing the nodes */
227 
228  nmddof=0;nmdnode=0;nmdboun=0;nmdmpc=0;nmdelem=0;
229 
230  NNEW(imddof,ITG,*nk*3);
231  NNEW(imdnode,ITG,*nk);
232  NNEW(imdboun,ITG,*nboun);
233  NNEW(imdmpc,ITG,*nmpc);
234  FORTRAN(createmddof,(imddof,&nmddof,istartset,iendset,
235  ialset,nactdof,ithermal,mi,imdnode,&nmdnode,
236  ikmpc,ilmpc,ipompc,nodempc,nmpc,
237  imdmpc,&nmdmpc,imdboun,&nmdboun,ikboun,
238  nboun,nset,ntie,tieset,set,lakon,kon,ipkon,labmpc,
239  ilboun,filab,prlab,prset,nprint,ne,&cyclicsymmetry));
240 
241  /* if results are requested in too many nodes, it is faster to
242  calculate the results in all nodes */
243 
244  if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
245  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
246  }
247 
248  if(nmdnode!=0){
249  if(!cyclicsymmetry){
250  for(i=0;i<*nload;i++){
251  iload=i+1;
252  FORTRAN(addimdnodedload,(nelemload,sideload,ipkon,kon,lakon,
253  &iload,imdnode,&nmdnode,ikmpc,ilmpc,ipompc,nodempc,nmpc,
254  imddof,&nmddof,nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
255  ikboun,nboun,ilboun,ithermal));
256  }
257  for(i=0;i<*nforc;i++){
258  iforc=i+1;
259  FORTRAN(addimdnodecload,(nodeforc,&iforc,imdnode,&nmdnode,xforc,
260  ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
261  nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
262  ikboun,nboun,ilboun,ithermal));
263  }
264  }
265 
266  /* determining the elements belonging to a given node */
267 
268  NNEW(iponoel,ITG,*nk);
269  NNEW(inoel,ITG,2**nkon);
270  FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
271  NNEW(imdelem,ITG,*ne);
272 
273  /* storing the elements in which integration point results
274  are needed; storing the nodes which are needed to
275  calculate these results */
276 
277  FORTRAN(createmdelem,(imdnode,&nmdnode,xforc,
278  ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
279  nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
280  ikboun,nboun,ilboun,ithermal,imdelem,&nmdelem,
281  iponoel,inoel,prlab,prset,nprint,lakon,set,nset,
282  ialset,ipkon,kon,istartset,iendset,nforc,
283  ikforc,ilforc));
284 
285  RENEW(imdelem,ITG,nmdelem);
286  SFREE(iponoel);SFREE(inoel);
287  }
288 
289  /* if results are requested in too many nodes, it is faster to
290  calculate the results in all nodes */
291 
292  if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
293  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
294  }
295 
296  /* subtracting 1 to comply with the C-convention */
297 
298  for(i=0;i<nmddof;i++){imddof[i]-=1;}
299  RENEW(imddof,ITG,nmddof);
300  RENEW(imdnode,ITG,nmdnode);
301  RENEW(imdboun,ITG,nmdboun);
302  RENEW(imdmpc,ITG,nmdmpc);
303 
304  nsectors=1;
305 
306  /* reading the eigenvalues / eigenmodes */
307 
308  if(!cyclicsymmetry){
309 
310  nkg=*nk;
311  neg=*ne;
312 
313  if(fread(&nev,sizeof(ITG),1,f1)!=1){
314  printf("*ERROR in dyna reading the number of eigenvalues in the eigenvalue file");
315  exit(0);
316  }
317 
318  NNEW(d,double,nev);
319 
320  if(fread(d,sizeof(double),nev,f1)!=nev){
321  printf("*ERROR in dyna reading the eigenvalues in the eigenvalue file");
322  exit(0);
323  }
324 
325  for(i=0;i<nev;i++){
326  if(d[i]>0){d[i]=sqrt(d[i]);}else{d[i]=0.;}
327  }
328 
329  NNEW(ad,double,neq[1]);
330  NNEW(adb,double,neq[1]);
331  NNEW(au,double,nzs[2]);
332  NNEW(aub,double,nzs[1]);
333 
334  if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
335  printf("*ERROR in dyna reading the diagonal of the stiffness matrix in the eigenvalue file");
336  exit(0);
337  }
338 
339  if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
340  printf("*ERROR in dyna reading the off-diagonals of the stiffness matrix in the eigenvalue file");
341  exit(0);
342  }
343 
344  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
345  printf("*ERROR in dyna reading the diagonal of the mass matrix in the eigenvalue file");
346  exit(0);
347  }
348 
349  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
350  printf("*ERROR in dyna reading the off-diagonals of the mass matrix in the eigenvalue file");
351  exit(0);
352  }
353 
354  NNEW(z,double,neq[1]*nev);
355 
356  if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
357  printf("*ERROR in dyna reading the eigenvectors in the eigenvalue file");
358  exit(0);
359  }
360  }
361  else{
362  nev=0;
363  do{
364  if(fread(&nmd,sizeof(ITG),1,f1)!=1){
365  break;
366  }
367  if(fread(&nevd,sizeof(ITG),1,f1)!=1){
368  printf("*ERROR in dyna reading the number of eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
369  exit(0);
370  }
371  if(nev==0){
372  NNEW(d,double,nevd);
373  NNEW(nm,ITG,nevd);
374  }else{
375  RENEW(d,double,nev+nevd);
376  RENEW(nm,ITG,nev+nevd);
377  }
378 
379  if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
380  printf("*ERROR in dyna reading the eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
381  exit(0);
382  }
383 
384  for(i=nev;i<nev+nevd;i++){
385  if(d[i]>0){d[i]=sqrt(d[i]);}else{d[i]=0.;}
386  }
387 
388  for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
389 
390  if(nev==0){
391  NNEW(adb,double,neq[1]);
392  NNEW(aub,double,nzs[1]);
393 
394  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
395  printf("*ERROR in dyna reading the diagonal of the mass matrix in the eigenvalue file");
396  exit(0);
397  }
398 
399  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
400  printf("*ERROR in dyna reading the off-diagonals of the mass matrix in the eigenvalue file");
401  exit(0);
402  }
403  }
404 
405  if(nev==0){
406  NNEW(z,double,neq[1]*nevd);
407  }else{
408  RENEW(z,double,(long long)neq[1]*(nev+nevd));
409  }
410 
411  if(fread(&z[(long long)neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
412  printf("*ERROR in dyna reading the eigenvectors for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
413  exit(0);
414  }
415  nev+=nevd;
416  }while(1);
417 
418  /* determining the maximum amount of segments */
419 
420  for(i=0;i<*mcs;i++){
421 // if(cs[17*i]>nsectors) nsectors=cs[17*i];
422  if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
423  }
424 
425  /* determining the maximum number of sectors to be plotted */
426 
427  for(j=0;j<*mcs;j++){
428  if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
429  }
430  nkg=*nk*ngraph;
431  neg=*ne*ngraph;
432 
433  /* allocating field for the expanded structure */
434 
435  RENEW(co,double,3**nk*nsectors);
436 
437  /* next line is necessary for multiple cyclic symmetry
438  conditions */
439 
440  for(i=3**nk;i<3**nk*nsectors;i++){co[i]=0.;}
441  if(*ithermal!=0){
442  RENEW(t0,double,*nk*nsectors);
443  RENEW(t1old,double,*nk*nsectors);
444  RENEW(t1,double,*nk*nsectors);
445  if(*nam>0) RENEW(iamt1,ITG,*nk*nsectors);
446  }
447  RENEW(nactdof,ITG,mt**nk*nsectors);
448  if(*ntrans>0) RENEW(inotr,ITG,2**nk*nsectors);
449  RENEW(kon,ITG,*nkon*nsectors);
450  RENEW(ipkon,ITG,*ne*nsectors);
451  for(i=*ne;i<*ne*nsectors;i++){ipkon[i]=-1;}
452  RENEW(lakon,char,8**ne*nsectors);
453  RENEW(ielmat,ITG,mi[2]**ne*nsectors);
454  if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne*nsectors);
455 // RENEW(z,double,(long long)neq[1]*nev*nsectors/2);
456 
457  RENEW(nodeboun,ITG,*nboun*nsectors);
458  RENEW(ndirboun,ITG,*nboun*nsectors);
459  if(*nam>0) RENEW(iamboun,ITG,*nboun*nsectors);
460  RENEW(xboun,double,*nboun*nsectors);
461  RENEW(xbounold,double,*nboun*nsectors);
462  RENEW(ikboun,ITG,*nboun*nsectors);
463  RENEW(ilboun,ITG,*nboun*nsectors);
464 
465  NNEW(ipompcold,ITG,*nmpc);
466  NNEW(nodempcold,ITG,3**mpcend);
467  NNEW(coefmpcold,double,*mpcend);
468  NNEW(labmpcold,char,20**nmpc);
469  NNEW(ikmpcold,ITG,*nmpc);
470  NNEW(ilmpcold,ITG,*nmpc);
471 
472  for(i=0;i<*nmpc;i++){ipompcold[i]=ipompc[i];}
473  for(i=0;i<3**mpcend;i++){nodempcold[i]=nodempc[i];}
474  for(i=0;i<*mpcend;i++){coefmpcold[i]=coefmpc[i];}
475  for(i=0;i<20**nmpc;i++){labmpcold[i]=labmpc[i];}
476  for(i=0;i<*nmpc;i++){ikmpcold[i]=ikmpc[i];}
477  for(i=0;i<*nmpc;i++){ilmpcold[i]=ilmpc[i];}
478  nmpcold=*nmpc;
479  mpcendold=*mpcend;
480 
481  RENEW(ipompc,ITG,*nmpc*nsectors);
482  RENEW(nodempc,ITG,3**mpcend*nsectors);
483  RENEW(coefmpc,double,*mpcend*nsectors);
484  RENEW(labmpc,char,20**nmpc*nsectors+1);
485  RENEW(ikmpc,ITG,*nmpc*nsectors);
486  RENEW(ilmpc,ITG,*nmpc*nsectors);
487  RENEW(fmpc,double,*nmpc*nsectors);
488 
489  /* determining the space needed to expand the
490  contact surfaces */
491 
492  NNEW(tchar1,char,81);
493  NNEW(tchar2,char,81);
494  NNEW(tchar3,char,81);
495  for(i=0; i<*ntie; i++){
496  if(tieset[i*(81*3)+80]=='C'){
497 
498  //a contact constraint was found, so increase nalset
499 
500  memcpy(tchar2,&tieset[i*(81*3)+81],81);
501  tchar2[80]='\0';
502  memcpy(tchar3,&tieset[i*(81*3)+81+81],81);
503  tchar3[80]='\0';
504  for(j=0; j<*nset; j++){
505  memcpy(tchar1,&set[j*81],81);
506  tchar1[80]='\0';
507  if(strcmp(tchar1,tchar2)==0){
508 
509  //dependent nodal surface was found
510 
511  (*nalset)+=(iendset[j]-istartset[j]+1)*(nsectors);
512  }
513  else if(strcmp(tchar1,tchar3)==0){
514 
515  //independent element face surface was found
516 
517  (*nalset)+=(iendset[j]-istartset[j]+1)*(nsectors);
518  }
519  }
520  }
521  }
522  SFREE(tchar1);
523  SFREE(tchar2);
524  SFREE(tchar3);
525 
526  RENEW(ialset,ITG,*nalset);
527 
528  /* save the information in istarset and isendset */
529 
530  NNEW(istartset_,ITG,*nset);
531  NNEW(iendset_,ITG,*nset);
532  for(j=0; j<*nset; j++){
533  istartset_[j]=istartset[j];
534  iendset_[j]=iendset[j];
535  }
536 
537  /* reallocating the fields for the nodes in which the
538  solution has to be calculated */
539 
540  RENEW(imddof,ITG,neq[1]/2*nsectors);
541  RENEW(imdnode,ITG,*nk*nsectors);
542  RENEW(imdboun,ITG,*nboun*nsectors);
543  RENEW(imdmpc,ITG,*nmpc*nsectors);
544 
545 //izdofNNEW(// izdof,ITG,1);
546 
547  expand(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
548  ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
549  nforc,nelemload,sideload,xload,nload,nactdof,neq,
550  nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
551  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
552  t0,ithermal,prestr,iprestr,vold,iperturb,sti,nzs,
553  adb,aub,filab,eme,plicon,nplicon,plkcon,nplkcon,
554  xstate,npmat_,matname,mi,ics,cs,mpcend,ncmat_,
555  nstate_,mcs,nkon,ener,jobnamec,output,set,nset,istartset,
556  iendset,ialset,nprint,prlab,prset,nener,trab,
557  inotr,ntrans,ttime,fmpc,&nev,&z,iamboun,xbounold,
558  &nsectors,nm,icol,irow,nzl,nam,ipompcold,nodempcold,coefmpcold,
559  labmpcold,&nmpcold,xloadold,iamload,t1old,t1,iamt1,xstiff,
560  &icole,&jqe,&irowe,isolver,nzse,&adbe,&aube,iexpl,
561  ibody,xbody,nbody,cocon,ncocon,tieset,ntie,imddof,&nmddof,
562  imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,&izdof,&nzdof,
563  &nherm,xmr,xmi,typeboun,ielprop,prop);
564 
565  RENEW(imddof,ITG,nmddof);
566  RENEW(imdnode,ITG,nmdnode);
567  RENEW(imdboun,ITG,nmdboun);
568  RENEW(imdmpc,ITG,nmdmpc);
569 
570  SFREE(vold);
571  NNEW(vold,double,mt**nk);
572  SFREE(veold);
573  NNEW(veold,double,mt**nk);
574  RENEW(eme,double,6*mi[0]**ne);
575  RENEW(sti,double,6*mi[0]**ne);
576 
577 // RENEW(xstiff,double,(long long)27*mi[0]**ne);
578  if(*nener==1) RENEW(ener,double,mi[0]**ne*2);
579  }
580 
581  fclose(f1);
582 
583  /* checking for steadystate calculations */
584 
585  if(*tper<0){
586  precision=-*tper;
587  *tper=1.e10;
588  isteadystate=1;
589  }else{
590  isteadystate=0;
591  }
592 
593  /* checking for nonlinear MPC's */
594 
595  for(i=0;i<*nmpc;i++){
596  if((strcmp1(&labmpc[20*i]," ")!=0)&&
597  (strcmp1(&labmpc[20*i],"CONTACT")!=0)&&
598  (strcmp1(&labmpc[20*i],"CYCLIC")!=0)&&
599  (strcmp1(&labmpc[20*i],"SUBCYCLIC")!=0)){
600  inonlinmpc=1;
601  iabsload=2;
602  break;
603  }
604  }
605 
606 
607  /* normalizing the time */
608 
609  FORTRAN(checktime,(itpamp,namta,tinc,ttime,amta,tmin,&inext,&itp,istep));
610  dtheta=(*tinc)/(*tper);
611  dthetaref=dtheta;
612  dthetaold=dtheta;
613 
614  *tmin=*tmin/(*tper);
615  *tmax=*tmax/(*tper);
616  theta=0.;
617 
618  /* check for rigid body modes
619  if there is a jump of 1.e4 in two subsequent eigenvalues
620  all eigenvalues preceding the jump are considered to
621  be rigid body modes and their frequency is set to zero */
622 
623  setnull=1.;
624  for(i=nev-2;i>-1;i--){
625  if(fabs(d[i])<0.0001*fabs(d[i+1])) setnull=0.;
626  d[i]*=setnull;
627  }
628 
629  /* check whether there are dashpot elements */
630 
631  dashpot=0;
632  for(i=0;i<*ne;i++){
633  if(ipkon[i]<0) continue;
634  if(strcmp1(&lakon[i*8],"ED")==0){
635  dashpot=1;break;}
636  }
637 
638  if(dashpot){
639 
640  if(cyclicsymmetry){
641  printf("*ERROR in dyna: dashpots are not allowed in combination with cyclic symmetry\n");
642  FORTRAN(stop,());
643  }
644 
645  liw=51;
646  NNEW(iwork,ITG,liw);
647  lrw=130+42*nev;
648  NNEW(rwork,double,lrw);
649  NNEW(xini,double,2*nev);
650  NNEW(adc,double,neq[1]);
651  NNEW(auc,double,nzs[1]);
652  FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
653  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
654  nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
655  adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
656  ikmpc,ilmpc,ikboun,ilboun,
657  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
658  ielorien,norien,orab,ntmat_,
659  t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
660  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
661  xstiff,npmat_,&dtime,matname,mi,ncmat_,
662  ttime,&time0,istep,&iinc,ibody,clearini,&mortar,springarea,
663  pslavsurf,pmastsurf,&reltime,&nasym));
664 
665  /* zc = damping matrix * eigenmodes */
666 
667  NNEW(zc,double,neq[1]*nev);
668  for(i=0;i<nev;i++){
669  FORTRAN(op,(&neq[1],&z[i*neq[1]],&zc[i*neq[1]],adc,auc,
670  jq,irow));
671  }
672 
673  /* cc is the reduced damping matrix (damping matrix mapped onto
674  space spanned by eigenmodes) */
675 
676  NNEW(cc,double,nev*nev);
677  for(i=0;i<nev;i++){
678  for(j=0;j<=i;j++){
679  for(k=0;k<neq[1];k++){
680  cc[i*nev+j]+=z[j*neq[1]+k]*zc[i*neq[1]+k];
681  }
682  }
683  }
684 
685  /* symmetric part of cc matrix */
686 
687  for(i=0;i<nev;i++){
688  for(j=i;j<nev;j++){
689  cc[i*nev+j]=cc[j*nev+i];
690  }
691  }
692  SFREE(zc);
693  }
694 
695  /* contact conditions */
696 
697  if(*nslavs==0){irenewxstate=1;}else{irenewxstate=0;}
698  inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
699  lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
700  &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
701  &mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,&ifacecount,
702  nmpc,&mpcfree,&memmpc_,
703  &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
704  iperturb,ikboun,nboun,co,istep,&xnoels);
705 
706  if(ncont!=0){
707 
708  if(mortar>0){
709  printf("*ERROR in dyna: modal dynamics cannot be combined with\n");
710  printf(" face-to-face penalty contact\n\n");
711  FORTRAN(stop,());
712  }
713 
714  if(dashpot){
715  printf("*ERROR in dyna: contact is not allowed in combination with dashpots\n");
716  FORTRAN(stop,());
717  }
718  RENEW(ipkon,ITG,*ne+*nslavs);
719  RENEW(lakon,char,8*(*ne+*nslavs));
720  if(*nener==1){
721  RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
722  }
723 
724  /* 11 instead of 10: last position is reserved for how
725  many dependent nodes are paired to this face */
726 
727  RENEW(kon,ITG,*nkon+11**nslavs);
728  if(*norien>0){
729  RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
730  for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
731  }
732  RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
733  for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
734  NNEW(cg,double,3*ncont);
735  NNEW(straight,double,16*ncont);
736 
737  /* internal state variables for contact */
738 
739  if((irenewxstate==1)&&(*nslavs!=0)&&(*nstate_>0)){
740  RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
741  for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
742  xstate[k]=0.;
743  }
744  }
745  if(*nstate_>0){
746  NNEW(xstateini,double,*nstate_*mi[0]*(*ne+*nslavs));
747  for(k=0;k<*nstate_*mi[0]*(*ne+*nslavs);++k){
748  xstateini[k]=xstate[k];
749  }
750  }
751 
752  NNEW(xmastnor,double,3*nmastnode[*ntie]);
753  NNEW(areaslav,double,ifacecount);
754  NNEW(springarea,double,2**nslavs);
755  NNEW(vini,double,mt**nk);
756  NNEW(bcontini,double,neq[1]);
757  NNEW(bcont,double,neq[1]);
758  NNEW(ikactcont,ITG,nactcont_);
759  }
760 
761  /* storing the element and topology information before introducing
762  contact elements */
763 
764  ne0=*ne;nkon0=*nkon;
765 
766  NNEW(zeta,double,nev);
767  NNEW(cstr,double,6);
768 
769  /* calculating the damping coefficients*/
770 
771  if(xmodal[10]<0){
772  for(i=0;i<nev;i++){
773  if(fabs(d[i])>(1.e-10)){
774  zeta[i]=(alpham+betam*d[i]*d[i])/(2.*d[i]);
775  }
776  else {
777  printf("*WARNING in dyna: one of the frequencies is zero\n");
778  printf(" no Rayleigh mass damping allowed\n");
779  zeta[i]=0.;
780  }
781  }
782  }
783  else{
784 
785  /*copy the damping coefficients for every eigenfrequency from xmodal[11....] */
786 
787  if(nev<(ITG)xmodal[10]){
788  imax=nev;
789  printf("*WARNING in dyna: too many modal damping coefficients applied\n");
790  printf(" damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
791  }
792  else{
793  imax=(ITG)xmodal[10];
794  }
795  for(i=0; i<imax; i++){
796  zeta[i]=xmodal[11+i];
797  }
798 
799  }
800 
801  /* modal decomposition of the initial conditions */
802  /* for cyclic symmetric structures the initial conditions
803  are assumed to be zero */
804 
805  NNEW(cd,double,nev);
806  NNEW(cv,double,nev);
807 
808  if(!cyclicsymmetry){
809  NNEW(temp_array1,double,neq[1]);
810  NNEW(temp_array2,double,neq[1]);
811  for(i=0;i<neq[1];i++){temp_array1[i]=0;temp_array2[i]=0;}
812 
813  /* displacement initial conditions */
814 
815  for(i=0;i<*nk;i++){
816  for(j=0;j<mt;j++){
817  if(nactdof[mt*i+j]!=0){
818  idof=nactdof[mt*i+j]-1;
819  temp_array1[idof]=vold[mt*i+j];
820  }
821  }
822  }
823 
824  FORTRAN(op,(&neq[1],temp_array1,temp_array2,adb,aub,jq,irow));
825 
826  for(i=0;i<neq[1];i++){
827  for(k=0;k<nev;k++){
828  cd[k]+=z[k*neq[1]+i]*temp_array2[i];
829  }
830  }
831 
832  /* velocity initial conditions */
833 
834  for(i=0;i<neq[1];i++){temp_array1[i]=0;temp_array2[i]=0;}
835  for(i=0;i<*nk;i++){
836  for(j=0;j<mt;j++){
837  if(nactdof[mt*i+j]!=0){
838  idof=nactdof[mt*i+j]-1;
839  temp_array1[idof]=veold[mt*i+j];
840  }
841  }
842  }
843 
844  FORTRAN(op,(&neq[1],temp_array1,temp_array2,adb,aub,jq,irow));
845 
846  for(i=0;i<neq[1];i++){
847  for(k=0;k<nev;k++){
848  cv[k]+=z[k*neq[1]+i]*temp_array2[i];
849  }
850  }
851 
852  SFREE(temp_array1);SFREE(temp_array2);
853 
854  }
855  NNEW(xforcact,double,*nforc);
856  NNEW(xforcdiff,double,*nforc);
857  NNEW(xloadact,double,2**nload);
858  NNEW(xloaddiff,double,2**nload);
859  NNEW(xbodyact,double,7**nbody);
860  NNEW(xbodydiff,double,7**nbody);
861 
862  /* copying the rotation axis and/or acceleration vector */
863 
864  for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
865  for(k=0;k<7**nbody;k++){xbodydiff[k]=xbody[k];}
866  NNEW(xbounact,double,*nboun);
867  NNEW(xboundiff,double,*nboun);
868  if(*ithermal==1) {NNEW(t1act,double,*nk);
869  NNEW(t1diff,double,*nk);}
870 
871  /* assigning the body forces to the elements */
872 
873  if(*nbody>0){
874  ifreebody=*ne+1;
875  NNEW(ipobody,ITG,2**ne);
876  for(k=1;k<=*nbody;k++){
877  FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
878  iendset,ialset,&inewton,nset,&ifreebody,&k));
879  RENEW(ipobody,ITG,2*(*ne+ifreebody));
880  }
881  RENEW(ipobody,ITG,2*(ifreebody-1));
882  }
883 
884  NNEW(b,double,neq[1]); /* load rhs vector and displacement solution vector */
885  NNEW(bp,double,neq[1]); /* velocity solution vector */
886  NNEW(bj,double,nev); /* response modal decomposition */
887  NNEW(bjp,double,nev); /* derivative of the response modal decomposition */
888  NNEW(ampli,double,*nam); /* instantaneous amplitude */
889 
890  /* constant coefficient of the linear amplitude function */
891 
892  NNEW(aa,double,nev);
893  NNEW(aanew,double,nev);
894  NNEW(aamech,double,nev);
895 
896  /* linear coefficient of the linear amplitude function */
897 
898  NNEW(bb,double,nev);
899 
900  NNEW(v,double,mt**nk);
901  NNEW(fn,double,mt**nk);
902  NNEW(stn,double,6**nk);
903  NNEW(inum,ITG,*nk);
904  strcpy1(&cflag[0],&filab[4],1);
905  FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
906  nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat));
907 
908  if(*ithermal>1) {NNEW(qfn,double,3**nk);
909  NNEW(qfx,double,3*mi[0]**ne);}
910 
911  if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
912  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
913  if(strcmp1(&filab[609],"SDV ")==0) NNEW(xstaten,double,*nstate_**nk);
914  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
915 
916  NNEW(eei,double,6*mi[0]**ne);
917  if(*nener==1){
918  NNEW(stiini,double,6*mi[0]**ne);
919  NNEW(enerini,double,mi[0]**ne);}
920 
921  /* check for nonzero SPC's */
922 
923  iprescribedboundary=0;
924  for(i=0;i<*nboun;i++){
925  if(fabs(xboun[i])>1.e-10){
926  iprescribedboundary=1;
927  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
928  break;
929  }
930  }
931 
932 /* calculating the instantaneous loads (forces, surface loading,
933  centrifugal and gravity loading or temperature) at time 0
934  setting iabsload to 2 if user subroutine dload is used */
935 
936 /* FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,
937  xload,xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,
938  xbodyact,t1old,t1,t1act,iamt1,nk,
939  amta,namta,nam,ampli,&time0,&reltime,ttime,&dtime,ithermal,nmethod,
940  xbounold,xboun,xbounact,iamboun,nboun,
941  nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
942  co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi));*/
943 
944  FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
945  xloadold,xload,xloadact,iamload,nload,ibody,xbody,
946  nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
947  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
948  nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
949  ndirboun,nodeforc,
950  ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
951  nelemload,sideload,mi,
952  xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
953  &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
954  fn));
955 
956  if(iabsload==2) NNEW(bold,double,neq[1]);
957 
958  /* calculating the instantaneous loading vector at time 0 */
959 
960  NNEW(ikactmech,ITG,neq[1]);
961  nactmech=0;
962  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
963  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
964  nforc,nelemload,sideload,xloadact,nload,xbodyact,ipobody,nbody,
965  cgr,b,nactdof,&neq[1],nmethod,
966  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,alcon,
967  nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,t0,t1act,
968  ithermal,iprestr,vold,iperturb,iexpl,plicon,
969  nplicon,plkcon,nplkcon,npmat_,ttime,&time0,istep,&iinc,&dtime,
970  physcon,ibody,xbodyold,&reltime,veold,matname,mi,ikactmech,
971  &nactmech,ielprop,prop));
972 
973  /* correction for nonzero SPC's */
974 
975  if(iprescribedboundary){
976 
977  if(cyclicsymmetry){
978  printf("*ERROR in dyna: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
979  FORTRAN(stop,());
980  }
981 
982  if(*idrct!=1){
983  printf("*ERROR in dyna: variable increment length is not allwed in combination with prescribed boundaries\n");
984  FORTRAN(stop,());
985  }
986 
987  /* LU decomposition of the stiffness matrix */
988 
989  if(*isolver==0){
990 #ifdef SPOOLES
991  spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
992  &symmetryflag,&inputformat,&nzs[2]);
993 #else
994  printf("*ERROR in dyna: the SPOOLES library is not linked\n\n");
995  FORTRAN(stop,());
996 #endif
997  }
998  else if(*isolver==4){
999 #ifdef SGI
1000  token=1;
1001  sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
1002 #else
1003  printf("*ERROR in dyna: the SGI library is not linked\n\n");
1004  FORTRAN(stop,());
1005 #endif
1006  }
1007  else if(*isolver==5){
1008 #ifdef TAUCS
1009  tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
1010 #else
1011  printf("*ERROR in dyna: the TAUCS library is not linked\n\n");
1012  FORTRAN(stop,());
1013 #endif
1014  }
1015  else if(*isolver==7){
1016 #ifdef PARDISO
1017  pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
1018  &symmetryflag,&inputformat,jq,&nzs[2]);
1019 #else
1020  printf("*ERROR in dyna: the PARDISO library is not linked\n\n");
1021  FORTRAN(stop,());
1022 #endif
1023  }
1024 
1025  NNEW(bact,double,neq[1]);
1026  NNEW(bmin,double,neq[1]);
1027  NNEW(bv,double,neq[1]);
1028  NNEW(bprev,double,neq[1]);
1029  NNEW(bdiff,double,neq[1]);
1030 
1031  init=1;
1032  dynboun(amta,namta,nam,ampli,&time0,ttime,&dtime,xbounold,xboun,
1033  xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1034  aub,icol,irow,neq,nzs,&sigma,b,isolver,
1035  &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1036  bprev,bdiff,&nactmech,&iabsload,&iprev);
1037  init=0;
1038  }
1039 
1040 /* creating contact elements and calculating the contact forces
1041  (normal and shear) */
1042 
1043  if(ncont!=0){
1044  DMEMSET(bcont,0,neq[1],0.);
1045  contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
1046  ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
1047  co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
1048  &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
1049  &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,iperturb,
1050  ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
1051  itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
1052  imastnode,nmastnode,xmastnor,filab,mcs,ics,
1053  &nasym,xnoels,&mortar,pslavsurf,pmastsurf,clearini,&theta);
1054 
1055  RENEW(ikactcont,ITG,nactcont_);
1056  DMEMSET(ikactcont,0,nactcont_,0.);
1057  nactcont=0;
1058 
1059  for(i=ne0;i<*ne;i++){
1060  indexe=ipkon[i];
1061  imat=ielmat[mi[2]*i];
1062  kodem=nelcon[2*imat-2];
1063  for(j=0;j<8;j++){lakonl[j]=lakon[8*i+j];}
1064  nope=atoi(&lakonl[7])+1;
1065  for(j=0;j<nope;j++){
1066  konl[j]=kon[indexe+j];
1067  for(j1=0;j1<3;j1++){
1068  xl[j*3+j1]=co[3*(konl[j]-1)+j1];
1069  voldl[mt*j+j1+1]=vold[mt*(konl[j]-1)+j1+1];
1070  veoldl[mt*j+j1+1]=veold[mt*(konl[j]-1)+j1+1];
1071  }
1072  }
1073  konl[nope]=kon[indexe+nope];
1074 
1075  FORTRAN(springforc_n2f,(xl,konl,voldl,&imat,elcon,nelcon,elas,
1076  fnl,ncmat_,ntmat_,&nope,lakonl,&t1l,&kodem,elconloc,
1077  plicon,nplicon,npmat_,&senergy,&iener,cstr,mi,
1078  &springarea[2*(konl[nope]-1)],nmethod,&ne0,nstate_,
1079  xstateini,xstate,&reltime,&ielas));
1080 
1081  storecontactdof(&nope,nactdof,&mt,konl,&ikactcont,&nactcont,
1082  &nactcont_,bcont,fnl,ikmpc,nmpc,ilmpc,ipompc,nodempc,
1083  coefmpc);
1084 
1085  }
1086  if(nactcont>100){nactcont_=nactcont;}else{nactcont_=100;}
1087  RENEW(ikactcont,ITG,nactcont_);
1088 
1089  }
1090 
1091  iit=1;
1092 
1093  /* load at the start of a new step:
1094  mechanical loading without contact */
1095 
1096  if(!cyclicsymmetry){
1097  for(i=0;i<nev;i++){
1098  i2=(long long)i*neq[1];
1099  aamech[i]=0.;
1100  if(nactmech<neq[1]/2){
1101  for(j=0;j<nactmech;j++){
1102  aamech[i]+=z[i2+ikactmech[j]]*b[ikactmech[j]];
1103  }
1104  }else{
1105  for(j=0;j<neq[1];j++){
1106  aamech[i]+=z[i2+j]*b[j];
1107  }
1108  }
1109  aanew[i]=aamech[i];
1110  if(ncont!=0){
1111  for(j=0;j<nactcont;j++){
1112  aanew[i]+=z[i2+ikactcont[j]]*bcont[ikactcont[j]];
1113  }
1114  }
1115  }
1116  }else{
1117  for(i=0;i<nev;i++){aamech[i]=0.;}
1118  for(j=0;j<nactmech;j++){
1119  FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
1120  if(id!=0){
1121  if(izdof[id-1]==ikactmech[j]){
1122  for(i=0;i<nev;i++){
1123  aamech[i]+=z[(long long)i*nzdof+id-1]*b[ikactmech[j]];
1124  }
1125  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1126  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1127  }
1128  memcpy(&aanew[0],&aamech[0],sizeof(double)*nev);
1129  if(ncont!=0){
1130  for(j=0;j<nactcont;j++){
1131  FORTRAN(nident,(izdof,&ikactcont[j],&nzdof,&id));
1132  if(id!=0){
1133  if(izdof[id-1]==ikactcont[j]){
1134  for(i=0;i<nev;i++){
1135  aanew[i]+=z[(long long)i*nzdof+id-1]*bcont[ikactcont[j]];
1136  }
1137  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1138  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1139  }
1140  }
1141  }
1142 
1143  /* check whether integration point values are requested; if not,
1144  the stress fields do not have to be allocated */
1145 
1146  intpointvar=0;
1147  if(*ithermal<=1){
1148 
1149  /* mechanical */
1150 
1151  if((strcmp1(&filab[174],"S")==0)||
1152  (strcmp1(&filab[261],"E")==0)||
1153  (strcmp1(&filab[348],"RF")==0)||
1154  (strcmp1(&filab[435],"PEEQ")==0)||
1155  (strcmp1(&filab[522],"ENER")==0)||
1156  (strcmp1(&filab[609],"SDV")==0)||
1157  (strcmp1(&filab[1044],"ZZS")==0)||
1158  (strcmp1(&filab[1044],"ERR")==0)||
1159  (strcmp1(&filab[1479],"PHS")==0)||
1160  (strcmp1(&filab[1653],"MAXS")==0)||
1161  (strcmp1(&filab[2175],"CONT")==0)||
1162  (strcmp1(&filab[2262],"CELS")==0)) intpointvar=1;
1163  for(i=0;i<*nprint;i++){
1164  if((strcmp1(&prlab[6*i],"S")==0)||
1165  (strcmp1(&prlab[6*i],"E")==0)||
1166  (strcmp1(&prlab[6*i],"PEEQ")==0)||
1167  (strcmp1(&prlab[6*i],"ENER")==0)||
1168  (strcmp1(&prlab[6*i],"SDV")==0)||
1169  (strcmp1(&prlab[6*i],"CDIS")==0)||
1170  (strcmp1(&prlab[6*i],"CSTR")==0)||
1171  (strcmp1(&prlab[6*i],"CELS")==0)||
1172  (strcmp1(&prlab[6*i],"RF")==0)) {intpointvar=1;break;}
1173  }
1174  }else{
1175 
1176  /* thermal */
1177 
1178  if((strcmp1(&filab[696],"HFL")==0)||
1179  (strcmp1(&filab[783],"RFL")==0)) intpointvar=1;
1180  for(i=0;i<*nprint;i++){
1181  if((strcmp1(&prlab[6*i],"HFL")==0)||
1182  (strcmp1(&prlab[6*i],"RFL")==0)) {intpointvar=1;break;}
1183  }
1184  }
1185 
1186  if((intpointvar==1)) NNEW(stx,double,6*mi[0]**ne);
1187 
1188  /* major loop */
1189 
1190  resultmaxprev=0.;
1191  resultmax=0.;
1192 
1193  while(1.-theta>1.e-6){
1194 
1195  time0=time;
1196 
1197 // printf("\nnew increment\n");
1198 
1199  if(*nener==1){
1200  memcpy(&enerini[0],&ener[0],sizeof(double)*mi[0]*ne0);
1201  if(*ithermal!=2){
1202  memcpy(&stiini[0],&sti[0],sizeof(double)*6*mi[0]*ne0);
1203  }
1204  }
1205 
1206  if(ncont!=0){
1207  if(nmdnode!=0){
1208  for(i=0;i<nmdnode;i++){
1209  i1=mt*(imdnode[i]-1);
1210  for(j=kmin;j<=kmax;j++){
1211  vini[i1+j]=vold[i1+j];
1212  }
1213  }
1214  }else{
1215  memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
1216  }
1217  if(*nstate_>0){
1218  for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
1219  xstateini[k]=xstate[k];
1220  }
1221  }
1222  }
1223  iinc++;
1224  jprint++;
1225 
1226  if(dashpot)RENEW(rpar,double,4+nev*(3+nev));
1227 
1228  /* check for max. # of increments */
1229 
1230  if(iinc>*jmax){
1231  printf(" *ERROR in dyna: max. # of increments reached\n\n");
1232  FORTRAN(stop,());
1233  }
1234 
1235  if(iinc>1){
1236  memcpy(&cd[0],&bj[0],sizeof(double)*nev);
1237  memcpy(&cv[0],&bjp[0],sizeof(double)*nev);
1238  }
1239 
1240 
1241  if((*idrct!=1)&&(iinc!=1)){
1242 
1243  /* increasing the increment size */
1244 
1245  dthetaold=dtheta;
1246  dtheta=dthetaref*dd;
1247 
1248  /* check increment length whether
1249  - it does not exceed tmax
1250  - the step length is not exceeded
1251  - a time point is not exceeded */
1252 
1253  dthetaref=dtheta;
1254  checkinclength(&time0,ttime,&theta,&dtheta,idrct,tper,tmax,
1255  tmin,ctrl, amta,namta,itpamp,&inext,&dthetaref,&itp,
1256  &jprint,jout);
1257  }
1258 
1259  reltime=theta+dtheta;
1260  time=reltime**tper;
1261  dtime=dtheta**tper;
1262 
1263  /* calculating the instantaneous loads (forces, surface loading,
1264  centrifugal and gravity loading or temperature) */
1265 
1266  FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
1267  xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1268  nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
1269  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
1270  nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
1271  ndirboun,nodeforc,
1272  ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
1273  nelemload,sideload,mi,
1274  xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
1275  &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
1276  fn));
1277 
1278  /* calculating the instantaneous loading vector */
1279 
1280  if(iabsload!=2){
1281  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1282  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcdiff,
1283  nforc,nelemload,sideload,xloaddiff,nload,xbodydiff,
1284  ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1285  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1286  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1287  t0,t1diff,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1288  nplicon,plkcon,nplkcon,
1289  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1290  xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1291  ielprop,prop));
1292  }else{
1293  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1294  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1295  nforc,nelemload,sideload,xloadact,nload,xbodyact,
1296  ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1297  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1298  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1299  t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1300  nplicon,plkcon,nplkcon,
1301  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1302  xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1303  ielprop,prop));
1304  }
1305 
1306  /* correction for nonzero SPC's */
1307 
1308  if(iprescribedboundary){
1309  dynboun(amta,namta,nam,ampli,&time,ttime,&dtime,
1310  xbounold,xboun,
1311  xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1312  aub,icol,irow,neq,nzs,&sigma,b,isolver,
1313  &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1314  bprev,bdiff,&nactmech,&iabsload,&iprev);
1315  }
1316 
1317  if(*idrct==0){
1318  bbmax=0.;
1319  if(iabsload!=2){
1320  if(nactmech<neq[1]/2){
1321  for(i=0;i<nactmech;i++){
1322  if(fabs(b[ikactmech[i]])>bbmax) bbmax=fabs(b[ikactmech[i]]);
1323  }
1324  }else{
1325  for(i=0;i<neq[1];i++){
1326  if(fabs(b[i])>bbmax) bbmax=fabs(b[i]);
1327  }
1328  }
1329  }else{
1330 
1331  /* bbmax is to be calculated from the difference of b and bold */
1332 
1333  if(nactmech<neq[1]/2){
1334  for(i=0;i<nactmech;i++){
1335  if(fabs(b[ikactmech[i]]-bold[ikactmech[i]])>bbmax)
1336  bbmax=fabs(b[ikactmech[i]]-bold[ikactmech[i]]);
1337  }
1338  }else{
1339  for(i=0;i<neq[1];i++){
1340  if(fabs(b[i])>bbmax) bbmax=fabs(b[i]-bold[i]);
1341  }
1342  }
1343 
1344  /* copy b into bold */
1345 
1346  if(nactmech<neq[1]/2){
1347  for(i=0;i<nactmech;i++){
1348  bold[ikactmech[i]]=b[ikactmech[i]];
1349  }
1350  }else{
1351  memcpy(&bold[0],&b[0],sizeof(double)*neq[1]);
1352  }
1353  }
1354 
1355  /* check for size of mechanical force */
1356 
1357  if((bbmax>deltmx)&&(((itp==1)&&(dtheta>*tmin))||(itp==0))){
1358 
1359  /* force increase too big: increment size is decreased */
1360 
1361  if(iabsload==0) iabsload=1;
1362  dtheta=dtheta*deltmx/bbmax;
1363  dthetaref=dtheta;
1364  if(itp==1){
1365  inext--;
1366  itp=0;
1367  }
1368 
1369  /* check whether the new increment size is not too small */
1370 
1371  if(dtheta<*tmin){
1372  dtheta=*tmin;
1373  }
1374 
1375  reltime=theta+dtheta;
1376  time=reltime**tper;
1377  dtime=dtheta**tper;
1378 
1379  /* calculating the instantaneous loads (forces, surface loading,
1380  centrifugal and gravity loading or temperature) */
1381 
1382  FORTRAN(temploaddiff,(xforcold,xforc,xforcact,iamforc,nforc,
1383  xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1384  nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,
1385  namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,
1386  nmethod,xbounold,xboun,xbounact,iamboun,nboun,nodeboun,
1387  ndirboun,nodeforc,
1388  ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,ikboun,ilboun,
1389  nelemload,sideload,mi,
1390  xforcdiff,xloaddiff,xbodydiff,t1diff,xboundiff,&iabsload,
1391  &iprescribedboundary,ntrans,trab,inotr,veold,nactdof,bcont,
1392  fn));
1393 
1394  /* calculating the instantaneous loading vector */
1395 
1396  if(iabsload!=2){
1397  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1398  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcdiff,
1399  nforc,nelemload,sideload,xloaddiff,nload,xbodydiff,
1400  ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1401  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1402  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1403  t0,t1diff,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1404  nplicon,plkcon,nplkcon,
1405  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1406  xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1407  ielprop,prop));
1408  }else{
1409  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1410  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcact,
1411  nforc,nelemload,sideload,xloadact,nload,xbodyact,
1412  ipobody,nbody,cgr,b,nactdof,&neq[1],nmethod,
1413  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1414  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1415  t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1416  nplicon,plkcon,nplkcon,
1417  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1418  xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
1419  ielprop,prop));
1420  }
1421 
1422  /* correction for nonzero SPC's */
1423 
1424  if(iprescribedboundary){
1425  dynboun(amta,namta,nam,ampli,&time,ttime,&dtime,
1426  xbounold,xboun,
1427  xbounact,iamboun,nboun,nodeboun,ndirboun,ad,au,adb,
1428  aub,icol,irow,neq,nzs,&sigma,b,isolver,
1429  &alpham,&betam,nzl,&init,bact,bmin,jq,amname,bv,
1430  bprev,bdiff,&nactmech,&iabsload,&iprev);
1431  }
1432  if(iabsload==1) iabsload=0;
1433  }
1434 
1435  if(ncont!=0){
1436  for(i=0;i<nactcont;i++){
1437  jdof=ikactcont[i];
1438  bcontini[jdof]=bcont[jdof];
1439  }
1440  }
1441 
1442  }
1443 
1444  /* step length is OK for mechanical load
1445  calculating equation for linearized loading */
1446 
1447  /* load: actual mechanical load +
1448  contact from last increment */
1449 
1450  if(!cyclicsymmetry){
1451  for(i=0;i<nev;i++){
1452  i2=(long long)i*neq[1];
1453  aa[i]=aanew[i];
1454  if(iabsload==2){aamech[i]=0.;}
1455  if(nactmech<neq[1]/2){
1456  for(j=0;j<nactmech;j++){
1457  aamech[i]+=z[i2+ikactmech[j]]*b[ikactmech[j]];
1458  }
1459  }else{
1460  for(j=0;j<neq[1];j++){
1461  aamech[i]+=z[i2+j]*b[j];
1462  }
1463  }
1464 
1465  aanew[i]=aamech[i];
1466  if(ncont!=0){
1467  for(j=0;j<nactcont;j++){
1468  aanew[i]+=z[i2+ikactcont[j]]*bcont[ikactcont[j]];
1469  }
1470  }
1471 
1472  bb[i]=(aanew[i]-aa[i])/dtime;
1473  aa[i]=aanew[i]-bb[i]*time;
1474  }
1475  }else{
1476  for(i=0;i<nev;i++){
1477  memcpy(&aa[0],&aanew[0],sizeof(double)*nev);
1478  if(iabsload==2){aamech[i]=0.;}
1479  }
1480  for(j=0;j<nactmech;j++){
1481  FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
1482  if(id!=0){
1483  if(izdof[id-1]==ikactmech[j]){
1484  for(i=0;i<nev;i++){
1485  aamech[i]+=z[(long long)i*nzdof+id-1]*b[ikactmech[j]];
1486  }
1487  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1488  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1489  }
1490  memcpy(&aanew[0],&aamech[0],sizeof(double)*nev);
1491  if(ncont!=0){
1492  for(j=0;j<nactcont;j++){
1493  FORTRAN(nident,(izdof,&ikactcont[j],&nzdof,&id));
1494  if(id!=0){
1495  if(izdof[id-1]==ikactcont[j]){
1496  for(i=0;i<nev;i++){
1497  aanew[i]+=z[(long long)i*nzdof+id-1]*bcont[ikactcont[j]];
1498  }
1499  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1500  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1501  }
1502  }
1503  for(i=0;i<nev;i++){
1504  bb[i]=(aanew[i]-aa[i])/(dtime);
1505  aa[i]=aanew[i]-bb[i]*time;
1506  }
1507  }
1508 
1509  /* calculating the response due to unchanged contact force during
1510  the increment */
1511 
1512  if(dashpot){
1513  FORTRAN(subspace,(d,aa,bb,cc,&alpham,&betam,&nev,xini,
1514  cd,cv,&time,rwork,&lrw,&iinc,jout,rpar,bj,
1515  iwork,&liw,&iddebdf,bjp));
1516  if(iddebdf==2){
1517  liw=56+2*nev;
1518  RENEW(iwork,ITG,liw);
1519  for(i=0;i<liw;i++){iwork[i]=0;}
1520  lrw=250+20*nev+4*nev*nev;
1521  RENEW(rwork,double,lrw);
1522  for(i=0;i<lrw;i++){rwork[i]=0.;}
1523  iddebdf=1;
1524  FORTRAN(subspace,(d,aa,bb,cc,&alpham,&betam,&nev,xini,
1525  cd,cv,&time,rwork,&lrw,&iinc,jout,rpar,bj,
1526  iwork,&liw,&iddebdf,bjp));
1527  }
1528  }
1529  else{
1530  for(l=0;l<nev;l++){
1531  zetaj=zeta[l];
1532  dj=d[l];
1533 
1534  /* zero eigenfrequency: rigid body mode */
1535 
1536  if(fabs(d[l])<=1.e-10){
1537  aai=aa[l];
1538  bbi=bb[l];
1539  tstart=time0;
1540  tend=time;
1541  sum=tend*(aai*time+
1542  tend*((bbi*time-aai)/2.-bbi*tend/3.))-
1543  tstart*(aai*time+
1544  tstart*((bbi*time-aai)/2.-bbi*tstart/3.));
1545  sump=tend*(aai+bbi*tend/2.)-tstart*(aai+bbi*tstart/2.);
1546  bj[l]=sum+cd[l]+dtime*cv[l];
1547  bjp[l]=sump+cv[l];
1548  }
1549 
1550  /* subcritical damping */
1551 
1552  else if(zetaj<1.-1.e-6){
1553  ddj=dj*sqrt(1.-zetaj*zetaj);
1554  h1=zetaj*dj;
1555  h2=h1*h1+ddj*ddj;
1556  h3=h1*h1-ddj*ddj;
1557  h4=2.*h1*ddj/h2;
1558  h14=h1/ddj;
1559  tstart=0.;
1560  FORTRAN(fsub,(&time,&dtime,&aa[l],&bb[l],&ddj,
1561  &h1,&h2,&h3,&h4,&func,&funcp));
1562  sum=func;sump=funcp;
1563  FORTRAN(fsub,(&time,&tstart,&aa[l],&bb[l],&ddj,
1564  &h1,&h2,&h3,&h4,&func,&funcp));
1565  sum-=func;sump-=funcp;
1566  fexp=exp(-h1*dtime);
1567  fsin=sin(ddj*dtime);
1568  fcos=cos(ddj*dtime);
1569 
1570  bj[l]=sum/ddj+fexp*(fcos+zetaj/sqrt(1.-zetaj*zetaj)*fsin)*cd[l]+
1571  fexp*fsin*cv[l]/ddj;
1572  bjp[l]=sump/ddj+fexp*((-h1+ddj*h14)*fcos+(-ddj-h1*h14)*fsin)*cd[l]
1573  +fexp*(-h1*fsin+ddj*fcos)*cv[l]/ddj;
1574 
1575  }
1576 
1577  /* supercritical damping */
1578 
1579  else if(zetaj>1.+1.e-6){
1580  ddj=dj*sqrt(zetaj*zetaj-1.);
1581  h1=ddj-zetaj*dj;
1582  h2=ddj+zetaj*dj;
1583  h3=1./h1;
1584  h4=1./h2;
1585  h5=h3*h3;
1586  h6=h4*h4;
1587  tstart=0.;
1588  FORTRAN(fsuper,(&time,&dtime,&aa[l],&bb[l],
1589  &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
1590  sum=func;sump=funcp;
1591  FORTRAN(fsuper,(&time,&tstart,&aa[l],&bb[l],
1592  &h1,&h2,&h3,&h4,&h5,&h6,&func,&funcp));
1593  sum-=func;sump-=funcp;
1594 
1595  fexm=exp(h1*dtime);
1596  fexp=exp(-h2*dtime);
1597  h14=zetaj*dj/ddj;
1598  bj[l]=sum/(2.*ddj)+(fexm+fexp)*cd[l]/2.+zetaj*(fexm-fexp)/(2.*
1599  sqrt(zetaj*zetaj-1.))*cd[l]+(fexm-fexp)*cv[l]/(2.*ddj);
1600  bjp[l]=sump/(2.*ddj)+(h1*fexm-h2*fexp)*cd[l]/2.
1601  +(h14*cd[l]+cv[l]/ddj)*(h1*fexm+h2*fexp)/2.;
1602  }
1603 
1604  /* critical damping */
1605 
1606  else{
1607  h1=zetaj*dj;
1608  h2=1./h1;
1609  h3=h2*h2;
1610  h4=h2*h3;
1611  tstart=0.;
1612  FORTRAN(fcrit,(&time,&dtime,&aa[l],&bb[l],&zetaj,&dj,
1613  &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
1614  sum=func;sump=funcp;
1615  FORTRAN(fcrit,(&time,&tstart,&aa[l],&bb[l],&zetaj,&dj,
1616  &ddj,&h1,&h2,&h3,&h4,&func,&funcp));
1617  sum-=func;sump-=funcp;
1618  fexp=exp(-h1*dtime);
1619  bj[l]=sum+fexp*((1.+h1*dtime)*cd[l]+dtime*cv[l]);
1620  bjp[l]=sump+fexp*(-h1*h1*dtime*cd[l]+
1621  (1.-h1*dtime)*cv[l]);
1622  }
1623  }
1624  }
1625 
1626  /* composing the response */
1627 
1628  if(iprescribedboundary){
1629  if(nmdnode==0){
1630  memcpy(&b[0],&bmin[0],sizeof(double)*neq[1]);
1631  memcpy(&bp[0],&bv[0],sizeof(double)*neq[1]);
1632  }else{
1633  for(i=0;i<nmddof;i++){
1634  b[imddof[i]]=bmin[imddof[i]];
1635  bp[imddof[i]]=bv[imddof[i]];
1636  }
1637  }
1638  }
1639  else{
1640  if(nmdnode==0){
1641  DMEMSET(b,0,neq[1],0.);
1642  DMEMSET(bp,0,neq[1],0.);
1643  }else{
1644  for(i=0;i<nmddof;i++){
1645  b[imddof[i]]=0.;
1646  bp[imddof[i]]=0.;
1647  }
1648  }
1649  }
1650 
1651  if(!cyclicsymmetry){
1652  if(nmdnode==0){
1653  for(i=0;i<neq[1];i++){
1654  for(j=0;j<nev;j++){
1655  b[i]+=bj[j]*z[(long long)j*neq[1]+i];
1656  bp[i]+=bjp[j]*z[(long long)j*neq[1]+i];
1657  }
1658  }
1659  }else{
1660  for(i=0;i<nmddof;i++){
1661  for(j=0;j<nev;j++){
1662  b[imddof[i]]+=bj[j]*z[(long long)j*neq[1]+imddof[i]];
1663  bp[imddof[i]]+=bjp[j]*z[(long long)j*neq[1]+imddof[i]];
1664  }
1665  }
1666  }
1667  }else{
1668  for(i=0;i<nmddof;i++){
1669  FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1670  if(id!=0){
1671  if(izdof[id-1]==imddof[i]){
1672  for(j=0;j<nev;j++){
1673  b[imddof[i]]+=bj[j]*z[(long long)j*nzdof+id-1];
1674  bp[imddof[i]]+=bjp[j]*z[(long long)j*nzdof+id-1];
1675  }
1676  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1677  }else{printf("*ERROR in dyna\n");FORTRAN(stop,());}
1678  }
1679  }
1680 
1681  /* update nonlinear MPC-coefficients (e.g. for rigid
1682  body MPC's */
1683 
1684  if(inonlinmpc==1){
1685  FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1686  nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1687  &maxlenmpc,ikmpc,ilmpc,&icascade,
1688  kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1689  &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1690  }
1691 
1692  /* calculating displacements/temperatures */
1693 
1694  FORTRAN(dynresults,(nk,v,ithermal,nactdof,vold,nodeboun,
1695  ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1696  b,bp,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1697  imdmpc,&nmdmpc,nmethod,&time));
1698 
1699  /* creating contact elements and calculating the contact forces
1700  based on the displacements at the end of the present increment */
1701 
1702  if(ncont!=0){
1703  dynacont(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
1704  ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
1705  nforc,nelemload,sideload,xload,nload,nactdof,neq,nzl,icol,
1706  irow,
1707  nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
1708  nrhcon,cocon,ncocon,alcon,nalcon,alzero,ielmat,ielorien,
1709  norien,orab,ntmat_,t0,t1,ithermal,prestr,iprestr,
1710  vold,iperturb,sti,nzs,tinc,tper,xmodal,veold,amname,amta,
1711  namta,nam,iamforc,iamload,iamt1,jout,filab,eme,xforcold,
1712  xloadold,t1old,iamboun,xbounold,iexpl,plicon,nplicon,plkcon,
1713  nplkcon,xstate,npmat_,matname,mi,ncmat_,nstate_,ener,jobnamec,
1714  ttime,set,nset,istartset,iendset,ialset,nprint,prlab,
1715  prset,nener,trab,inotr,ntrans,fmpc,cbody,ibody,xbody,nbody,
1716  xbodyold,istep,isolver,jq,output,mcs,nkon,mpcend,ics,cs,ntie,
1717  tieset,idrct,jmax,tmin,tmax,ctrl,itpamp,tietol,&iit,
1718  &ncont,&ne0,&reltime,&dtime,bcontini,bj,aux,iaux,bcont,
1719  &nev,v,&nkon0,&deltmx,&dtheta,&theta,&iprescribedboundary,
1720  &mpcfree,&memmpc_,itietri,koncont,cg,straight,&iinc,
1721  vini,aa,bb,aanew,d,z,zeta,b,&time0,&time,ipobody,
1722  xforcact,xloadact,t1act,xbounact,xbodyact,cd,cv,ampli,
1723  &dthetaref,bjp,bp,cstr,imddof,&nmddof,
1724  &ikactcont,&nactcont,&nactcont_,aamech,bprev,&iprev,&inonlinmpc,
1725  &ikactmech,&nactmech,imdnode,&nmdnode,imdboun,&nmdboun,
1726  imdmpc,&nmdmpc,&itp,&inext,imastop,
1727  nslavnode,islavnode,islavsurf,itiefac,areaslav,iponoels,
1728  inoels,springarea,izdof,&nzdof,fn,imastnode,nmastnode,xmastnor,
1729  xstateini,nslavs,&cyclicsymmetry,xnoels,&ielas,ielprop,prop);
1730  }
1731 
1732  theta+=dtheta;
1733 // (*ttime)+=dtime;
1734 
1735  /* check whether a time point was reached */
1736 
1737  if((*itpamp>0)&&(*idrct==0)){
1738  if(itp==1){
1739  jprint=*jout;
1740  }else{
1741  jprint=*jout+1;
1742  }
1743  }
1744 
1745  /* check whether output is needed */
1746 
1747  if((*jout==jprint)||(1.-theta<=1.e-6)){
1748  iout=2;
1749  jprint=0;
1750  }else if((*nener==1)){
1751  iout=-2;
1752  }else{
1753  iout=0;
1754  }
1755 
1756  if((iout==2)||(iout==-2)){
1757 
1758  /* deactivating the elements for which the stresses are not
1759  needed */
1760 
1761  if(nmdnode>0){
1762  if((intpointvar==1)){
1763  for(k=0;k<ne0;k++){
1764  if(ipkon[k]<-1){
1765  printf("*ERROR in dyna: contact remeshing of quadratic elements is not allowed\n\n");
1766  FORTRAN(stop,());
1767  }else if(ipkon[k]!=-1){
1768  ipkon[k]=-ipkon[k]-2;
1769  }
1770  }
1771  for(k=0;k<nmdelem;k++){
1772  ielem=imdelem[k]-1;
1773  ipkon[ielem]=-2-ipkon[ielem];
1774  }
1775  }
1776  }
1777 
1778  results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1779  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1780  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1781  ithermal,prestr,iprestr,filab,eme,emn,een,
1782  iperturb,f,fn,nactdof,&iout,qa,
1783  vold,b,nodeboun,ndirboun,xbounact,nboun,
1784  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1785  veold,accold,&bet,&gam,&dtime,&time,ttime,
1786  plicon,nplicon,plkcon,nplkcon,
1787  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1788  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1789  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1790  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1791  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,ikmpc,
1792  ilmpc,istep,&iinc,springarea,&reltime,&ne0,xforc,nforc,
1793  thicke,shcon,nshcon,
1794  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1795  &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1796  islavsurf,ielprop,prop);
1797 
1798  /* restoring */
1799 
1800  if(nmdnode>0){
1801  if((intpointvar==1)){
1802  for(k=0;k<ne0;k++){
1803  if(ipkon[k]<-1){ipkon[k]=-2-ipkon[k];}
1804  }
1805  }
1806  }
1807 
1808  if((*ithermal!=2)&&(intpointvar==1)){
1809  for(k=0;k<6*mi[0]*ne0;++k){
1810  sti[k]=stx[k];
1811  }
1812  }
1813  }
1814  if(iout==2){
1815  (*kode)++;
1816  if(strcmp1(&filab[1044],"ZZS")==0){
1817  NNEW(neigh,ITG,40**ne);
1818  NNEW(ipneigh,ITG,*nk);
1819  }
1820 
1821  ptime=*ttime+time;
1822  frd(co,&nkg,kon,ipkon,lakon,&neg,v,stn,inum,nmethod,
1823  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1824  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1825  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1826  mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1827  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1828  thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1829 
1830  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1831  }
1832 
1833  if(isteadystate==1){
1834 
1835  /* calculate maximum displacement/temperature */
1836 
1837  resultmax=0.;
1838  if(*ithermal<2){
1839  for(i=1;i<mt**nk;i=i+mt){
1840  if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1841  for(i=2;i<mt**nk;i=i+mt){
1842  if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1843  for(i=3;i<mt**nk;i=i+mt){
1844  if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1845  }else if(*ithermal==2){
1846  for(i=0;i<mt**nk;i=i+mt){
1847  if(fabs(v[i])>resultmax) resultmax=fabs(v[i]);}
1848  }else{
1849  printf("*ERROR in dyna: coupled temperature-displacement calculations are not allowed\n");
1850  }
1851  if(fabs((resultmax-resultmaxprev)/resultmax)<precision){
1852  break;
1853  }else{resultmaxprev=resultmax;}
1854  }
1855 
1856  }
1857 
1858  if((intpointvar==1)) SFREE(stx);
1859 
1860  /* calculating the displacements and velocities in all nodes as
1861  initial condition for the next step; only needed if
1862  - nonzero initial conditions are allowed (-> no cyclic symmetry)
1863  - the output was restricted (-> nmdnode nonzero) */
1864 
1865  if((nmdnode!=0)&&(!cyclicsymmetry)){
1866 
1867  /* determining the solution in the independent nodes */
1868 
1869  DMEMSET(b,0,neq[1],0.);
1870  DMEMSET(bp,0,neq[1],0.);
1871 
1872  for(i=0;i<neq[1];i++){
1873  for(j=0;j<nev;j++){
1874  b[i]+=bj[j]*z[(long long)j*neq[1]+i];
1875  bp[i]+=bjp[j]*z[(long long)j*neq[1]+i];
1876  }
1877  }
1878 
1879  /* update nonlinear MPC-coefficients (e.g. for rigid
1880  body MPC's */
1881 
1882  if(inonlinmpc==1){
1883  FORTRAN(nonlinmpc,(co,vold,ipompc,nodempc,coefmpc,labmpc,
1884  nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,
1885  &maxlenmpc,ikmpc,ilmpc,&icascade,
1886  kon,ipkon,lakon,ne,&reltime,&newstep,xboun,fmpc,
1887  &iit,&idiscon,&ncont,trab,ntrans,ithermal,mi));
1888  }
1889 
1890  /* calculating displacements/temperatures */
1891 
1892  nmdnode=0;
1893  FORTRAN(dynresults,(nk,v,ithermal,nactdof,vold,nodeboun,
1894  ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1895  b,bp,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1896  imdmpc,&nmdmpc,nmethod,&time));
1897  }
1898 
1899  SFREE(eei);
1900  SFREE(vbounact);
1901  SFREE(abounact);
1902 
1903  if(*nener==1){SFREE(stiini);SFREE(enerini);}
1904 
1905  if(strcmp1(&filab[261],"E ")==0) SFREE(een);
1906  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1907  if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstaten);
1908  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
1909  if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
1910 
1911  /* updating the loading at the end of the step;
1912  important in case the amplitude at the end of the step
1913  is not equal to one */
1914 
1915  for(k=0;k<*nboun;++k){xboun[k]=xbounact[k];}
1916  for(k=0;k<*nforc;++k){xforc[k]=xforcact[k];}
1917  for(k=0;k<2**nload;++k){xload[k]=xloadact[k];}
1918  for(k=0;k<7**nbody;k=k+7){xbody[k]=xbodyact[k];}
1919  if(*ithermal==1){
1920  for(k=0;k<*nk;++k){t1[k]=t1act[k];}
1921  }
1922 
1923  SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(adb);SFREE(d);
1924  SFREE(aub);SFREE(z);SFREE(b);SFREE(zeta);SFREE(bj);SFREE(cd);SFREE(cv);
1925  SFREE(xforcact);SFREE(xloadact);SFREE(xbounact);SFREE(aa);SFREE(bb);SFREE(aanew);
1926  SFREE(ampli);SFREE(xbodyact);SFREE(bjp);SFREE(bp);SFREE(aamech);SFREE(ikactmech);
1927  SFREE(xforcdiff);SFREE(xloaddiff);SFREE(xboundiff),SFREE(xbodydiff);
1928 
1929  if(*ithermal==1) {SFREE(t1act);SFREE(t1diff);}
1930 
1931  if(iprescribedboundary){
1932  if(*isolver==0){
1933 #ifdef SPOOLES
1934  spooles_cleanup();
1935 #endif
1936  }
1937  else if(*isolver==4){
1938 #ifdef SGI
1939  sgi_cleanup(token);
1940 #endif
1941  }
1942  else if(*isolver==5){
1943 #ifdef TAUCS
1944  tau_cleanup();
1945 #endif
1946  }
1947  else if(*isolver==7){
1948 #ifdef PARDISO
1949  pardiso_cleanup(&neq[1],&symmetryflag);
1950 #endif
1951  }
1952  SFREE(bact);SFREE(bmin);SFREE(bv);SFREE(bprev);SFREE(bdiff);
1953  }
1954 
1955  /* deleting the contact information */
1956 
1957 // *ne=ne0; *nkon=nkon0;
1958  if(ncont!=0){
1959  *ne=ne0; *nkon=nkon0;
1960  if(*nener==1){
1961  RENEW(ener,double,mi[0]**ne*2);
1962  }
1963  RENEW(ipkon,ITG,*ne);
1964  RENEW(lakon,char,8**ne);
1965  RENEW(kon,ITG,*nkon);
1966  if(*norien>0){
1967  RENEW(ielorien,ITG,mi[2]**ne);
1968  }
1969  RENEW(ielmat,ITG,mi[2]**ne);
1970  SFREE(cg);SFREE(straight);
1971 
1972  SFREE(vini);SFREE(bcont);SFREE(bcontini);SFREE(ikactcont);
1973 
1974  SFREE(imastop);SFREE(itiefac);SFREE(islavsurf);SFREE(islavnode);
1975  SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
1976  SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
1977  SFREE(springarea);SFREE(xmastnor);
1978 
1979  SFREE(areaslav);
1980 
1981  if(*nstate_>0){SFREE(xstateini);}
1982 
1983  }
1984 
1985  if(!cyclicsymmetry){
1986  SFREE(ad);SFREE(au);
1987  }else{
1988  SFREE(adbe); SFREE(aube);SFREE(icole); SFREE(irowe); SFREE(jqe);SFREE(izdof);
1989  SFREE(nm);
1990 
1991  *nk/=nsectors;
1992  *ne/=nsectors;
1993  *nkon/=nsectors;
1994  *nboun/=nsectors;
1995  neq[1]=neq[1]*2/nsectors;
1996 
1997  RENEW(ialset,ITG,nalset_);
1998 
1999  /* restore the infomration in istartset and iendset */
2000 
2001  for(j=0; j<*nset; j++){
2002  istartset[j]=istartset_[j];
2003  iendset[j]=iendset_[j];
2004  }
2005  SFREE(istartset_);
2006  SFREE(iendset_);
2007 
2008  RENEW(co,double,3**nk);
2009  if((*ithermal!=0)&&(*nam>0)) RENEW(iamt1,ITG,*nk);
2010  RENEW(nactdof,ITG,mt**nk);
2011  if(*ntrans>0) RENEW(inotr,ITG,2**nk);
2012  RENEW(kon,ITG,*nkon);
2013  RENEW(ipkon,ITG,*ne);
2014  RENEW(lakon,char,8**ne);
2015  RENEW(ielmat,ITG,mi[2]**ne);
2016  if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne);
2017  RENEW(nodeboun,ITG,*nboun);
2018  RENEW(ndirboun,ITG,*nboun);
2019  if(*nam>0) RENEW(iamboun,ITG,*nboun);
2020  RENEW(xboun,double,*nboun);
2021  RENEW(xbounold,double,*nboun);
2022  RENEW(ikboun,ITG,*nboun);
2023  RENEW(ilboun,ITG,*nboun);
2024 
2025  /* recovering the original multiple point constraints */
2026 
2027  RENEW(ipompc,ITG,*nmpc);
2028  RENEW(nodempc,ITG,3**mpcend);
2029  RENEW(coefmpc,double,*mpcend);
2030  RENEW(labmpc,char,20**nmpc+1);
2031  RENEW(ikmpc,ITG,*nmpc);
2032  RENEW(ilmpc,ITG,*nmpc);
2033  RENEW(fmpc,double,*nmpc);
2034 
2035  *nmpc=nmpcold;
2036  *mpcend=mpcendold;
2037  for(i=0;i<*nmpc;i++){ipompc[i]=ipompcold[i];}
2038  for(i=0;i<3**mpcend;i++){nodempc[i]=nodempcold[i];}
2039  for(i=0;i<*mpcend;i++){coefmpc[i]=coefmpcold[i];}
2040  for(i=0;i<20**nmpc;i++){labmpc[i]=labmpcold[i];}
2041  for(i=0;i<*nmpc;i++){ikmpc[i]=ikmpcold[i];}
2042  for(i=0;i<*nmpc;i++){ilmpc[i]=ilmpcold[i];}
2043  SFREE(ipompcold);SFREE(nodempcold);SFREE(coefmpcold);
2044  SFREE(labmpcold);SFREE(ikmpcold);SFREE(ilmpcold);
2045 
2046  RENEW(vold,double,mt**nk);
2047  RENEW(veold,double,mt**nk);
2048  RENEW(eme,double,6*mi[0]**ne);
2049  RENEW(sti,double,6*mi[0]**ne);
2050  if(*nener==1)RENEW(ener,double,mi[0]**ne*2);
2051 
2052 /* distributed loads */
2053 
2054  for(i=0;i<*nload;i++){
2055  if(nelemload[2*i]<nsectors){
2056  nelemload[2*i]-=*ne*nelemload[2*i+1];
2057  }else{
2058  nelemload[2*i]-=*ne*(nelemload[2*i+1]-nsectors);
2059  }
2060  }
2061 
2062  /* sorting the elements with distributed loads */
2063 
2064  if(*nload>0){
2065  if(*nam>0){
2066  FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
2067  }else{
2068  FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
2069  }
2070  }
2071 
2072 /* point loads */
2073 
2074  for(i=0;i<*nforc;i++){
2075  if(nodeforc[2*i+1]<nsectors){
2076  nodeforc[2*i]-=*nk*nodeforc[2*i+1];
2077  }else{
2078  nodeforc[2*i]-=*nk*(nodeforc[2*i+1]-nsectors);
2079  }
2080  }
2081  }
2082 
2083  SFREE(xstiff);if(*nbody>0) SFREE(ipobody);
2084 
2085  if(dashpot){
2086  SFREE(xini);SFREE(rwork);SFREE(adc);SFREE(auc);SFREE(cc);
2087  SFREE(rpar);SFREE(iwork);}
2088 
2089  SFREE(cstr);
2090 
2091  SFREE(imddof);SFREE(imdnode);SFREE(imdboun);SFREE(imdmpc);SFREE(imdelem);
2092 
2093  if(iabsload==2) SFREE(bold);
2094 
2095  *ialsetp=ialset;
2096  *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
2097  *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
2098  *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;
2099  *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
2100  *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
2101  *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2102  *fmpcp=fmpc;*veoldp=veold;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
2103  *stip=sti;*xstatep=xstate;
2104 
2105  (*ttime)+=(*tper);
2106 
2107  return;
2108 }
void storecontactdof(ITG *nope, ITG *nactdof, ITG *mt, ITG *konl, ITG **ikactcontp, ITG *nactcont, ITG *nactcont_, double *bcont, double *fnl, ITG *ikmpc, ITG *nmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, double *coefmpc)
#define ITGFORMAT
Definition: CalculiX.h:52
subroutine subspace(d, aa, bb, cc, alpham, betam, nev, xini, cd, cv, time, rwork, lrw, m, jout, rpar, bj, iwork, liw, iddebdf, bjp)
Definition: subspace.f:19
subroutine init(nktet, inodfa, ipofa, netet_)
Definition: init.f:19
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
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))
ITG strcpy1(char *s1, const char *s2, ITG length)
Definition: strcpy1.c:24
subroutine op(n, x, y, ad, au, jq, irow)
Definition: op.f:25
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
subroutine fsub(time, t, a, b, dd, h1, h2, h3, h4, func, funcp)
Definition: fsub.f:19
void pardiso_cleanup(ITG *neq, ITG *symmetryflag)
void dynacont(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, ITG *nmpc, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, ITG *nactdof, ITG *neq, ITG *nzl, ITG *icol, ITG *irow, ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *cocon, ITG *ncocon, 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, double *vold, ITG *iperturb, double *sti, ITG *nzs, double *tinc, double *tper, double *xmodal, double *veold, char *amname, double *amta, ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload, ITG *iamt1, ITG *jout, char *filab, double *eme, double *xforcold, double *xloadold, double *t1old, ITG *iamboun, double *xbounold, ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *xstate, ITG *npmat_, char *matname, ITG *mi, ITG *ncmat_, ITG *nstate_, double *ener, char *jobnamec, double *ttime, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, ITG *nener, double *trab, ITG *inotr, ITG *ntrans, double *fmpc, char *cbody, ITG *ibody, double *xbody, ITG *nbody, double *xbodyold, ITG *istep, ITG *isolver, ITG *jq, char *output, ITG *mcs, ITG *nkon, ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset, ITG *idrct, ITG *jmax, double *tmin, double *tmax, double *ctrl, ITG *itpamp, double *tietol, ITG *iit, ITG *ncont, ITG *ne0, double *reltime, double *dtime, double *bcontini, double *bj, double *aux, ITG *iaux, double *bcont, ITG *nev, double *v, ITG *nkon0, double *deltmx, double *dtheta, double *theta, ITG *iprescribedboundary, ITG *mpcfree, ITG *memmpc_, ITG *itietri, ITG *koncont, double *cg, double *straight, ITG *iinc, double *vini, double *aa, double *bb, double *aanew, double *d, double *z, double *zeta, double *b, double *time0, double *time1, ITG *ipobody, double *xforcact, double *xloadact, double *t1act, double *xbounact, double *xbodyact, double *cd, double *cv, double *ampli, double *dthetaref, double *bjp, double *bp, double *cstr, ITG *imddof, ITG *nmddof, ITG **ikactcontp, ITG *nactcont, ITG *nactcont_, double *aamech, double *bprev, ITG *iprev, ITG *inonlinmpc, ITG **ikactmechp, ITG *nactmech, ITG *imdnode, ITG *nmdnode, ITG *imdboun, ITG *nmdboun, ITG *imdmpc, ITG *nmdmpc, ITG *itp, ITG *inext, ITG *imastop, ITG *nslavnode, ITG *islavnode, ITG *islavsurf, ITG *itiefac, double *areaslav, ITG *iponoels, ITG *inoels, double *springarea, ITG *izdof, ITG *nzdof, double *fn, ITG *imastnode, ITG *nmastnode, double *xmastnor, double *xstateini, ITG *nslavs, ITG *cyclicsymmetry, double *xnoels, ITG *ielas, ITG *ielprop, double *prop)
Definition: dynacont.c:38
#define DMEMSET(a, b, c, d)
Definition: CalculiX.h:45
subroutine fsuper(time, t, a, b, h1, h2, h3, h4, h5, h6, func, funcp)
Definition: fsuper.f:19
subroutine addimdnodedload(nelemload, sideload, ipkon, kon, lakon, iload, imdnode, nmdnode, ikmpc, ilmpc, ipompc, nodempc, nmpc, imddof, nmddof, nactdof, mi, imdmpc, nmdmpc, imdboun, nmdboun, ikboun, nboun, ilboun, ithermal)
Definition: addimdnodedload.f:19
subroutine stop()
Definition: stop.f:19
void sgi_cleanup(ITG token)
void tau_cleanup()
void expand(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, ITG *nmpc, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, ITG *nactdof, ITG *neq, ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun, 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, ITG *ithermal, double *prestr, ITG *iprestr, double *vold, ITG *iperturb, double *sti, ITG *nzs, double *adb, double *aub, char *filab, double *eme, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *xstate, ITG *npmat_, char *matname, ITG *mi, ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_, ITG *nstate_, ITG *mcs, ITG *nkon, double *ener, char *jobnamec, char *output, char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, ITG *nprint, char *prlab, char *prset, ITG *nener, double *trab, ITG *inotr, ITG *ntrans, double *ttime, double *fmpc, ITG *nev, double **z, ITG *iamboun, double *xbounold, ITG *nsectors, ITG *nm, ITG *icol, ITG *irow, ITG *nzl, ITG *nam, ITG *ipompcold, ITG *nodempcold, double *coefmpcold, char *labmpcold, ITG *nmpcold, double *xloadold, ITG *iamload, double *t1old, double *t1, ITG *iamt1, double *xstiff, ITG **icolep, ITG **jqep, ITG **irowep, ITG *isolver, ITG *nzse, double **adbep, double **aubep, ITG *iexpl, ITG *ibody, double *xbody, ITG *nbody, double *cocon, ITG *ncocon, char *tieset, ITG *ntie, ITG *imddof, ITG *nmddof, ITG *imdnode, ITG *nmdnode, ITG *imdboun, ITG *nmdboun, ITG *imdmpc, ITG *nmdmpc, ITG **izdofp, ITG *nzdof, ITG *nherm, double *xmr, double *xmi, char *typeboun, ITG *ielprop, double *prop)
Definition: expand.c:33
void sgi_factor(double *ad, double *au, double *adb, double *aub, double *sigma, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG token)
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 pardiso_factor(double *ad, double *au, double *adb, double *aub, double *sigma, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG *symmetryflag, ITG *inputformat, ITG *jq, ITG *nzs3)
subroutine createmdelem(imdnode, nmdnode, xforc, ikmpc, ilmpc, ipompc, nodempc, nmpc, imddof, nmddof, nactdof, mi, imdmpc, nmdmpc, imdboun, nmdboun, ikboun, nboun, ilboun, ithermal, imdelem, nmdelem, iponoel, inoel, prlab, prset, nprint, lakon, set, nset, ialset, ipkon, kon, istartset, iendset, nforc, ikforc, ilforc)
Definition: createmdelem.f:19
#define RENEW(a, b, c)
Definition: CalculiX.h:40
#define SFREE(a)
Definition: CalculiX.h:41
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 temploaddiff(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, xforcdiff, xloaddiff, xbodydiff, t1diff, xboundiff, iabsload, iprescribedboundary, ntrans, trab, inotr, veold, nactdof, bcont, fn)
Definition: temploaddiff.f:19
subroutine checktime(itpamp, namta, tinc, ttime, amta, tmin, inext, itp, istep)
Definition: checktime.f:19
void checkinclength(double *time, double *ttime, double *theta, double *dtheta, ITG *idrct, double *tper, double *tmax, double *tmin, double *ctrl, double *amta, ITG *namta, ITG *itpamp, ITG *inext, double *dthetaref, ITG *itp, ITG *jprint, ITG *jout)
subroutine nident(x, px, n, id)
Definition: nident.f:25
subroutine fcrit(time, t, a, b, ze, d, dd, h1, h2, h3, h4, func, funcp)
Definition: fcrit.f:19
subroutine isortiiddc(ix1, ix2, dy1, dy2, cy, n, kflag)
Definition: isortiiddc.f:5
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
subroutine isortiddc(ix, dy1, dy2, cy, n, kflag)
Definition: isortiddc.f:5
subroutine createinum(ipkon, inum, kon, lakon, nk, ne, cflag, nelemload, nload, nodeboun, nboun, ndirboun, ithermal, co, vold, mi, ielmat)
Definition: createinum.f:19
void spooles_factor(double *ad, double *au, double *adb, double *aub, double *sigma, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, ITG *symmetryflag, ITG *inputformat, ITG *nzs3)
subroutine addimdnodecload(nodeforc, iforc, imdnode, nmdnode, xforc, ikmpc, ilmpc, ipompc, nodempc, nmpc, imddof, nmddof, nactdof, mi, imdmpc, nmdmpc, imdboun, nmdboun, ikboun, nboun, ilboun, ithermal)
Definition: addimdnodecload.f:19
void spooles_cleanup()
#define ITG
Definition: CalculiX.h:51
subroutine elementpernode(iponoel, inoel, lakon, ipkon, kon, ne)
Definition: elementpernode.f:19
void tau_factor(double *ad, double **aup, double *adb, double *aub, double *sigma, ITG *icol, ITG **irowp, ITG *neq, ITG *nzs)
void dynboun(double *amta, ITG *namta, ITG *nam, double *ampli, double *time, double *ttime, double *dtime, double *xbounold, double *xboun, double *xbounact, ITG *iamboun, ITG *nboun, ITG *nodeboun, ITG *ndirboun, double *ad, double *au, double *adb, double *aub, ITG *icol, ITG *irow, ITG *neq, ITG *nzs, double *sigma, double *b, ITG *isolver, double *alpham, double *betam, ITG *nzl, ITG *init, double *bact, double *bmin, ITG *jq, char *amname, double *bv, double *bprev, double *bdiff, ITG *nactmech, ITG *icorrect, ITG *iprev)
Definition: dynboun.c:37
subroutine createmddof(imddof, nmddof, istartset, iendset, ialset, nactdof, ithermal, mi, imdnode, nmdnode, ikmpc, ilmpc, ipompc, nodempc, nmpc, imdmpc, nmdmpc, imdboun, nmdboun, ikboun, nboun, nset, ntie, tieset, set, lakon, kon, ipkon, labmpc, ilboun, filab, prlab, prset, nprint, ne, cyclicsymmetry)
Definition: createmddof.f:19
#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
subroutine springforc_n2f(xl, konl, vl, imat, elcon, nelcon, elas, fnl, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, senergy, iener, cstr, mi, springarea, nmethod, ne0, nstate_, xstateini, xstate, reltime, ielas)
Definition: springforc_n2f.f:19
subroutine dynresults(nk, v, ithermal, nactdof, vold, nodeboun, ndirboun, xboun, nboun, ipompc, nodempc, coefmpc, labmpc, nmpc, b, bp, veold, dtime, mi, imdnode, nmdnode, imdboun, nmdboun, imdmpc, nmdmpc, nmethod, time)
Definition: dynresults.f:19
Hosted by OpenAircraft.com, (Michigan UAV, LLC)