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

Go to the source code of this file.

Functions

void steadystate (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 *sti, 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 *xstate, 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 *ialset, 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 *ics, double *cs, ITG *mpcend, double *ctrl, ITG *ikforc, ITG *ilforc, double *thicke, ITG *nmat, char *typeboun, ITG *ielprop, double *prop)
 

Function Documentation

void steadystate ( 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 *  sti,
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 *  xstate,
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 ialset,
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 ics,
double *  cs,
ITG mpcend,
double *  ctrl,
ITG ikforc,
ITG ilforc,
double *  thicke,
ITG nmat,
char *  typeboun,
ITG ielprop,
double *  prop 
)

Definition at line 37 of file steadystate.c.

70  {
71 
72  char fneig[132]="",description[13]=" ",*lakon=NULL,*labmpc=NULL,
73  *labmpcold=NULL,cflag[1]=" ";
74 
75  ITG nev,i,j,k, *inum=NULL,*ipobody=NULL,inewton=0,nsectors,im,
76  iinc=0,l,iout,ielas,icmd,iprescribedboundary,ndata,nmd,nevd,
77  ndatatot,*iphaseforc=NULL,*iphaseload=NULL,*iphaseboun=NULL,
78  *isave=NULL,nfour,ii,ir,ic,mode,noddiam=-1,*nm=NULL,*islavact=NULL,
79  *kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,*inotr=NULL,
80  *nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,jj,
81  *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
82  *ilmpc=NULL,*ipompcold=NULL,*nodempcold=NULL,*ikmpcold=NULL,
83  *ilmpcold=NULL,nmpcold,mpcendold,kflag=2,*iamt1=NULL,ifreebody,
84  *itg=NULL,ntg=0,symmetryflag=0,inputformat=0,dashpot,nrhs=1,
85  *ipiv=NULL,info,nev2,ngraph=1,nkg,neg,iflag=1,idummy=1,imax,
86  nzse[3],mt=mi[1]+1,*ikactmech=NULL,nactmech,id,nasym=0,
87  *imddof=NULL,nmddof,*imdnode=NULL,nmdnode,*imdboun=NULL,nmdboun,
88  *imdmpc=NULL,nmdmpc,*izdof=NULL,nzdof,cyclicsymmetry,mortar=0,
89  *ikactmechr=NULL,*ikactmechi=NULL,nactmechr,nactmechi,intpointvar,
90  iforc,iload,ne0,*iponoel=NULL,*inoel=NULL,*imdelem=NULL,
91  nmdelem,*integerglob=NULL,*nshcon=NULL,nherm,icfd=0,*inomat=NULL,
92  *islavnode=NULL,*nslavnode=NULL,*islavsurf=NULL;
93 
94  long long i2;
95 
96  double *d=NULL, *z=NULL,*stiini=NULL,*vini=NULL,*freqnh=NULL,
97  *xforcact=NULL, *xloadact=NULL,y,*fr=NULL,*fi=NULL,*cc=NULL,
98  *t1act=NULL,*ampli=NULL, *aa=NULL,*bb=NULL,*vr=NULL,*vi=NULL,
99  *stn=NULL,*stx=NULL,*een=NULL,*adb=NULL,*xstiff=NULL,ptime,
100  *aub=NULL,*bjr=NULL, *bji=NULL,*xbodyr=NULL,*cdn=NULL,
101  *f=NULL, *fn=NULL, *xbounact=NULL,*epn=NULL,*xstateini=NULL,
102  *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
103  *qfx=NULL, *xbodyact=NULL, *cgr=NULL, *au=NULL,*xbodyi=NULL,
104  time,dtime,reltime,*co=NULL,*xboun=NULL,*xbounold=NULL,
105  physcon[1],qa[3],cam[5],accold[1],bet,gam,*emn=NULL,timem,
106  *ad=NULL,sigma=0.,alpham,betam,*fnr=NULL,*fni=NULL,*emeini=NULL,
107  fmin,fmax,bias,*freq=NULL,*xforcr=NULL,dd,pi,vreal,constant,
108  *xforci=NULL,*xloadr=NULL,*xloadi=NULL,*xbounr=NULL,*xbouni=NULL,
109  *br=NULL,*bi=NULL,*ubr=NULL,*ubi=NULL,*mubr=NULL,*mubi=NULL,
110  *wsave=NULL,*r=NULL,*xbounacttime=NULL,*btot=NULL,breal,tmin,tmax,
111  *vold=NULL,*eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,
112  *coefmpcold=NULL,*t0=NULL,*t1=NULL,*t1old=NULL,*adc=NULL,*auc=NULL,
113  *am=NULL,*bm=NULL,*zc=NULL,*e=NULL,*stnr=NULL,*stni=NULL,
114  *vmax=NULL,*stnmax=NULL,*va=NULL,*vp=NULL,*fric=NULL,*springarea=NULL,
115  *stna=NULL,*stnp=NULL,*bp=NULL,*eenmax=NULL,*clearini=NULL,
116  *doubleglob=NULL,*shcon=NULL,*veold=NULL,*xmr=NULL,*xmi=NULL,*eig=NULL,
117  *ax=NULL,*bx=NULL,*pslavsurf=NULL,*pmastsurf=NULL,xnull=0.,
118  *cdnr=NULL,*cdni=NULL,*tinc,*tper;
119 
120  /* dummy arguments for the call of expand*/
121 
122  char* tieset=NULL;
123  ITG *jqe=NULL,*icole=NULL,*irowe=NULL,ntie=0;
124  double *adbe=NULL,*aube=NULL;
125 
126  FILE *f1;
127 
128  ITG *ipneigh=NULL,*neigh=NULL;
129 
130 #ifdef SGI
131  ITG token;
132 #endif
133 
134  co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
135  ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
136  ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;veold=*veoldp;
137  xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
138  vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
139  coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
140  fmpc=*fmpcp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
141 
142  tinc=&timepar[0];
143  tper=&timepar[1];
144 
145  pi=4.*atan(1.);
146  iout=2;
147 
148  alpham=xmodal[0];
149  betam=xmodal[1];
150 
151  fmin=2.*pi*xmodal[2];
152  fmax=2.*pi*xmodal[3];
153  ndata=floor(xmodal[4]);
154  bias=xmodal[5];
155  nfour=floor(xmodal[6]);
156  if(nfour>0){
157  tmin=xmodal[7];
158  tmax=xmodal[8];
159  }
160 
161  /* determining nzl */
162 
163  *nzl=0;
164  for(i=neq[1];i>0;i--){
165  if(icol[i-1]>0){
166  *nzl=i;
167  break;
168  }
169  }
170 
171  /* opening the eigenvalue file and checking for cyclic symmetry */
172 
173  strcpy(fneig,jobnamec);
174  strcat(fneig,".eig");
175 
176  if((f1=fopen(fneig,"rb"))==NULL){
177  printf("*ERROR in steadystate: cannot open eigenvalue file for reading");
178  exit(0);
179  }
180 
181  printf(" *INFO in steadystate: if there are problems reading the .eig file this may be due to:\n");
182  printf(" 1) the nonexistence of the .eig file\n");
183  printf(" 2) other boundary conditions than in the input deck\n");
184  printf(" which created the .eig file\n\n");
185 
186  if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
187  printf("*ERROR in steadystate reading the cyclic symmetry flag in the eigenvalue file");
188  exit(0);
189  }
190 
191  if(fread(&nherm,sizeof(ITG),1,f1)!=1){
192  printf("*ERROR in steadystate reading the Hermitian flag in the eigenvalue file");
193  exit(0);
194  }
195 
196  /* creating imddof containing the degrees of freedom
197  retained by the user and imdnode containing the nodes */
198 
199  nmddof=0;nmdnode=0;nmdboun=0;nmdmpc=0;nmdelem=0;
200 
201  NNEW(imddof,ITG,*nk*3);
202  NNEW(imdnode,ITG,*nk);
203  NNEW(imdboun,ITG,*nboun);
204  NNEW(imdmpc,ITG,*nmpc);
205  FORTRAN(createmddof,(imddof,&nmddof,istartset,iendset,
206  ialset,nactdof,ithermal,mi,imdnode,&nmdnode,
207  ikmpc,ilmpc,ipompc,nodempc,nmpc,
208  imdmpc,&nmdmpc,imdboun,&nmdboun,ikboun,
209  nboun,nset,&ntie,tieset,set,lakon,kon,ipkon,labmpc,
210  ilboun,filab,prlab,prset,nprint,ne,&cyclicsymmetry));
211 
212  /* if results are requested in too many nodes, it is faster to
213  calculate the results in all nodes */
214 
215  if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
216  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
217  }
218 
219  if(nmdnode!=0){
220  if(!cyclicsymmetry){
221  for(i=0;i<*nload;i++){
222  iload=i+1;
223  FORTRAN(addimdnodedload,(nelemload,sideload,ipkon,kon,lakon,
224  &iload,imdnode,&nmdnode,ikmpc,ilmpc,ipompc,nodempc,nmpc,
225  imddof,&nmddof,nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
226  ikboun,nboun,ilboun,ithermal));
227  }
228 
229  for(i=0;i<*nforc;i++){
230  iforc=i+1;
231  FORTRAN(addimdnodecload,(nodeforc,&iforc,imdnode,&nmdnode,xforc,
232  ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
233  nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
234  ikboun,nboun,ilboun,ithermal));
235  }
236  }
237 
238  /* determining the elements belonging to a given node */
239 
240  NNEW(iponoel,ITG,*nk);
241  NNEW(inoel,ITG,2**nkon);
242  FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
243  NNEW(imdelem,ITG,*ne);
244 
245  /* storing the elements in which integration point results
246  are needed; storing the nodes which are needed to
247  calculate these results */
248 
249  FORTRAN(createmdelem,(imdnode,&nmdnode,xforc,
250  ikmpc,ilmpc,ipompc,nodempc,nmpc,imddof,&nmddof,
251  nactdof,mi,imdmpc,&nmdmpc,imdboun,&nmdboun,
252  ikboun,nboun,ilboun,ithermal,imdelem,&nmdelem,
253  iponoel,inoel,prlab,prset,nprint,lakon,set,nset,
254  ialset,ipkon,kon,istartset,iendset,nforc,
255  ikforc,ilforc));
256 
257  RENEW(imdelem,ITG,nmdelem);
258  SFREE(iponoel);SFREE(inoel);
259  }
260 
261  /* if results are requested in too many nodes, it is faster to
262  calculate the results in all nodes */
263 
264  if((nmdnode>*nk/2)&&(!cyclicsymmetry)){
265  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
266  }
267 
268  /* subtracting 1 to comply with the C-convention */
269 
270  for(i=0;i<nmddof;i++){imddof[i]-=1;}
271 
272  RENEW(imddof,ITG,nmddof);
273  RENEW(imdnode,ITG,nmdnode);
274  RENEW(imdboun,ITG,nmdboun);
275  RENEW(imdmpc,ITG,nmdmpc);
276 
277  nsectors=1;
278 
279  if(!cyclicsymmetry){
280 
281  nkg=*nk;
282  neg=*ne;
283 
284  if(fread(&nev,sizeof(ITG),1,f1)!=1){
285  printf("*ERROR in steadystate reading the number of eigenvalues in the eigenvalue file");
286  exit(0);
287  }
288 
289  if(nherm==1){
290  NNEW(d,double,nev);
291  if(fread(d,sizeof(double),nev,f1)!=nev){
292  printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
293  exit(0);
294  }
295 
296  for(i=0;i<nev;i++){
297  if(d[i]<0){d[i]=0.;}
298  }
299  }else{
300  NNEW(d,double,2*nev);
301  if(fread(d,sizeof(double),2*nev,f1)!=2*nev){
302  printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
303  exit(0);
304  }
305  }
306 
307  NNEW(ad,double,neq[1]);
308  NNEW(adb,double,neq[1]);
309  NNEW(au,double,nzs[2]);
310  NNEW(aub,double,nzs[1]);
311 
312  /* reading the stiffness matrix */
313 
314  if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
315  printf("*ERROR in steadystate reading the diagonal of the stiffness matrix in the eigenvalue file");
316  exit(0);
317  }
318 
319  if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
320  printf("*ERROR in steadystate reading the off-diagonals of the stiffness matrix in the eigenvalue file");
321  exit(0);
322  }
323 
324  /* reading the mass matrix */
325 
326  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
327  printf("*ERROR in steadystate reading the diagonal of the mass matrix in the eigenvalue file");
328  exit(0);
329  }
330 
331  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
332  printf("*ERROR in steadystate reading the off-diagonals of the mass matrix in the eigenvalue file");
333  exit(0);
334  }
335 
336  /* reading the eigenvectors */
337 
338  if(nherm==1){
339  NNEW(z,double,neq[1]*nev);
340  if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
341  printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
342  exit(0);
343  }
344  }else{
345  NNEW(z,double,2*neq[1]*nev);
346  if(fread(z,sizeof(double),2*neq[1]*nev,f1)!=2*neq[1]*nev){
347  printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
348  exit(0);
349  }
350  }
351 
352  /* reading the orthogonality matrices */
353 
354  if(nherm!=1){
355  NNEW(xmr,double,nev*nev);
356  NNEW(xmi,double,nev*nev);
357  if(fread(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
358  printf("*ERROR in steadystate reading the real orthogonality matrix to the eigenvalue file...");
359  exit(0);
360  }
361 
362  if(fread(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
363  printf("*ERROR in steadystate reading the imaginary orthogonality matrix to the eigenvalue file...");
364  exit(0);
365  }
366  }
367  }
368  else{
369  nev=0;
370  do{
371  if(fread(&nmd,sizeof(ITG),1,f1)!=1){
372  break;
373  }
374 
375  if(fread(&nevd,sizeof(ITG),1,f1)!=1){
376  printf("*ERROR in steadystate reading the number of eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
377  exit(0);
378  }
379 
380  /* reading the eigenvalues (complex for non-Hermitian systems) */
381 
382  if(nev==0){
383  if(nherm==1){NNEW(d,double,nevd);
384  }else{NNEW(d,double,2*nevd);}
385  NNEW(nm,ITG,nevd);
386  }else{
387  if(nherm!=1){
388  printf("*ERROR in steadystate: non-Hermitian systems cannot\n");
389  printf(" be combined with multiple modal diameters\n");
390  printf(" in cyclic symmetry calculations\n\n");
391  FORTRAN(stop,());
392  }
393  RENEW(d,double,nev+nevd);
394  RENEW(nm,ITG,nev+nevd);
395  }
396 
397  if(nherm==1){
398  if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
399  printf("*ERROR in steadystate reading the eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
400  exit(0);
401  }
402  for(i=nev;i<nev+nevd;i++){
403  if(d[i]<0){d[i]=0.;}
404  }
405  }else{
406  if(fread(&d[nev],sizeof(double),2*nevd,f1)!=2*nevd){
407  printf("*ERROR in steadystate reading the eigenvalues in the eigenvalue file...");
408  exit(0);
409  }
410  }
411 
412  for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
413 
414  if(nev==0){
415  NNEW(adb,double,neq[1]);
416  NNEW(aub,double,nzs[1]);
417 
418  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
419  printf("*ERROR in steadystate reading the diagonal of the mass matrix in the eigenvalue file");
420  exit(0);
421  }
422 
423  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
424  printf("*ERROR in steadystate reading the off-diagonals of the mass matrix in the eigenvalue file");
425  exit(0);
426  }
427  }
428 
429  /* reading the eigenvectors */
430 
431  if(nev==0){
432  NNEW(z,double,neq[1]*nevd);
433  }else{
434  RENEW(z,double,(long long)neq[1]*(nev+nevd));
435  }
436 
437  if(fread(&z[(long long)neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
438  printf("*ERROR in steadystate reading the eigenvectors for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
439  exit(0);
440  }
441 
442  /* reading the orthogonality matrices */
443 
444  if(nherm!=1){
445  NNEW(xmr,double,nev*nev);
446  NNEW(xmi,double,nev*nev);
447  if(fread(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
448  printf("*ERROR in steadystate reading the real orthogonality matrix to the eigenvalue file...");
449  exit(0);
450  }
451 
452  if(fread(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
453  printf("*ERROR in steadystate reading the imaginary orthogonality matrix to the eigenvalue file...");
454  exit(0);
455  }
456  }
457 
458  nev+=nevd;
459  }while(1);
460 
461  /* determining the maximum amount of sectors */
462 
463  for(i=0;i<*mcs;i++){
464  if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
465  }
466 
467  /* determining the maximum number of sectors to be plotted */
468 
469  for(j=0;j<*mcs;j++){
470  if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
471  }
472  nkg=*nk*ngraph;
473  neg=*ne*ngraph;
474 
475  /* allocating field for the expanded structure */
476 
477  RENEW(co,double,3**nk*nsectors);
478 
479  /* next line is necessary for multiple cyclic symmetry
480  conditions */
481 
482  for(i=3**nk;i<3**nk*nsectors;i++){co[i]=0.;}
483  if(*ithermal!=0){
484  RENEW(t0,double,*nk*nsectors);
485  RENEW(t1old,double,*nk*nsectors);
486  RENEW(t1,double,*nk*nsectors);
487  if(*nam>0) RENEW(iamt1,ITG,*nk*nsectors);
488  }
489  RENEW(nactdof,ITG,mt**nk*nsectors);
490  if(*ntrans>0) RENEW(inotr,ITG,2**nk*nsectors);
491  RENEW(kon,ITG,*nkon*nsectors);
492  RENEW(ipkon,ITG,*ne*nsectors);
493  for(i=*ne;i<*ne*nsectors;i++){ipkon[i]=-1;}
494  RENEW(lakon,char,8**ne*nsectors);
495  RENEW(ielmat,ITG,mi[2]**ne*nsectors);
496  if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne*nsectors);
497 // RENEW(z,double,(long long)neq[1]*nev*nsectors/2);
498 
499  RENEW(nodeboun,ITG,*nboun*nsectors);
500  RENEW(ndirboun,ITG,*nboun*nsectors);
501  if(*nam>0) RENEW(iamboun,ITG,*nboun*nsectors);
502  RENEW(xboun,double,*nboun*nsectors);
503  RENEW(xbounold,double,*nboun*nsectors);
504  RENEW(ikboun,ITG,*nboun*nsectors);
505  RENEW(ilboun,ITG,*nboun*nsectors);
506 
507  NNEW(ipompcold,ITG,*nmpc);
508  NNEW(nodempcold,ITG,3**mpcend);
509  NNEW(coefmpcold,double,*mpcend);
510  NNEW(labmpcold,char,20**nmpc);
511  NNEW(ikmpcold,ITG,*nmpc);
512  NNEW(ilmpcold,ITG,*nmpc);
513 
514  for(i=0;i<*nmpc;i++){ipompcold[i]=ipompc[i];}
515  for(i=0;i<3**mpcend;i++){nodempcold[i]=nodempc[i];}
516  for(i=0;i<*mpcend;i++){coefmpcold[i]=coefmpc[i];}
517  for(i=0;i<20**nmpc;i++){labmpcold[i]=labmpc[i];}
518  for(i=0;i<*nmpc;i++){ikmpcold[i]=ikmpc[i];}
519  for(i=0;i<*nmpc;i++){ilmpcold[i]=ilmpc[i];}
520  nmpcold=*nmpc;
521  mpcendold=*mpcend;
522 
523  RENEW(ipompc,ITG,*nmpc*nsectors);
524  RENEW(nodempc,ITG,3**mpcend*nsectors);
525  RENEW(coefmpc,double,*mpcend*nsectors);
526  RENEW(labmpc,char,20**nmpc*nsectors+1);
527  RENEW(ikmpc,ITG,*nmpc*nsectors);
528  RENEW(ilmpc,ITG,*nmpc*nsectors);
529  RENEW(fmpc,double,*nmpc*nsectors);
530 
531  /* reallocating the fields for the nodes in which the
532  solution has to be calculated */
533 
534  RENEW(imddof,ITG,neq[1]/2*nsectors);
535  RENEW(imdnode,ITG,*nk*nsectors);
536  RENEW(imdboun,ITG,*nboun*nsectors);
537  RENEW(imdmpc,ITG,*nmpc*nsectors);
538 
539 //izdofNNEW(// izdof,ITG,1);
540 
541  expand(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
542  ipompc,nodempc,coefmpc,labmpc,nmpc,nodeforc,ndirforc,xforc,
543  nforc,nelemload,sideload,xload,nload,nactdof,neq,
544  nmethod,ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,
545  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
546  t0,ithermal,prestr,iprestr,vold,iperturb,sti,nzs,
547  adb,aub,filab,eme,plicon,nplicon,plkcon,nplkcon,
548  xstate,npmat_,matname,mi,ics,cs,mpcend,ncmat_,
549  nstate_,mcs,nkon,ener,jobnamec,output,set,nset,istartset,
550  iendset,ialset,nprint,prlab,prset,nener,trab,
551  inotr,ntrans,ttime,fmpc,&nev,&z,iamboun,xbounold,
552  &nsectors,nm,icol,irow,nzl,nam,ipompcold,nodempcold,coefmpcold,
553  labmpcold,&nmpcold,xloadold,iamload,t1old,t1,iamt1,xstiff,&icole,&jqe,
554  &irowe,isolver,nzse,&adbe,&aube,iexpl,
555  ibody,xbody,nbody,cocon,ncocon,tieset,&ntie,imddof,&nmddof,
556  imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,&izdof,&nzdof,
557  &nherm,xmr,xmi,typeboun,ielprop,prop);
558 
559  RENEW(imddof,ITG,nmddof);
560  RENEW(imdnode,ITG,nmdnode);
561  RENEW(imdboun,ITG,nmdboun);
562  RENEW(imdmpc,ITG,nmdmpc);
563 
564  SFREE(vold);
565  NNEW(vold,double,mt**nk);
566  SFREE(veold);
567  NNEW(veold,double,mt**nk);
568  RENEW(eme,double,6*mi[0]**ne);
569 
570  if(*nener==1) RENEW(ener,double,mi[0]**ne);
571  }
572 
573  fclose(f1);
574 
575  /* allocating space for the friction coefficients */
576 
577  if(nherm==1){
578  NNEW(fric,double,nev);
579  }else{
580  NNEW(fric,double,2*nev);
581  }
582 
583  /* check whether there are dashpot elements */
584 
585  dashpot=0;
586  for(i=0;i<*ne;i++){
587  if(ipkon[i]==-1) continue;
588  if(strcmp1(&lakon[i*8],"ED")==0){
589  dashpot=1;break;}
590  }
591  if(dashpot){
592 
593  if(cyclicsymmetry){
594  printf("*ERROR in steadystate: dashpots are not allowed in combination with cyclic symmetry\n");
595  FORTRAN(stop,());
596  }
597 
598  if(nherm!=1){
599  printf("ERROR in steadystate: dashpots cannot be combined with non-Hermitian systems (in the present version of CalculiX)\n");
600  FORTRAN(stop,());
601  }
602 
603  /* cc is the reduced damping matrix (damping matrix mapped onto
604  space spanned by eigenmodes) */
605 
606  NNEW(cc,double,nev*nev);
607 /* nev2=2*nev;
608  NNEW(am,double,nev2*nev2);
609  NNEW(bm,double,nev2);
610  NNEW(ipiv,ITG,nev2);*/
611  }
612 
613  NNEW(inum,ITG,*nk);
614  strcpy1(&cflag[0],&filab[4],1);
615  FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
616  nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat));
617 
618  /* check whether integration point values are requested; if not,
619  the stress fields do not have to be allocated */
620 
621  intpointvar=0;
622  if(*ithermal<=1){
623 
624  /* mechanical */
625 
626  if((strcmp1(&filab[174],"S")==0)||
627  (strcmp1(&filab[261],"E")==0)||
628  (strcmp1(&filab[348],"RF")==0)||
629  (strcmp1(&filab[435],"PEEQ")==0)||
630  (strcmp1(&filab[522],"ENER")==0)||
631  (strcmp1(&filab[609],"SDV")==0)||
632  (strcmp1(&filab[1044],"ZZS")==0)||
633  (strcmp1(&filab[1044],"ERR")==0)||
634  (strcmp1(&filab[1479],"PHS")==0)||
635  (strcmp1(&filab[1653],"MAXS")==0)||
636  (strcmp1(&filab[2175],"CONT")==0)||
637  (strcmp1(&filab[2262],"CELS")==0)) intpointvar=1;
638  for(i=0;i<*nprint;i++){
639  if((strcmp1(&prlab[6*i],"S")==0)||
640  (strcmp1(&prlab[6*i],"E")==0)||
641  (strcmp1(&prlab[6*i],"PEEQ")==0)||
642  (strcmp1(&prlab[6*i],"ENER")==0)||
643  (strcmp1(&prlab[6*i],"SDV")==0)||
644  (strcmp1(&prlab[6*i],"RF")==0)) {intpointvar=1;break;}
645  }
646  }else{
647 
648  /* thermal */
649 
650  if((strcmp1(&filab[696],"HFL")==0)||
651  (strcmp1(&filab[783],"RFL")==0)) intpointvar=1;
652  for(i=0;i<*nprint;i++){
653  if((strcmp1(&prlab[6*i],"HFL")==0)||
654  (strcmp1(&prlab[6*i],"RFL")==0)) {intpointvar=1;break;}
655  }
656  }
657 
658  if(nfour<=0){
659 
660  /* harmonic excitation */
661 
662  NNEW(ikactmechr,ITG,neq[1]);
663  NNEW(ikactmechi,ITG,neq[1]);
664  nactmechr=0;nactmechi=0;
665 
666  /* result fields */
667 
668  if(intpointvar==1){
669  NNEW(fn,double,mt**nk);
670  NNEW(stnr,double,6**nk);
671  NNEW(stni,double,6**nk);
672  NNEW(stx,double,6*mi[0]**ne);
673  NNEW(eei,double,6*mi[0]**ne);
674 
675  if(*ithermal>1) {NNEW(qfn,double,3**nk);
676  NNEW(qfx,double,3*mi[0]**ne);}
677 
678  if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
679  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
680  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
681 
682  if(*nener==1){
683  NNEW(stiini,double,6*mi[0]**ne);
684  NNEW(enerini,double,mi[0]**ne);}
685  }
686 
687  /* determining the frequency data points */
688 
689  NNEW(freq,double,ndata*(nev+1));
690 
691  ndatatot=0.;
692  freq[0]=fmin;
693  if(fabs(fmax-fmin)<1.e-10){
694  ndatatot=1;
695  }else{
696 
697  /* copy the eigenvalues and sort them in ascending order
698  (important for values from distinct nodal diameters */
699 
700  NNEW(e,double,nev);
701  if(nherm==1){
702  for(i=0;i<nev;i++){e[i]=sqrt(d[i]);}
703  }else{
704 
705  /* for complex eigenvalues: sorting the real part */
706 
707  for(i=0;i<nev;i++){
708  e[i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/sqrt(2.);
709  }
710  }
711  FORTRAN(dsort,(e,&idummy,&nev,&iflag));
712 
713  for(i=0;i<nev;i++){
714  if(i!=0){
715  if(fabs(e[i]-e[i-1])<1.e-5){continue;}
716  }
717  if(e[i]>=fmin){
718  if(e[i]<=fmax){
719  for(j=1;j<ndata;j++){
720  y=-1.+2.*j/((double)(ndata-1));
721  if(fabs(y)<1.e-10){freq[ndatatot+j]=
722  (freq[ndatatot]+e[i])/2.;}
723  else{
724  freq[ndatatot+j]=(freq[ndatatot]+e[i])/2.+
725  (e[i]-freq[ndatatot])*pow(fabs(y),1./bias)*
726  y/(2.*fabs(y));
727  }
728  }
729  ndatatot+=ndata-1;
730  }
731  else{break;}
732  }
733  }
734 
735  SFREE(e);
736 
737  for(j=1;j<ndata;j++){
738  y=-1.+2.*j/((double)(ndata-1));
739  if(fabs(y)<1.e-10){freq[ndatatot+j]=(freq[ndatatot]+fmax)/2.;}
740  else{
741  freq[ndatatot+j]=(freq[ndatatot]+fmax)/2.+
742  (fmax-freq[ndatatot])*pow(fabs(y),1./bias)*
743  y/(2.*fabs(y));
744  }
745  }
746  ndatatot+=ndata;
747  }
748  RENEW(freq,double,ndatatot);
749 
750  /* check for nonzero SPC's */
751 
752  iprescribedboundary=0;
753  for(i=0;i<*nboun;i++){
754  if(fabs(xboun[i])>1.e-10){
755  iprescribedboundary=1;
756  nmdnode=0;nmddof=0;nmdboun=0;nmdmpc=0;
757  break;
758  }
759  }
760 
761  if((iprescribedboundary)&&(cyclicsymmetry)){
762  printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
763  FORTRAN(stop,());
764  }
765 
766  /* calculating the damping coefficients = friction coefficient*2*eigenvalue */
767 
768  if(xmodal[10]<0){
769  for(i=0;i<nev;i++){
770  if(nherm==1){
771  if(sqrt(d[i])>(1.e-10)){
772  fric[i]=(alpham+betam*d[i]);
773  }
774  else {
775  printf("*WARNING in steadystate: one of the frequencies is zero\n");
776  printf(" no Rayleigh mass damping allowed\n");
777  fric[i]=0.;
778  }
779  }else{
780  fric[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/
781  sqrt(2.);
782  fric[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])/
783  sqrt(2.);
784  if(d[2*i+1]<0.) fric[2*i+1]=-fric[2*i+1];
785  fric[2*i]=alpham+betam*fric[2*i];
786  fric[2*i+1]=betam*fric[2*i+1];
787  }
788  }
789  }
790  else{
791  if(iprescribedboundary){
792  printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with direct modal damping\n");
793  FORTRAN(stop,());
794  }
795 
796  /*copy the damping coefficients for every eigenfrequencie from xmodal[11....] */
797  if(nev<(ITG)xmodal[10]){
798  imax=nev;
799  printf("*WARNING in steadystate: too many modal damping coefficients applied\n");
800  printf(" damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
801  }
802  else{
803  imax=(ITG)xmodal[10];
804  }
805  for(i=0; i<imax; i++){
806  if(nherm==1){
807  fric[i]=2.*sqrt(d[i])*xmodal[11+i];
808  }else{
809  fric[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])/
810  sqrt(2.);
811  fric[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])/
812  sqrt(2.);
813  if(d[2*i+1]<0.) fric[2*i+1]=-fric[2*i+1];
814  fric[2*i]=2.*fric[2*i]*xmodal[11+i];
815  fric[2*i+1]=2.*fric[2*i+1]*xmodal[11+i];
816  }
817  }
818 
819  }
820 
821  /* check whether the loading is real or imaginary */
822 
823  NNEW(iphaseforc,ITG,*nforc);
824  for (i=0;i<*nforc;i++){
825  if(nodeforc[2*i+1]>=nsectors){
826  iphaseforc[i]=1;
827  }
828  }
829 
830  NNEW(iphaseload,ITG,*nload);
831  for (i=0;i<*nload;i++){
832  if(nelemload[2*i+1]>=nsectors){
833  iphaseload[i]=1;
834  }
835  }
836 
837  if(iprescribedboundary){
838  NNEW(iphaseboun,ITG,*nboun);
839  for (i=0;i<*nboun;i++){
840  if(nodeboun[i]>*nk){
841  iphaseboun[i]=1;
842  nodeboun[i]=nodeboun[i]-*nk;
843  }
844  }
845  }
846 
847  /* allocating actual loading fields */
848 
849  NNEW(xforcact,double,*nforc);
850  NNEW(xforcr,double,*nforc);
851  NNEW(xforci,double,*nforc);
852 
853  NNEW(xloadact,double,2**nload);
854  NNEW(xloadr,double,2**nload);
855  NNEW(xloadi,double,2**nload);
856 
857  NNEW(xbodyact,double,7**nbody);
858  NNEW(xbodyr,double,7**nbody);
859  NNEW(xbodyi,double,7**nbody);
860  /* copying the rotation axis and/or acceleration vector */
861  for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
862 
863  NNEW(xbounact,double,*nboun);
864 
865  if(*ithermal==1) NNEW(t1act,double,*nk);
866 
867  /* assigning the body forces to the elements */
868 
869  if(*nbody>0){
870  ifreebody=*ne+1;
871  NNEW(ipobody,ITG,2*ifreebody**nbody);
872  for(k=1;k<=*nbody;k++){
873  FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
874  iendset,ialset,&inewton,nset,&ifreebody,&k));
875  RENEW(ipobody,ITG,2*(*ne+ifreebody));
876  }
877  RENEW(ipobody,ITG,2*(ifreebody-1));
878  }
879 
880  NNEW(br,double,neq[1]); /* load rhs vector */
881  NNEW(bi,double,neq[1]); /* load rhs vector */
882 
883  if(iprescribedboundary){
884  NNEW(xbounr,double,*nboun);
885  NNEW(xbouni,double,*nboun);
886 
887  NNEW(fr,double,neq[1]); /* force corresponding to real particular solution */
888  NNEW(fi,double,neq[1]); /* force corresponding to imaginary particular solution */
889 
890  NNEW(ubr,double,neq[1]); /* real particular solution */
891  NNEW(ubi,double,neq[1]); /* imaginary particular solution */
892 
893  NNEW(mubr,double,neq[1]); /* mass times real particular solution */
894  NNEW(mubi,double,neq[1]); /* mass times imaginary particular solution */
895  }
896 
897  NNEW(bjr,double,nev); /* real response modal decomposition */
898  NNEW(bji,double,nev); /* imaginary response modal decomposition */
899 
900  NNEW(ampli,double,*nam); /* instantaneous amplitude */
901 
902  if(nherm==1){
903  NNEW(aa,double,nev); /* modal coefficients of the real loading */
904  NNEW(bb,double,nev); /* modal coefficients of the imaginary loading */
905  }else{
906  NNEW(aa,double,2*nev); /* modal coefficients of the real loading */
907  NNEW(bb,double,2*nev); /* modal coefficients of the imaginary loading */
908  }
909 
910  /* result fields */
911 
912  NNEW(vr,double,mt**nk);
913  NNEW(vi,double,mt**nk);
914 
915  if(iprescribedboundary){
916 
917  /* LU decomposition of the stiffness matrix */
918 
919  if(*isolver==0){
920 #ifdef SPOOLES
921  spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
922  &symmetryflag,&inputformat,&nzs[2]);
923 #else
924  printf("*ERROR in steadystate: the SPOOLES library is not linked\n\n");
925  FORTRAN(stop,());
926 #endif
927  }
928  else if(*isolver==4){
929 #ifdef SGI
930  token=1;
931  sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
932 #else
933  printf("*ERROR in steadystate: the SGI library is not linked\n\n");
934  FORTRAN(stop,());
935 #endif
936  }
937  else if(*isolver==5){
938 #ifdef TAUCS
939  tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
940 #else
941  printf("*ERROR in steadystate: the TAUCS library is not linked\n\n");
942  FORTRAN(stop,());
943 #endif
944  }
945  else if(*isolver==7){
946 #ifdef PARDISO
947  pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
948  &symmetryflag,&inputformat,jq,&nzs[2]);
949 #else
950  printf("*ERROR in steadystate: the PARDISO library is not linked\n\n");
951  FORTRAN(stop,());
952 #endif
953  }
954  }
955 
956  for(l=0;l<ndatatot;l=l+*jout){
957  for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
958  time=freq[l]/(2.*pi);
959  timem=-time;
960  ptime=time;
961 
962  /* calculating cc */
963 
964  if(dashpot){
965  NNEW(adc,double,neq[1]);
966  NNEW(auc,double,nzs[1]);
967  FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
968  ndirboun,xboun,nboun,
969  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
970  nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
971  adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
972  ikmpc,ilmpc,ikboun,ilboun,
973  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
974  ielorien,norien,orab,ntmat_,
975  t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
976  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
977  xstiff,npmat_,&dtime,matname,mi,ncmat_,
978  ttime,&time,istep,&iinc,ibody,clearini,&mortar,springarea,
979  pslavsurf,pmastsurf,&reltime,&nasym));
980 
981  /* zc = damping matrix * eigenmodes */
982 
983  NNEW(zc,double,neq[1]*nev);
984  for(i=0;i<nev;i++){
985  FORTRAN(op,(&neq[1],&z[(long long)i*neq[1]],&zc[i*neq[1]],adc,auc,
986  jq,irow));
987  }
988 
989  /* cc is the reduced damping matrix (damping matrix mapped onto
990  space spanned by eigenmodes) */
991 
992  for(i=0;i<nev*nev;i++){cc[i]=0.;}
993  for(i=0;i<nev;i++){
994  for(j=0;j<=i;j++){
995  for(k=0;k<neq[1];k++){
996  cc[i*nev+j]+=z[(long long)j*neq[1]+k]*zc[i*neq[1]+k];
997  }
998  }
999  }
1000 
1001  /* symmetric part of cc matrix */
1002 
1003  for(i=0;i<nev;i++){
1004  for(j=i;j<nev;j++){
1005  cc[i*nev+j]=cc[j*nev+i];
1006  }
1007  }
1008  SFREE(zc);SFREE(adc);SFREE(auc);
1009  }
1010 
1011  /* calculating the instantaneous loads (forces, surface loading,
1012  centrifugal and gravity loading or temperature) */
1013 
1014  FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,
1015  xloadold,xload,xloadact,iamload,nload,ibody,xbody,
1016  nbody,xbodyold,xbodyact,t1old,t1,t1act,iamt1,nk,amta,namta,
1017  nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
1018  xbounold,xboun,xbounact,iamboun,nboun,nodeboun,ndirboun,
1019  nodeforc,ndirforc,istep,&iinc,co,vold,itg,&ntg,amname,
1020  ikboun,ilboun,nelemload,sideload,mi,
1021  ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
1022  iendset,ialset,&ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
1023 
1024  /* real part of forces */
1025 
1026  for (i=0;i<*nforc;i++){
1027  xforcr[i]=xforcact[i]*(1-iphaseforc[i]);
1028  }
1029 
1030  for (i=0;i<*nload;i++){
1031  for(j=0;j<2;j++){
1032  xloadr[2*i+j]=xloadact[2*i+j]*(1-iphaseload[i]);
1033  }
1034  }
1035 
1036  for(i=0;i<*nbody;i++){
1037  for(j=0;j<7;j++){
1038  xbodyr[7*i+j]=xbodyact[7*i+j];
1039  }
1040  if(ibody[3*i+2]==2){
1041  xbodyr[7*i]=0.;
1042  }
1043  }
1044 
1045  /* calculating the instantaneous loading vector */
1046 
1047  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1048  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcr,
1049  nforc,nelemload,sideload,xloadr,nload,xbodyr,
1050  ipobody,nbody,cgr,br,nactdof,&neq[1],nmethod,
1051  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1052  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1053  t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1054  nplicon,plkcon,nplkcon,
1055  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1056  xbodyold,&reltime,veold,matname,mi,ikactmechr,
1057  &nactmechr,ielprop,prop));
1058 
1059  /* real modal coefficients */
1060 
1061  if(!iprescribedboundary){
1062 
1063  if(nherm==1){
1064  if(!cyclicsymmetry){
1065  for(i=0;i<nev;i++){
1066  i2=(long long)i*neq[1];
1067  aa[i]=0.;
1068  if(nactmechr<neq[1]/2){
1069  for(j=0;j<nactmechr;j++){
1070  aa[i]+=z[i2+ikactmechr[j]]*br[ikactmechr[j]];
1071  }
1072  }else{
1073  for(j=0;j<neq[1];j++){
1074  aa[i]+=z[i2+j]*br[j];
1075  }
1076  }
1077  }
1078  }else{
1079  for(i=0;i<nev;i++){aa[i]=0.;}
1080  for(j=0;j<nactmechr;j++){
1081  for(i=0;i<nev;i++){
1082  FORTRAN(nident,(izdof,&ikactmechr[j],&nzdof,&id));
1083  if(id!=0){
1084  if(izdof[id-1]==ikactmechr[j]){
1085  aa[i]+=z[(long long)i*nzdof+id-1]*br[ikactmechr[j]];
1086  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1087  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1088  }
1089  }
1090  }
1091  }else{
1092  if(!cyclicsymmetry){
1093  for(i=0;i<nev;i++){
1094  aa[2*i]=0.;aa[2*i+1]=0.;
1095  if(nactmechr<neq[1]/2){
1096  for(j=0;j<nactmechr;j++){
1097  aa[2*i]+=z[(long long)2*i*neq[1]+ikactmechr[j]]*br[ikactmechr[j]];
1098  aa[2*i+1]+=z[(long long)(2*i+1)*neq[1]+ikactmechr[j]]*br[ikactmechr[j]];
1099  }
1100  }else{
1101  for(j=0;j<neq[1];j++){
1102  aa[2*i]+=z[(long long)2*i*neq[1]+j]*br[j];
1103  aa[2*i+1]+=z[(long long)(2*i+1)*neq[1]+j]*br[j];
1104  }
1105  }
1106  }
1107  }else{
1108  for(i=0;i<nev;i++){aa[2*i]=0.;aa[2*i+1]=0.;}
1109  for(j=0;j<nactmechr;j++){
1110  for(i=0;i<nev;i++){
1111  FORTRAN(nident,(izdof,&ikactmechr[j],&nzdof,&id));
1112  if(id!=0){
1113  if(izdof[id-1]==ikactmechr[j]){
1114  aa[i]+=z[(long long)2*i*nzdof+id-1]*br[ikactmechr[j]];
1115  aa[i]+=z[(long long)(2*i+1)*nzdof+id-1]*br[ikactmechr[j]];
1116  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1117  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1118  }
1119  }
1120  }
1121  }
1122 
1123  }else{
1124 
1125  /* correction for nonzero SPC's */
1126 
1127  /* next statement makes sure that br is reset to zero at the
1128  start of rhs.f */
1129  nactmechr=neq[1];
1130 
1131  /* real part of boundary conditions */
1132 
1133  for (i=0;i<*nboun;i++){
1134  xbounr[i]=xbounact[i]*(1-iphaseboun[i]);
1135  }
1136 
1137  for(j=0;j<neq[1];j++){fr[j]=0.;ubr[j]=0.;}
1138  for(i=0;i<*nboun;i++){
1139  ic=neq[1]+i;
1140  for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
1141  ir=irow[j]-1;
1142  fr[ir]=fr[ir]-au[j]*xbounr[i];
1143  ubr[ir]=fr[ir];
1144  }
1145  }
1146  if(*isolver==0){
1147 #ifdef SPOOLES
1148  spooles_solve(ubr,&neq[1]);
1149 #endif
1150  }
1151  else if(*isolver==4){
1152 #ifdef SGI
1153  sgi_solve(ubr,token);
1154 #endif
1155  }
1156  else if(*isolver==5){
1157 #ifdef TAUCS
1158  tau_solve(ubr,&neq[1]);
1159 #endif
1160  }
1161  else if(*isolver==7){
1162 #ifdef PARDISO
1163  pardiso_solve(ubr,&neq[1],&symmetryflag);
1164 #endif
1165  }
1166  FORTRAN(op,(&neq[1],ubr,mubr,adb,aub,jq,irow));
1167  }
1168 
1169  /* imaginary part of forces */
1170 
1171  for (i=0;i<*nforc;i++){
1172  xforci[i]=xforcact[i]*iphaseforc[i];
1173  }
1174 
1175  for (i=0;i<*nload;i++){
1176  for(j=0;j<2;j++){
1177  xloadi[2*i+j]=xloadact[2*i+j]*iphaseload[i];
1178  }
1179  }
1180 
1181  for(i=0;i<*nbody;i++){
1182  for(j=0;j<7;j++){
1183  xbodyi[7*i+j]=xbodyact[7*i+j];
1184  }
1185  if(ibody[3*i+2]==1){
1186  xbodyi[7*i]=0.;
1187  }
1188  }
1189 
1190  /* calculating the instantaneous loading vector */
1191 
1192  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
1193  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforci,
1194  nforc,nelemload,sideload,xloadi,nload,xbodyi,
1195  ipobody,nbody,cgr,bi,nactdof,&neq[1],nmethod,
1196  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
1197  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
1198  t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
1199  nplicon,plkcon,nplkcon,
1200  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
1201  xbodyold,&reltime,veold,matname,mi,ikactmechi,
1202  &nactmechi,ielprop,prop));
1203 
1204  /* imaginary modal coefficients */
1205 
1206  if(!iprescribedboundary){
1207 
1208  if(nherm==1){
1209  if(!cyclicsymmetry){
1210  for(i=0;i<nev;i++){
1211  i2=(long long)i*neq[1];
1212  bb[i]=0.;
1213  if(nactmechi<neq[1]/2){
1214  for(j=0;j<nactmechi;j++){
1215  bb[i]+=z[i2+ikactmechi[j]]*bi[ikactmechi[j]];
1216  }
1217  }else{
1218  for(j=0;j<neq[1];j++){
1219  bb[i]+=z[i2+j]*bi[j];
1220  }
1221  }
1222  }
1223  }else{
1224  for(i=0;i<nev;i++){bb[i]=0.;}
1225  for(j=0;j<nactmechi;j++){
1226  for(i=0;i<nev;i++){
1227  FORTRAN(nident,(izdof,&ikactmechi[j],&nzdof,&id));
1228  if(id!=0){
1229  if(izdof[id-1]==ikactmechi[j]){
1230  bb[i]+=z[(long long)i*nzdof+id-1]*bi[ikactmechi[j]];
1231  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1232  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1233  }
1234  }
1235  }
1236  }else{
1237  if(!cyclicsymmetry){
1238  for(i=0;i<nev;i++){
1239  bb[2*i]=0.;bb[2*i+1]=0.;
1240  if(nactmechi<neq[1]/2){
1241  for(j=0;j<nactmechi;j++){
1242  bb[2*i]+=z[(long long)2*i*neq[1]+ikactmechi[j]]*bi[ikactmechi[j]];
1243  bb[2*i+1]+=z[(long long)(2*i+1)*neq[1]+ikactmechi[j]]*bi[ikactmechi[j]];
1244  }
1245  }else{
1246  for(j=0;j<neq[1];j++){
1247  bb[2*i]+=z[(long long)2*i*neq[1]+j]*bi[j];
1248  bb[2*i+1]+=z[(long long)(2*i+1)*neq[1]+j]*bi[j];
1249  }
1250  }
1251  }
1252  }else{
1253  for(i=0;i<nev;i++){bb[2*i]=0.;bb[2*i+1]=0.;}
1254  for(j=0;j<nactmechi;j++){
1255  for(i=0;i<nev;i++){
1256  FORTRAN(nident,(izdof,&ikactmechi[j],&nzdof,&id));
1257  if(id!=0){
1258  if(izdof[id-1]==ikactmechi[j]){
1259  bb[2*i]+=z[(long long)2*i*nzdof+id-1]*bi[ikactmechi[j]];
1260  bb[2*i+1]+=z[(long long)(2*i+1)*nzdof+id-1]*bi[ikactmechi[j]];
1261  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1262  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1263  }
1264  }
1265  }
1266  }
1267 
1268  }else{
1269 
1270  /* correction for nonzero SPC's */
1271 
1272  /* next statement makes sure that bi is reset to zero at the
1273  start of rhs.f */
1274 
1275  nactmechi=neq[1];
1276 
1277  /* imaginary part of boundary conditions */
1278 
1279  for (i=0;i<*nboun;i++){
1280  xbouni[i]=xbounact[i]*iphaseboun[i];
1281  }
1282 
1283  for(j=0;j<neq[1];j++){fi[j]=0.;ubi[j]=0.;}
1284  for(i=0;i<*nboun;i++){
1285  ic=neq[1]+i;
1286  for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
1287  ir=irow[j]-1;
1288  fi[ir]=fi[ir]-au[j]*xbouni[i];
1289  ubi[ir]=fi[ir];
1290  }
1291  }
1292  if(*isolver==0){
1293 #ifdef SPOOLES
1294  spooles_solve(ubi,&neq[1]);
1295 #endif
1296  }
1297  else if(*isolver==4){
1298 #ifdef SGI
1299  sgi_solve(ubi,token);
1300 #endif
1301  }
1302  else if(*isolver==5){
1303 #ifdef TAUCS
1304  tau_solve(ubi,&neq[1]);
1305 #endif
1306  }
1307  else if(*isolver==7){
1308 #ifdef PARDISO
1309  pardiso_solve(ubi,&neq[1],&symmetryflag);
1310 #endif
1311  }
1312  FORTRAN(op,(&neq[1],ubi,mubi,adb,aub,jq,irow));
1313 
1314  /* correction for prescribed boundary conditions */
1315 
1316  for(i=0;i<neq[1];i++){
1317  br[i]+=freq[l]*(freq[l]*mubr[i]+alpham*mubi[i]+betam*fi[i]);
1318  bi[i]+=freq[l]*(freq[l]*mubi[i]-alpham*mubr[i]-betam*fr[i]);
1319  }
1320 
1321  /* real and imaginary modal coefficients */
1322 
1323  for(i=0;i<nev;i++){
1324  aa[i]=0.;
1325  for(j=0;j<neq[1];j++){
1326  aa[i]+=z[(long long)i*neq[1]+j]*br[j];
1327  }
1328  }
1329 
1330  for(i=0;i<nev;i++){
1331  bb[i]=0.;
1332  for(j=0;j<neq[1];j++){
1333  bb[i]+=z[(long long)i*neq[1]+j]*bi[j];
1334  }
1335  }
1336 
1337  }
1338 
1339  /* calculating the modal coefficients */
1340 
1341  if(nherm==1){
1342  if(dashpot==0){
1343  for(i=0;i<nev;i++){
1344  dd=pow(d[i]-pow(freq[l],2),2)+
1345  pow(fric[i],2)*pow(freq[l],2);
1346  bjr[i]=(aa[i]*(d[i]-freq[l]*freq[l])+
1347  bb[i]*fric[i]*freq[l])/dd;
1348  bji[i]=(bb[i]*(d[i]-freq[l]*freq[l])-
1349  aa[i]*fric[i]*freq[l])/dd;
1350  }
1351  /* printf("old l=%" ITGFORMAT ",bjr=%f,bji=%f\n",l,bjr[0],bji[0]);*/
1352  }else{
1353  nev2=2*nev;
1354  NNEW(am,double,nev2*nev2);
1355  NNEW(bm,double,nev2);
1356  NNEW(ipiv,ITG,nev2);
1357 
1358  for(i=0;i<nev2;i++){
1359  for(j=0;j<nev2;j++){
1360  am[i*nev2+j]=0.;
1361  }
1362  bm[i]=0.;
1363  }
1364  for(i=0;i<nev;i++){
1365  am[i*nev2+i]=d[i]-freq[l]*freq[l];
1366  am[(i+nev)*nev2+i]=-fric[i]*freq[l];
1367  bm[i]=aa[i];
1368  am[i*nev2+nev+i]=-am[(i+nev)*nev2+i];
1369  am[(i+nev)*nev2+nev+i]=am[i*nev2+i];
1370  bm[nev+i]=bb[i];
1371  for(j=0;j<nev;j++){
1372  am[(j+nev)*nev2+i]=am[(j+nev)*nev2+i]
1373  -cc[i*nev+j]*freq[l];
1374  am[j*nev2+nev+i]=am[j*nev2+nev+i]
1375  +cc[i*nev+j]*freq[l];
1376  }
1377  }
1378 
1379  /* solving the system of equations */
1380 
1381  FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
1382  if(info!=0){
1383  printf("*ERROR in steadystate: fatal termination of dgesv\n");
1384  printf(" info=%" ITGFORMAT "\n",info);
1385  FORTRAN(stop,());
1386  }
1387 
1388  /* storing the solution in bjr and bji */
1389 
1390  for(i=0;i<nev;i++){
1391  bjr[i]=bm[i];
1392  bji[i]=bm[nev+i];
1393  }
1394 
1395  SFREE(am);SFREE(bm);SFREE(ipiv);
1396  }
1397  }else{
1398  nev2=2*nev;
1399  NNEW(am,double,nev2*nev2);
1400  NNEW(bm,double,nev2);
1401  NNEW(ipiv,ITG,nev2);
1402 
1403  NNEW(ax,double,nev);
1404  NNEW(bx,double,nev);
1405  for(i=0;i<nev;i++){
1406  ax[i]=-pow(freq[l],2)-freq[l]*fric[2*i+1]+d[2*i];
1407  bx[i]=-freq[l]*fric[2*i]-d[2*i+1];
1408  }
1409  for(i=0;i<nev;i++){
1410  for(j=0;j<nev;j++){
1411  am[j*nev2+i]=xmr[j*nev+i]*ax[j]+xmi[j*nev+i]*bx[j];
1412  am[(j+nev)*nev2+i]=xmr[j*nev+i]*bx[j]-xmi[j*nev+i]*bx[j];
1413  am[j*nev2+nev+i]=xmi[j*nev+i]*ax[j]-xmr[j*nev+i]*bx[j];
1414  am[(j+nev)*nev2+nev+i]=xmi[j*nev+i]*bx[j]+xmr[j*nev+i]*ax[j];
1415  }
1416  bm[i]=aa[i]-bb[nev+i];
1417  bm[nev+i]=bb[i]+aa[nev+i];
1418  }
1419  SFREE(ax);SFREE(bx);
1420 
1421  /* solving the system of equations */
1422 
1423  FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
1424  if(info!=0){
1425  printf("*ERROR in steadystate: fatal termination of dgesv\n");
1426  printf(" info=%" ITGFORMAT "\n",info);
1427  FORTRAN(stop,());
1428  }
1429 
1430  /* storing the solution in bjr and bji */
1431 
1432  for(i=0;i<nev;i++){
1433  bjr[i]=bm[i];
1434  bji[i]=bm[nev+i];
1435  }
1436 
1437  SFREE(am);SFREE(bm);SFREE(ipiv);
1438  }
1439 
1440  /* storing the participation factors */
1441 
1442  if(nherm==1){
1443  NNEW(eig,double,nev);
1444  for(i=0;i<nev;i++){
1445  eig[i]=sqrt(d[i]);
1446  }
1447  }else{
1448 
1449  NNEW(eig,double,2*nev);
1450  for(i=0;i<nev;i++){
1451  eig[2*i]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])+d[2*i])
1452  /sqrt(2.);
1453  eig[2*i+1]=sqrt(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])-d[2*i])
1454  /sqrt(2.);
1455  if(d[2*i+1]<0.) eig[2*i+1]=-eig[2*i+1];
1456  }
1457  }
1458 
1459  mode=0;
1460  FORTRAN(writepf,(eig,bjr,bji,&time,&nev,&mode,&nherm));
1461  SFREE(eig);
1462 
1463  /* calculating the real response */
1464 
1465  if(iprescribedboundary){
1466  if(nmdnode==0){
1467  memcpy(&br[0],&ubr[0],sizeof(double)*neq[1]);
1468  }else{
1469  for(i=0;i<nmddof;i++){
1470  br[imddof[i]]=ubr[imddof[i]];
1471  }
1472  }
1473  }
1474  else{
1475  if(nmdnode==0){
1476  DMEMSET(br,0,neq[1],0.);
1477  }else{
1478  for(i=0;i<nmddof;i++){
1479  br[imddof[i]]=0.;
1480  }
1481  }
1482  }
1483 
1484  if(!cyclicsymmetry){
1485  if(nmdnode==0){
1486  for(i=0;i<neq[1];i++){
1487  for(j=0;j<nev;j++){
1488  br[i]+=bjr[j]*z[(long long)j*neq[1]+i];
1489  }
1490  }
1491  }else{
1492  for(i=0;i<nmddof;i++){
1493  for(j=0;j<nev;j++){
1494  br[imddof[i]]+=bjr[j]*z[(long long)j*neq[1]+imddof[i]];
1495  }
1496  }
1497  }
1498  }else{
1499  for(i=0;i<nmddof;i++){
1500  FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1501  if(id!=0){
1502  if(izdof[id-1]==imddof[i]){
1503  for(j=0;j<nev;j++){
1504  br[imddof[i]]+=bjr[j]*z[(long long)j*nzdof+id-1];
1505  }
1506  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1507  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1508  }
1509  }
1510 
1511  if(nmdnode==0){
1512  DMEMSET(vr,0,mt**nk,0.);
1513  }else{
1514  for(jj=0;jj<nmdnode;jj++){
1515  i=imdnode[jj]-1;
1516  for(j=1;j<4;j++){
1517  vr[mt*i+j]=0.;
1518  }
1519  }
1520  }
1521 
1522  /* calculating the imaginary response */
1523 
1524  if(iprescribedboundary){
1525  if(nmdnode==0){
1526  memcpy(&bi[0],&ubi[0],sizeof(double)*neq[1]);
1527  }else{
1528  for(i=0;i<nmddof;i++){
1529  bi[imddof[i]]=ubi[imddof[i]];
1530  }
1531  }
1532  }
1533  else{
1534  if(nmdnode==0){
1535  DMEMSET(bi,0,neq[1],0.);
1536  }else{
1537  for(i=0;i<nmddof;i++){
1538  bi[imddof[i]]=0.;
1539  }
1540  }
1541  }
1542 
1543  if(!cyclicsymmetry){
1544  if(nmdnode==0){
1545  for(i=0;i<neq[1];i++){
1546  for(j=0;j<nev;j++){
1547  bi[i]+=bji[j]*z[(long long)j*neq[1]+i];
1548  }
1549  }
1550  }else{
1551  for(i=0;i<nmddof;i++){
1552  for(j=0;j<nev;j++){
1553  bi[imddof[i]]+=bji[j]*z[(long long)j*neq[1]+imddof[i]];
1554  }
1555  }
1556  }
1557  }else{
1558  for(i=0;i<nmddof;i++){
1559  FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
1560  if(id!=0){
1561  if(izdof[id-1]==imddof[i]){
1562  for(j=0;j<nev;j++){
1563  bi[imddof[i]]+=bji[j]*z[(long long)j*nzdof+id-1];
1564  }
1565  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1566  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
1567  }
1568  }
1569 
1570 
1571  if(nmdnode==0){
1572  DMEMSET(vi,0,mt**nk,0.);
1573  }else{
1574  for(jj=0;jj<nmdnode;jj++){
1575  i=imdnode[jj]-1;
1576  for(j=1;j<4;j++){
1577  vi[mt*i+j]=0.;
1578  }
1579  }
1580  }
1581 
1582  /* real response */
1583 
1584  if(iprescribedboundary){
1585 
1586  /* calculating displacements/temperatures */
1587 
1588  FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
1589  ndirboun,xbounr,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1590  br,bi,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1591  imdmpc,&nmdmpc,nmethod,&timem));
1592 
1593  results(co,nk,kon,ipkon,lakon,ne,vr,stnr,inum,
1594  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1595  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1596  ithermal,prestr,iprestr,filab,eme,emn,een,
1597  iperturb,f,fn,nactdof,&iout,qa,
1598  vold,br,nodeboun,ndirboun,xbounr,nboun,
1599  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1600  veold,accold,&bet,&gam,&dtime,&time,&xnull,
1601  plicon,nplicon,plkcon,nplkcon,
1602  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1603  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1604  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1605  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1606  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1607  ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1608  xforc,nforc,thicke,shcon,nshcon,
1609  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1610  &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
1611  islavsurf,ielprop,prop);}
1612  else{
1613 
1614  /* calculating displacements/temperatures */
1615 
1616  FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
1617  ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1618  br,bi,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1619  imdmpc,&nmdmpc,nmethod,&timem));
1620 
1621  results(co,nk,kon,ipkon,lakon,ne,vr,stnr,inum,
1622  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1623  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1624  ithermal,prestr,iprestr,filab,eme,emn,een,
1625  iperturb,f,fn,nactdof,&iout,qa,
1626  vold,br,nodeboun,ndirboun,xbounact,nboun,
1627  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1628  veold,accold,&bet,&gam,&dtime,&time,&xnull,
1629  plicon,nplicon,plkcon,nplkcon,
1630  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1631  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1632  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1633  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1634  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1635  ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1636  xforc,nforc,thicke,shcon,nshcon,
1637  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1638  &mortar,islavact,cdn,islavnode,nslavnode,&ntie,
1639  clearini,islavsurf,ielprop,prop);
1640 
1641  if(nmdnode==0){
1642  DMEMSET(br,0,neq[1],0.);
1643  }else{
1644  for(i=0;i<nmddof;i++){
1645  br[imddof[i]]=0.;
1646  }
1647  }
1648  }
1649 
1650  (*kode)++;
1651 
1652  mode=-1;
1653  if(strcmp1(&filab[1044],"ZZS")==0){
1654  NNEW(neigh,ITG,40**ne);
1655  NNEW(ipneigh,ITG,*nk);
1656  }
1657 
1658  frd(co,&nkg,kon,ipkon,lakon,&neg,vr,stnr,inum,nmethod,
1659  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1660  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1661  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1662  mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&neg,
1663  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1664  thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1665 
1666  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1667 
1668  /* imaginary response */
1669 
1670  if(iprescribedboundary){
1671 
1672  /* calculating displacements/temperatures */
1673 
1674  FORTRAN(dynresults,(nk,vi,ithermal,nactdof,vold,nodeboun,
1675  ndirboun,xbouni,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1676  bi,br,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1677  imdmpc,&nmdmpc,nmethod,&time));
1678 
1679  results(co,nk,kon,ipkon,lakon,ne,vi,stni,inum,
1680  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1681  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1682  ithermal,prestr,iprestr,filab,eme,emn,een,
1683  iperturb,f,fn,nactdof,&iout,qa,
1684  vold,bi,nodeboun,ndirboun,xbouni,nboun,
1685  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1686  veold,accold,&bet,&gam,&dtime,&time,&xnull,
1687  plicon,nplicon,plkcon,nplkcon,
1688  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1689  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1690  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1691  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1692  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1693  ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1694  xforc,nforc,thicke,shcon,nshcon,
1695  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1696  &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
1697  islavsurf,ielprop,prop);}
1698  else{
1699 
1700  /* calculating displacements/temperatures */
1701 
1702  FORTRAN(dynresults,(nk,vi,ithermal,nactdof,vold,nodeboun,
1703  ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,
1704  bi,br,veold,&dtime,mi,imdnode,&nmdnode,imdboun,&nmdboun,
1705  imdmpc,&nmdmpc,nmethod,&time));
1706 
1707  results(co,nk,kon,ipkon,lakon,ne,vi,stni,inum,
1708  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
1709  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
1710  ithermal,prestr,iprestr,filab,eme,emn,een,
1711  iperturb,f,fn,nactdof,&iout,qa,
1712  vold,bi,nodeboun,ndirboun,xbounact,nboun,
1713  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1714  veold,accold,&bet,&gam,&dtime,&time,&xnull,
1715  plicon,nplicon,plkcon,nplkcon,
1716  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1717  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
1718  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
1719  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
1720  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,
1721  ikmpc,ilmpc,istep,&iinc,springarea,&reltime,&ne0,
1722  xforc,nforc,thicke,shcon,nshcon,
1723  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1724  &mortar,islavact,cdn,islavnode,nslavnode,&ntie,
1725  clearini,islavsurf,ielprop,prop);
1726 
1727  if(nmdnode==0){
1728  DMEMSET(bi,0,neq[1],0.);
1729  }else{
1730  for(i=0;i<nmddof;i++){
1731  bi[imddof[i]]=0.;
1732  }
1733  }
1734  }
1735 
1736  /* calculating the magnitude and phase */
1737 
1738  if(strcmp1(&filab[870],"PU")==0){
1739 
1740  constant=180./pi;
1741  NNEW(va,double,mt**nk);
1742  NNEW(vp,double,mt**nk);
1743 
1744  if(*ithermal<=1){
1745  if(nmdnode==0){
1746  for(i=0;i<*nk;i++){
1747  for(j=1;j<4;j++){
1748  vreal=vr[mt*i+j];
1749  va[mt*i+j]=sqrt(vr[mt*i+j]*vr[mt*i+j]+vi[mt*i+j]*vi[mt*i+j]);
1750  if(fabs(vreal)<1.e-10){
1751  if(vi[mt*i+j]>0.){vp[mt*i+j]=90.;}
1752  else{vp[mt*i+j]=-90.;}
1753  }
1754  else{
1755  vp[mt*i+j]=atan(vi[mt*i+j]/vreal)*constant;
1756  if(vreal<0.) vp[mt*i+j]+=180.;
1757  }
1758  }
1759  }
1760  }else{
1761  for(jj=0;jj<nmdnode;jj++){
1762  i=imdnode[jj]-1;
1763  for(j=1;j<4;j++){
1764  vreal=vr[mt*i+j];
1765  va[mt*i+j]=sqrt(vr[mt*i+j]*vr[mt*i+j]+vi[mt*i+j]*vi[mt*i+j]);
1766  if(fabs(vreal)<1.e-10){
1767  if(vi[mt*i+j]>0.){vp[mt*i+j]=90.;}
1768  else{vp[mt*i+j]=-90.;}
1769  }
1770  else{
1771  vp[mt*i+j]=atan(vi[mt*i+j]/vreal)*constant;
1772  if(vreal<0.) vp[mt*i+j]+=180.;
1773  }
1774  }
1775  }
1776  }
1777  }
1778  else{
1779  if(nmdnode==0){
1780  for(i=0;i<*nk;i++){
1781  vreal=vr[mt*i];
1782  va[mt*i]=sqrt(vr[mt*i]*vr[mt*i]+vi[mt*i]*vi[mt*i]);
1783  if(fabs(vreal)<1.e-10){
1784  if(vi[mt*i]>0){vp[mt*i]=90.;}
1785  else{vp[mt*i]=-90.;}
1786  }
1787  else{
1788  vp[mt*i]=atan(vi[mt*i]/vreal)*constant;
1789  if(vreal<0.) vp[mt*i]+=180.;
1790  }
1791  }
1792  }else{
1793  for(jj=0;jj<nmdnode;jj++){
1794  i=imdnode[jj]-1;
1795  vreal=vr[mt*i];
1796  va[mt*i]=sqrt(vr[mt*i]*vr[mt*i]+vi[mt*i]*vi[mt*i]);
1797  if(fabs(vreal)<1.e-10){
1798  if(vi[mt*i]>0){vp[mt*i]=90.;}
1799  else{vp[mt*i]=-90.;}
1800  }
1801  else{
1802  vp[mt*i]=atan(vi[mt*i]/vreal)*constant;
1803  if(vreal<0.) vp[mt*i]+=180.;
1804  }
1805  }
1806  }
1807  }
1808  }
1809 
1810  if(strcmp1(&filab[1479],"PHS")==0){
1811 
1812  constant=180./pi;
1813  NNEW(stna,double,6**nk);
1814  NNEW(stnp,double,6**nk);
1815 
1816  if(*ithermal<=1){
1817  if(nmdnode==0){
1818  for(i=0;i<*nk;i++){
1819  for(j=0;j<6;j++){
1820  vreal=stnr[6*i+j];
1821  stna[6*i+j]=sqrt(stnr[6*i+j]*stnr[6*i+j]+stni[6*i+j]*stni[6*i+j]);
1822  if(fabs(vreal)<1.e-10){
1823  if(stni[6*i+j]>0.){stnp[6*i+j]=90.;}
1824  else{stnp[6*i+j]=-90.;}
1825  }
1826  else{
1827  stnp[6*i+j]=atan(stni[6*i+j]/vreal)*constant;
1828  if(vreal<0.) stnp[6*i+j]+=180.;
1829  }
1830  }
1831  }
1832  }else{
1833  for(jj=0;jj<nmdnode;jj++){
1834  i=imdnode[jj]-1;
1835  for(j=0;j<6;j++){
1836  vreal=stnr[6*i+j];
1837  stna[6*i+j]=sqrt(stnr[6*i+j]*stnr[6*i+j]+stni[6*i+j]*stni[6*i+j]);
1838  if(fabs(vreal)<1.e-10){
1839  if(stni[6*i+j]>0.){stnp[6*i+j]=90.;}
1840  else{stnp[6*i+j]=-90.;}
1841  }
1842  else{
1843  stnp[6*i+j]=atan(stni[6*i+j]/vreal)*constant;
1844  if(vreal<0.) stnp[6*i+j]+=180.;
1845  }
1846  }
1847  }
1848  }
1849  }
1850  }
1851 
1852  (*kode)++;
1853  mode=0;
1854 
1855  if(strcmp1(&filab[1044],"ZZS")==0){
1856  NNEW(neigh,ITG,40**ne);
1857  NNEW(ipneigh,ITG,*nk);
1858  }
1859 
1860  frd(co,&nkg,kon,ipkon,lakon,&neg,vi,stni,inum,nmethod,
1861  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
1862  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
1863  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1864  mi,stx,va,vp,stna,stnp,vmax,stnmax,&ngraph,veold,ener,&neg,
1865  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1866  thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1867 
1868  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1869 
1870  SFREE(va);SFREE(vp);SFREE(stna);SFREE(stnp);
1871 
1872  }
1873 
1874  /* restoring the imaginary loading */
1875 
1876  SFREE(iphaseforc);SFREE(xforcr);SFREE(xforci);
1877 
1878  SFREE(iphaseload);SFREE(xloadr);SFREE(xloadi);
1879 
1880  SFREE(xbodyr);SFREE(xbodyi);
1881 
1882  if(iprescribedboundary){
1883  for (i=0;i<*nboun;i++){
1884  if(iphaseboun[i]==1){
1885  nodeboun[i]=nodeboun[i]+*nk;
1886  }
1887  }
1888  SFREE(iphaseboun);
1889  }
1890 
1891  /* updating the loading at the end of the step;
1892  important in case the amplitude at the end of the step
1893  is not equal to one */
1894 
1895  for(k=0;k<*nboun;++k){xboun[k]=xbounact[k];}
1896  for(k=0;k<*nforc;++k){xforc[k]=xforcact[k];}
1897  for(k=0;k<2**nload;++k){xload[k]=xloadact[k];}
1898  for(k=0;k<7**nbody;k=k+7){xbody[k]=xbodyact[k];}
1899  if(*ithermal==1){
1900  for(k=0;k<*nk;++k){t1[k]=t1act[k];}
1901  }
1902 
1903  SFREE(br);SFREE(bi);SFREE(bjr);SFREE(bji),SFREE(freq);
1904  SFREE(xforcact);SFREE(xloadact);SFREE(xbounact);SFREE(aa);SFREE(bb);
1905  SFREE(ampli);SFREE(xbodyact);SFREE(vr);SFREE(vi);if(*nbody>0) SFREE(ipobody);
1906 
1907  if(*ithermal==1) SFREE(t1act);
1908 
1909  if(iprescribedboundary){
1910  if(*isolver==0){
1911 #ifdef SPOOLES
1912  spooles_cleanup();
1913 #endif
1914  }
1915  else if(*isolver==4){
1916 #ifdef SGI
1917  sgi_cleanup(token);
1918 #endif
1919  }
1920  else if(*isolver==5){
1921 #ifdef TAUCS
1922  tau_cleanup();
1923 #endif
1924  }
1925  else if(*isolver==7){
1926 #ifdef PARDISO
1927  pardiso_cleanup(&neq[1],&symmetryflag);
1928 #endif
1929  }
1930  SFREE(xbounr);SFREE(xbouni);SFREE(fr);SFREE(fi);SFREE(ubr);SFREE(ubi);
1931  SFREE(mubr);SFREE(mubi);
1932  }
1933 
1934  SFREE(ikactmechr);SFREE(ikactmechi);
1935 
1936  if(intpointvar==1){
1937  SFREE(fn);
1938  SFREE(stnr);SFREE(stni);SFREE(stx);SFREE(eei);
1939 
1940  if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
1941 
1942  if(strcmp1(&filab[261],"E ")==0) SFREE(een);
1943  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1944  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
1945 
1946  if(*nener==1){SFREE(stiini);SFREE(enerini);}
1947  }
1948 
1949  }else{
1950 
1951  /* steady state response to a nonharmonic periodic loading */
1952 
1953  NNEW(ikactmech,ITG,neq[1]);
1954  nactmech=0;
1955 
1956  NNEW(xforcact,double,nfour**nforc);
1957  NNEW(xloadact,double,nfour*2**nload);
1958  NNEW(xbodyact,double,nfour*7**nbody);
1959  NNEW(xbounact,double,nfour**nboun);
1960  NNEW(xbounacttime,double,nfour**nboun);
1961  if(*ithermal==1) NNEW(t1act,double,*nk);
1962 
1963  NNEW(r,double,nfour);
1964  NNEW(wsave,double,2*nfour);
1965  NNEW(isave,ITG,15);
1966 
1967  /* check for nonzero SPC's */
1968 
1969  iprescribedboundary=0;
1970  for(i=0;i<*nboun;i++){
1971  if(fabs(xboun[i])>1.e-10){
1972  iprescribedboundary=1;
1973  break;
1974  }
1975  }
1976 
1977  if((iprescribedboundary)&&(cyclicsymmetry)){
1978  printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with cyclic symmetry\n");
1979  FORTRAN(stop,());
1980  }
1981 
1982  /* calculating the damping coefficients = friction coefficient*2*eigenvalue */
1983 
1984  if(xmodal[10]<0){
1985  for(i=0;i<nev;i++){
1986  if(sqrt(d[i])>(1.e-10)){
1987  fric[i]=(alpham+betam*d[i]);
1988  }
1989  else {
1990  printf("*WARNING in steadystate: one of the frequencies is zero\n");
1991  printf(" no Rayleigh mass damping allowed\n");
1992  fric[i]=0.;
1993  }
1994  }
1995  }
1996  else{
1997  if(iprescribedboundary){
1998  printf("*ERROR in steadystate: prescribed boundaries are not allowed in combination with direct modal damping\n");
1999  FORTRAN(stop,());
2000  }
2001 
2002  /*copy the damping coefficients for every eigenfrequencie from xmodal[11....] */
2003  if(nev<(ITG)xmodal[10]){
2004  imax=nev;
2005  printf("*WARNING in steadystate: too many modal damping coefficients applied\n");
2006  printf(" damping coefficients corresponding to nonexisting eigenvalues are ignored\n");
2007  }
2008  else{
2009  imax=(ITG)xmodal[10];
2010  }
2011  for(i=0; i<imax; i++){
2012  fric[i]=2.*sqrt(d[i])*xmodal[11+i];
2013  }
2014 
2015  }
2016 
2017  /* determining the load time history */
2018 
2019  NNEW(ampli,double,*nam); /* instantaneous amplitude */
2020 
2021  for(l=0;l<nfour;l++){
2022 
2023  time=tmin+(tmax-tmin)*(double)l/(double)nfour;
2024 
2025  FORTRAN(tempload,(xforcold,xforc,&xforcact[l**nforc],iamforc,nforc,
2026  xloadold,xload,&xloadact[l*2**nload],iamload,nload,ibody,xbody,
2027  nbody,xbodyold,&xbodyact[l*7**nbody],t1old,t1,t1act,
2028  iamt1,nk,amta,namta,nam,ampli,&time,&reltime,ttime,&dtime,
2029  ithermal,nmethod,xbounold,xboun,&xbounact[l**nboun],iamboun,nboun,
2030  nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,co,vold,itg,&ntg,
2031  amname,ikboun,ilboun,nelemload,sideload,mi,
2032  ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
2033  iendset,ialset,&ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc));
2034 
2035  }
2036 
2037  SFREE(ampli);
2038 
2039  for(i=0;i<l**nboun;i++){xbounacttime[i]=xbounact[i];}
2040 
2041  /* determining the load frequency history:
2042  frequency transform of the load time history */
2043 
2044  for(i=0;i<*nforc;i++){
2045  for(l=0;l<nfour;l++){
2046  r[l]=xforcact[l**nforc+i];
2047  }
2048  FORTRAN(drffti,(&nfour,wsave,isave));
2049  FORTRAN(drfftf,(&nfour,r,wsave,isave));
2050  for(l=0;l<nfour;l++){
2051  xforcact[l**nforc+i]=r[l]/nfour*2.;
2052  }
2053  xforcact[i]=xforcact[i]/2.;
2054  }
2055 
2056  for(i=0;i<*nload;i++){
2057  for(l=0;l<nfour;l++){
2058  r[l]=xloadact[l*2**nload+2*i];
2059  }
2060  FORTRAN(drffti,(&nfour,wsave,isave));
2061  FORTRAN(drfftf,(&nfour,r,wsave,isave));
2062  for(l=0;l<nfour;l++){
2063  xloadact[l*2**nload+2*i]=r[l]/nfour*2.;
2064  }
2065  xloadact[2*i]=xloadact[2*i]/2.;
2066  }
2067 
2068  for(i=0;i<*nbody;i++){
2069  for(l=0;l<nfour;l++){
2070  r[l]=xbodyact[l**nbody+7*i];
2071  }
2072  FORTRAN(drffti,(&nfour,wsave,isave));
2073  FORTRAN(drfftf,(&nfour,r,wsave,isave));
2074  for(l=0;l<nfour;l++){
2075  xbodyact[l**nbody+7*i]=r[l]/nfour*2.;
2076  }
2077  xbodyact[7*i]=xbodyact[7*i]/2.;
2078  }
2079 
2080  if(iprescribedboundary){
2081  for(i=0;i<*nboun;i++){
2082  for(l=0;l<nfour;l++){
2083  r[l]=xbounact[l**nboun+i];
2084  }
2085  FORTRAN(drffti,(&nfour,wsave,isave));
2086  FORTRAN(drfftf,(&nfour,r,wsave,isave));
2087  for(l=0;l<nfour;l++){
2088  xbounact[l**nboun+i]=r[l]/nfour*2.;
2089  }
2090  xbounact[i]=xbounact[i]/2.;
2091  }
2092 
2093  /* LU decomposition of the stiffness matrix */
2094 
2095  if(*isolver==0){
2096 #ifdef SPOOLES
2097  spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
2098  &symmetryflag,&inputformat,&nzs[2]);
2099 #else
2100  printf("*ERROR in steadystate: the SPOOLES library is not linked\n\n");
2101  FORTRAN(stop,());
2102 #endif
2103  }
2104  else if(*isolver==4){
2105 #ifdef SGI
2106  token=1;
2107  sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
2108 #else
2109  printf("*ERROR in steadystate: the SGI library is not linked\n\n");
2110  FORTRAN(stop,());
2111 #endif
2112  }
2113  else if(*isolver==5){
2114 #ifdef TAUCS
2115  tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
2116 #else
2117  printf("*ERROR in steadystate: the TAUCS library is not linked\n\n");
2118  FORTRAN(stop,());
2119 #endif
2120  }
2121  else if(*isolver==7){
2122 #ifdef PARDISO
2123  pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
2124  &symmetryflag,&inputformat,jq,&nzs[2]);
2125 #else
2126  printf("*ERROR in steadystate: the PARDISO library is not linked\n\n");
2127  FORTRAN(stop,());
2128 #endif
2129  }
2130 
2131  }
2132 
2133  SFREE(r);SFREE(wsave);SFREE(isave);
2134 
2135  /* determining the frequency data points */
2136 
2137  NNEW(freqnh,double,ndata*(nev+1));
2138 
2139  ndatatot=0.;
2140  freqnh[0]=fmin;
2141  if(fabs(fmax-fmin)<1.e-10){
2142  ndatatot=1;
2143  }else{
2144 
2145  /* copy the eigenvalues and sort them in ascending order
2146  (important for values from distinct nodal diameters */
2147 
2148  NNEW(e,double,nev);
2149  for(i=0;i<nev;i++){e[i]=sqrt(d[i]);}
2150  FORTRAN(dsort,(e,&idummy,&nev,&iflag));
2151 
2152  for(i=0;i<nev;i++){
2153  if(i!=0){
2154  if(fabs(e[i]-e[i-1])<1.e-5){continue;}
2155  }
2156  if(e[i]>=fmin){
2157  if(e[i]<=fmax){
2158  for(j=1;j<ndata;j++){
2159  y=-1.+2.*j/((double)(ndata-1));
2160  if(fabs(y)<1.e-10){freqnh[ndatatot+j]=
2161  (freqnh[ndatatot]+e[i])/2.;}
2162  else{
2163  freqnh[ndatatot+j]=(freqnh[ndatatot]+e[i])/2.+
2164  (e[i]-freqnh[ndatatot])*pow(fabs(y),1./bias)
2165  *y/(2.*fabs(y));
2166  }
2167  }
2168  ndatatot+=ndata-1;
2169  }
2170  else{break;}
2171  }
2172  }
2173  SFREE(e);
2174  for(j=1;j<ndata;j++){
2175  y=-1.+2.*j/((double)(ndata-1));
2176  if(fabs(y)<1.e-10){freqnh[ndatatot+j]=
2177  (freqnh[ndatatot]+fmax)/2.;}
2178  else{
2179  freqnh[ndatatot+j]=(freqnh[ndatatot]+fmax)/2.+
2180  (fmax-freqnh[ndatatot])*pow(fabs(y),1./bias)*
2181  y/(2.*fabs(y));
2182  }
2183  }
2184  ndatatot+=ndata;
2185  }
2186  RENEW(freqnh,double,ndatatot);
2187 
2188  for(ii=0;ii<ndatatot;ii++){
2189  for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
2190 
2191  sprintf(description,"%12f",freqnh[ii]/(2.*pi));
2192 
2193  NNEW(xforcr,double,*nforc);
2194  NNEW(xloadr,double,2**nload);
2195  NNEW(xbodyr,double,7**nbody);
2196  for(k=0;k<7**nbody;k++){xbodyr[k]=xbody[k];}
2197  if(iprescribedboundary){
2198  NNEW(xbounr,double,*nboun);
2199  NNEW(fr,double,neq[1]); /* force corresponding to real particular solution */
2200  NNEW(ubr,double,neq[1]); /* real particular solution */
2201  NNEW(mubr,double,neq[1]); /* mass times real particular solution */
2202  }
2203 
2204  /* assigning the body forces to the elements */
2205 
2206  if(*nbody>0){
2207  ifreebody=*ne+1;
2208  NNEW(ipobody,ITG,2*ifreebody**nbody);
2209  for(k=1;k<=*nbody;k++){
2210  FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
2211  iendset,ialset,&inewton,nset,&ifreebody,&k));
2212  RENEW(ipobody,ITG,2*(*ne+ifreebody));
2213  }
2214  RENEW(ipobody,ITG,2*(ifreebody-1));
2215  }
2216 
2217  NNEW(br,double,neq[1]); /* load rhs vector (real part) */
2218  NNEW(bi,double,neq[1]); /* load rhs vector (imaginary part) */
2219  NNEW(btot,double,nfour*neq[1]);
2220  NNEW(bp,double,nfour*neq[1]);
2221 
2222  NNEW(bjr,double,nev); /* real response modal decomposition */
2223  NNEW(bji,double,nev); /* imaginary response modal decomposition */
2224 
2225  NNEW(aa,double,nev); /* modal coefficients of the real loading */
2226  NNEW(bb,double,nev); /* modal coefficients of the imaginary loading */
2227 
2228  /* loop over all Fourier frequencies */
2229 
2230  NNEW(freq,double,nfour);
2231 
2232  for(l=0;l<nfour;l++){
2233 
2234  /* frequency */
2235 
2236  freq[l]=freqnh[ii]*floor((l+1.)/2.+0.1);
2237 
2238  /* calculating cc */
2239 
2240  if(dashpot){
2241  NNEW(adc,double,neq[1]);
2242  NNEW(auc,double,nzs[1]);
2243  FORTRAN(mafilldm,(co,nk,kon,ipkon,lakon,ne,nodeboun,
2244  ndirboun,xboun,nboun,
2245  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
2246  nforc,nelemload,sideload,xload,nload,xbody,ipobody,
2247  nbody,cgr,adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
2248  ikmpc,ilmpc,ikboun,ilboun,
2249  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
2250  ielorien,norien,orab,ntmat_,
2251  t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
2252  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
2253  xstiff,npmat_,&dtime,matname,mi,ncmat_,
2254  ttime,&freq[l],istep,&iinc,ibody,clearini,&mortar,springarea,
2255  pslavsurf,pmastsurf,&reltime,&nasym));
2256 
2257  /* zc = damping matrix * eigenmodes */
2258 
2259  NNEW(zc,double,neq[1]*nev);
2260  for(i=0;i<nev;i++){
2261  FORTRAN(op,(&neq[1],&z[(long long)i*neq[1]],&zc[i*neq[1]],
2262  adc,auc,jq,irow));
2263  }
2264 
2265  /* cc is the reduced damping matrix (damping matrix mapped onto
2266  space spanned by eigenmodes) */
2267 
2268  for(i=0;i<nev*nev;i++){cc[i]=0.;}
2269  for(i=0;i<nev;i++){
2270  for(j=0;j<=i;j++){
2271  for(k=0;k<neq[1];k++){
2272  cc[i*nev+j]+=z[(long long)j*neq[1]+k]*zc[i*neq[1]+k];
2273  }
2274  }
2275  }
2276 
2277  /* symmetric part of cc matrix */
2278 
2279  for(i=0;i<nev;i++){
2280  for(j=i;j<nev;j++){
2281  cc[i*nev+j]=cc[j*nev+i];
2282  }
2283  }
2284  SFREE(zc);SFREE(adc);SFREE(auc);
2285  }
2286 
2287  /* loading for this frequency */
2288 
2289  for(i=0;i<*nforc;i++){
2290  xforcr[i]=xforcact[l**nforc+i];
2291  }
2292 
2293  for(i=0;i<*nload;i++){
2294  xloadr[2*i]=xloadact[l*2**nload+2*i];
2295  }
2296 
2297  for(i=0;i<*nbody;i++){
2298  xbodyr[7*i]=xbodyact[l**nbody+7*i];
2299  }
2300 
2301  /* calculating the instantaneous loading vector */
2302 
2303  FORTRAN(rhs,(co,nk,kon,ipkon,lakon,ne,
2304  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforcr,
2305  nforc,nelemload,sideload,xloadr,nload,xbodyr,
2306  ipobody,nbody,cgr,br,nactdof,&neq[1],nmethod,
2307  ikmpc,ilmpc,elcon,nelcon,rhcon,nrhcon,
2308  alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
2309  t0,t1act,ithermal,iprestr,vold,iperturb,iexpl,plicon,
2310  nplicon,plkcon,nplkcon,
2311  npmat_,ttime,&time,istep,&iinc,&dtime,physcon,ibody,
2312  xbodyold,&reltime,veold,matname,mi,ikactmech,&nactmech,
2313  ielprop,prop));
2314 
2315  /* real modal coefficients */
2316 
2317  if(!iprescribedboundary){
2318  if(!cyclicsymmetry){
2319  for(i=0;i<nev;i++){
2320  i2=(long long)i*neq[1];
2321  aa[i]=0.;
2322  if(nactmech<neq[1]/2){
2323  for(j=0;j<nactmech;j++){
2324  aa[i]+=z[i2+ikactmech[j]]*br[ikactmech[j]];
2325  }
2326  }else{
2327  for(j=0;j<neq[1];j++){
2328  aa[i]+=z[i2+j]*br[j];
2329  }
2330  }
2331  }
2332  }else{
2333  for(i=0;i<nev;i++){aa[i]=0.;}
2334  for(j=0;j<nactmech;j++){
2335  for(i=0;i<nev;i++){
2336  FORTRAN(nident,(izdof,&ikactmech[j],&nzdof,&id));
2337  if(id!=0){
2338  if(izdof[id-1]==ikactmech[j]){
2339  aa[i]+=z[(long long)i*nzdof+id-1]*br[ikactmech[j]];
2340  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2341  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2342  }
2343  }
2344  }
2345 
2346  /* imaginary modal coefficients */
2347 
2348  for(i=0;i<nev;i++){
2349  bb[i]=0.;
2350  }
2351 
2352  }else{
2353 
2354  /* prescribed boundary conditions */
2355 
2356  /* next statement makes sure that br is reset to zero at the
2357  start of rhs.f */
2358 
2359  nactmech=neq[1];
2360 
2361  for(i=0;i<neq[1];i++){bi[i]=0.;}
2362 
2363  for(i=0;i<*nboun;i++){
2364  xbounr[i]=xbounact[l**nboun+i];
2365  }
2366 
2367  for(j=0;j<neq[1];j++){fr[j]=0.;ubr[j]=0.;}
2368  for(i=0;i<*nboun;i++){
2369  ic=neq[1]+i;
2370  for(j=jq[ic]-1;j<jq[ic+1]-1;j++){
2371  ir=irow[j]-1;
2372  fr[ir]=fr[ir]-au[j]*xbounr[i];
2373  ubr[ir]=fr[ir];
2374  }
2375  }
2376  if(*isolver==0){
2377 #ifdef SPOOLES
2378  spooles_solve(ubr,&neq[1]);
2379 #endif
2380  }
2381  else if(*isolver==4){
2382 #ifdef SGI
2383  sgi_solve(ubr,token);
2384 #endif
2385  }
2386  else if(*isolver==5){
2387 #ifdef TAUCS
2388  tau_solve(ubr,&neq[1]);
2389 #endif
2390  }
2391  else if(*isolver==7){
2392 #ifdef PARDISO
2393  pardiso_solve(ubr,&neq[1],&symmetryflag);
2394 #endif
2395  }
2396  FORTRAN(op,(&neq[1],ubr,mubr,adb,aub,jq,irow));
2397 
2398  for(i=0;i<neq[1];i++){
2399  br[i]+=freq[l]*(freq[l]*mubr[i]);
2400  bi[i]+=freq[l]*(-alpham*mubr[i]-betam*fr[i]);
2401  }
2402 
2403  /* real and imaginary modal coefficients */
2404 
2405  for(i=0;i<nev;i++){
2406  aa[i]=0.;
2407  for(j=0;j<neq[1];j++){
2408  aa[i]+=z[(long long)i*neq[1]+j]*br[j];
2409  }
2410  }
2411 
2412  for(i=0;i<nev;i++){
2413  bb[i]=0.;
2414  for(j=0;j<neq[1];j++){
2415  bb[i]+=z[(long long)i*neq[1]+j]*bi[j];
2416  }
2417  }
2418  }
2419 
2420  /* calculating the modal coefficients */
2421 
2422  if(dashpot==0){
2423  for(i=0;i<nev;i++){
2424  dd=pow(d[i]-pow(freq[l],2),2)+
2425  pow(fric[i],2)*pow(freq[l],2);
2426  bjr[i]=(aa[i]*(d[i]-freq[l]*freq[l])+
2427  bb[i]*fric[i]*freq[l])/dd;
2428  bji[i]=(bb[i]*(d[i]-freq[l]*freq[l])-
2429  aa[i]*fric[i]*freq[l])/dd;
2430  }
2431  }else{
2432  nev2=2*nev;
2433  NNEW(am,double,nev2*nev2);
2434  NNEW(bm,double,nev2);
2435  NNEW(ipiv,ITG,nev2);
2436 
2437  for(i=0;i<nev2;i++){
2438  for(j=0;j<nev2;j++){
2439  am[i*nev2+j]=0.;
2440  }
2441  bm[i]=0.;
2442  }
2443  for(i=0;i<nev;i++){
2444  am[i*nev2+i]=d[i]-freq[l]*freq[l];
2445  am[(i+nev)*nev2+i]=-fric[i]*freq[l];
2446  bm[i]=aa[i];
2447  am[i*nev2+nev+i]=-am[(i+nev)*nev2+i];
2448  am[(i+nev)*nev2+nev+i]=am[i*nev2+i];
2449  bm[nev+i]=bb[i];
2450  for(j=0;j<nev;j++){
2451  am[(j+nev)*nev2+i]=am[(j+nev)*nev2+i]
2452  -cc[i*nev+j]*freq[l];
2453  am[j*nev2+nev+i]=am[j*nev2+nev+i]
2454  +cc[i*nev+j]*freq[l];
2455  }
2456  }
2457 
2458  /* solving the system of equations */
2459 
2460  FORTRAN(dgesv,(&nev2,&nrhs,am,&nev2,ipiv,bm,&nev2,&info));
2461  if(info!=0){
2462  printf("*ERROR in steadystate: fatal termination of dgesv\n");
2463  printf(" info=%" ITGFORMAT "\n",info);
2464 /* FORTRAN(stop,());*/
2465  }
2466 
2467  /* storing the solution in bjr and bji */
2468 
2469  for(i=0;i<nev;i++){
2470  bjr[i]=bm[i];
2471  bji[i]=bm[nev+i];
2472  }
2473 
2474  SFREE(am);SFREE(bm);SFREE(ipiv);
2475  }
2476 
2477  /* calculating the real response */
2478 
2479  if(iprescribedboundary){
2480  if(nmdnode==0){
2481  memcpy(&br[0],&ubr[0],sizeof(double)*neq[1]);
2482  }else{
2483  for(i=0;i<nmddof;i++){
2484  br[imddof[i]]=ubr[imddof[i]];
2485  }
2486  }
2487  }
2488  else{
2489  if(nmdnode==0){
2490  DMEMSET(br,0,neq[1],0.);
2491  }else{
2492  for(i=0;i<nmddof;i++){
2493  br[imddof[i]]=0.;
2494  }
2495  }
2496  }
2497 
2498  if(!cyclicsymmetry){
2499  if(nmdnode==0){
2500  for(i=0;i<neq[1];i++){
2501  for(j=0;j<nev;j++){
2502  br[i]+=bjr[j]*z[(long long)j*neq[1]+i];
2503  }
2504  }
2505  }else{
2506  for(i=0;i<nmddof;i++){
2507  for(j=0;j<nev;j++){
2508  br[imddof[i]]+=bjr[j]*z[(long long)j*neq[1]+imddof[i]];
2509  }
2510  }
2511  }
2512  }else{
2513  for(i=0;i<nmddof;i++){
2514  FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
2515  if(id!=0){
2516  if(izdof[id-1]==imddof[i]){
2517  for(j=0;j<nev;j++){
2518  br[imddof[i]]+=bjr[j]*z[(long long)j*nzdof+id-1];
2519  }
2520  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2521  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2522  }
2523  }
2524 
2525  /* calculating the imaginary response */
2526 
2527  if(nmdnode==0){
2528  DMEMSET(bi,0,neq[1],0.);
2529  }else{
2530  for(i=0;i<nmddof;i++){
2531  bi[imddof[i]]=0.;
2532  }
2533  }
2534 
2535  if(!cyclicsymmetry){
2536  if(nmdnode==0){
2537  for(i=0;i<neq[1];i++){
2538  for(j=0;j<nev;j++){
2539  bi[i]+=bji[j]*z[(long long)j*neq[1]+i];
2540  }
2541  }
2542  }else{
2543  for(i=0;i<nmddof;i++){
2544  for(j=0;j<nev;j++){
2545  bi[imddof[i]]+=bji[j]*z[(long long)j*neq[1]+imddof[i]];
2546  }
2547  }
2548  }
2549  }else{
2550  for(i=0;i<nmddof;i++){
2551  FORTRAN(nident,(izdof,&imddof[i],&nzdof,&id));
2552  if(id!=0){
2553  if(izdof[id-1]==imddof[i]){
2554  for(j=0;j<nev;j++){
2555  bi[imddof[i]]+=bji[j]*z[(long long)j*nzdof+id-1];
2556  }
2557  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2558  }else{printf("*ERROR in steadystate\n");FORTRAN(stop,());}
2559  }
2560  }
2561 
2562  if(nmdnode==0){
2563 
2564  /* magnitude and phase of the response */
2565 
2566  for(i=0;i<neq[1];i++){
2567  breal=br[i];
2568  br[i]=sqrt(br[i]*br[i]+bi[i]*bi[i]);
2569  if(fabs(breal)<1.e-10){
2570  if(bi[i]>0.){bi[i]=pi/2.;}
2571  else{bi[i]=-pi/2.;}
2572  }
2573  else{
2574  bi[i]=atan(bi[i]/breal);
2575  if(breal<0.){bi[i]+=pi;}
2576  }
2577  }
2578 
2579  /* correction for the sinus terms */
2580 
2581  if((l!=0)&&(2*(ITG)floor(l/2.+0.1)==l)){
2582  for(i=0;i<neq[1];i++){
2583 // bi[i]-=pi/2.;}
2584  bi[i]+=pi/2.;}
2585  }
2586 
2587  /* contribution to the time response */
2588 
2589  for(j=0;j<nfour;j++){
2590  time=tmin+2.*pi/freqnh[ii]*(double)j/(double)nfour;
2591  for(i=0;i<neq[1];i++){
2592  btot[j*neq[1]+i]+=br[i]*cos(freq[l]*time+bi[i]);
2593  bp[j*neq[1]+i]-=freq[l]*br[i]*sin(freq[l]*time+bi[i]);
2594  }
2595  }
2596  }else{
2597 
2598  /* magnitude and phase of the response */
2599 
2600  for(jj=0;jj<nmddof;jj++){
2601  i=imddof[jj];
2602  breal=br[i];
2603  br[i]=sqrt(br[i]*br[i]+bi[i]*bi[i]);
2604  if(fabs(breal)<1.e-10){
2605  if(bi[i]>0.){bi[i]=pi/2.;}
2606  else{bi[i]=-pi/2.;}
2607  }
2608  else{
2609  bi[i]=atan(bi[i]/breal);
2610  if(breal<0.){bi[i]+=pi;}
2611  }
2612  }
2613 
2614  /* correction for the sinus terms */
2615 
2616  if((l!=0)&&(2*(ITG)floor(l/2.+0.1)==l)){
2617  for(jj=0;jj<nmddof;jj++){
2618  i=imddof[jj];
2619 // bi[i]-=pi/2.;}
2620  bi[i]+=pi/2.;}
2621  }
2622 
2623  /* contribution to the time response */
2624 
2625  for(j=0;j<nfour;j++){
2626  time=tmin+2.*pi/freqnh[ii]*(double)j/(double)nfour;
2627  for(jj=0;jj<nmddof;jj++){
2628  i=imddof[jj];
2629  btot[j*neq[1]+i]+=br[i]*cos(freq[l]*time+bi[i]);
2630  bp[j*neq[1]+i]-=freq[l]*br[i]*sin(freq[l]*time+bi[i]);
2631  }
2632  }
2633  }
2634 
2635  /* resetting the part of br occupied by the variables to be printed
2636  to zero */
2637 
2638  if(!iprescribedboundary){
2639  if(nmdnode==0){
2640  DMEMSET(br,0,neq[1],0.);
2641  }else{
2642  for(i=0;i<nmddof;i++){
2643  br[imddof[i]]=0.;
2644  }
2645  }
2646  }
2647 
2648  }
2649 
2650  SFREE(xforcr);SFREE(xloadr);SFREE(xbodyr);SFREE(br);SFREE(bi);SFREE(freq);
2651  SFREE(bjr);SFREE(bji);SFREE(aa);SFREE(bb);
2652 
2653  if(*nbody>0) SFREE(ipobody);
2654  if(iprescribedboundary) {SFREE(xbounr);SFREE(fr);SFREE(ubr);SFREE(mubr);}
2655 
2656 
2657  /* result fields */
2658 
2659  NNEW(vr,double,mt**nk);
2660 
2661  if(intpointvar==1){
2662  NNEW(fn,double,mt**nk);
2663  NNEW(stn,double,6**nk);
2664  NNEW(stx,double,6*mi[0]**ne);
2665 
2666  if(*ithermal>1) {
2667  NNEW(qfn,double,3**nk);
2668  NNEW(qfx,double,3*mi[0]**ne);}
2669 
2670  if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
2671  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
2672  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
2673 
2674  NNEW(eei,double,6*mi[0]**ne);
2675  if(*nener==1){
2676  NNEW(stiini,double,6*mi[0]**ne);
2677  NNEW(enerini,double,mi[0]**ne);}
2678  }
2679 
2680  /* storing the results */
2681 
2682  for(l=0;l<nfour;l++){
2683  time=tmin+2.*pi/freqnh[ii]*(double)l/(double)nfour;
2684  ptime=time;
2685 
2686  if(nmdnode==0){
2687  DMEMSET(vr,0,mt**nk,0.);
2688  }else{
2689  for(jj=0;jj<nmdnode;jj++){
2690  i=imdnode[jj]-1;
2691  for(j=1;j<4;j++){
2692  vr[mt*i+j]=0.;
2693  }
2694  }
2695  }
2696 
2697  /* calculating displacements/temperatures */
2698 
2699  *nmethod=4;
2700  FORTRAN(dynresults,(nk,vr,ithermal,nactdof,vold,nodeboun,
2701  ndirboun,&xbounacttime[l**nboun],nboun,ipompc,nodempc,
2702  coefmpc,labmpc,nmpc,&btot[l*neq[1]],&bp[l*neq[1]],veold,&dtime,mi,
2703  imdnode,&nmdnode,imdboun,&nmdboun,imdmpc,&nmdmpc,nmethod,&time));
2704  *nmethod=5;
2705 
2706  results(co,nk,kon,ipkon,lakon,ne,vr,stn,inum,
2707  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
2708  ielmat,ielorien,norien,orab,ntmat_,t0,t1,
2709  ithermal,prestr,iprestr,filab,eme,emn,een,
2710  iperturb,f,fn,nactdof,&iout,qa,
2711  vold,&btot[l*neq[1]],nodeboun,ndirboun,
2712  &xbounacttime[l**nboun],nboun,
2713  ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
2714  veold,accold,&bet,&gam,&dtime,&time,&xnull,
2715  plicon,nplicon,plkcon,nplkcon,
2716  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
2717  &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,
2718  enern,emeini,xstaten,eei,enerini,cocon,ncocon,
2719  set,nset,istartset,iendset,ialset,nprint,prlab,prset,
2720  qfx,qfn,trab,inotr,ntrans,fmpc,nelemload,nload,ikmpc,
2721  ilmpc,istep,&iinc,springarea,&reltime,&ne0,xforc,nforc,
2722  thicke,shcon,nshcon,
2723  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
2724  &mortar,islavact,cdn,islavnode,nslavnode,&ntie,clearini,
2725  islavsurf,ielprop,prop);
2726 
2727  (*kode)++;
2728  mode=-1;
2729 
2730  if(strcmp1(&filab[1044],"ZZS")==0){
2731  NNEW(neigh,ITG,40**ne);
2732  NNEW(ipneigh,ITG,*nk);
2733  }
2734 
2735  frd(co,&nkg,kon,ipkon,lakon,&neg,vr,stn,inum,nmethod,
2736  kode,filab,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
2737  nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
2738  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2739  mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&neg,
2740  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
2741  thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
2742 
2743  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2744 
2745  }
2746 
2747  SFREE(vr);SFREE(btot);SFREE(bp);
2748 
2749  if(intpointvar==1){
2750  SFREE(fn);SFREE(stn);SFREE(stx);SFREE(eei);
2751  if(*ithermal>1) {SFREE(qfn);SFREE(qfx);}
2752 
2753  if(strcmp1(&filab[261],"E ")==0) SFREE(een);
2754  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2755  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
2756 
2757  if(*nener==1){SFREE(stiini);SFREE(enerini);}
2758  }
2759 
2760  }
2761  SFREE(xforcact);SFREE(xloadact);SFREE(xbodyact);SFREE(xbounact);
2762  SFREE(xbounacttime);SFREE(freqnh);
2763  if(*ithermal==1) SFREE(t1act);
2764  if(iprescribedboundary){
2765  if(*isolver==0){
2766 #ifdef SPOOLES
2767  spooles_cleanup();
2768 #endif
2769  }
2770  else if(*isolver==4){
2771 #ifdef SGI
2772  sgi_cleanup(token);
2773 #endif
2774  }
2775  else if(*isolver==5){
2776 #ifdef TAUCS
2777  tau_cleanup();
2778 #endif
2779  }
2780  else if(*isolver==7){
2781 #ifdef PARDISO
2782  pardiso_cleanup(&neq[1],&symmetryflag);
2783 #endif
2784  }
2785  }
2786 
2787  SFREE(ikactmech);
2788 
2789  }
2790 
2791  SFREE(adb);SFREE(aub);SFREE(z);SFREE(d);SFREE(inum);
2792 
2793  if(!cyclicsymmetry){
2794  SFREE(ad);SFREE(au);
2795  }else{
2796  SFREE(izdof);SFREE(nm);
2797 
2798  *nk/=nsectors;
2799  *ne/=nsectors;
2800  *nboun/=nsectors;
2801  neq[1]=neq[1]*2/nsectors;
2802 
2803  RENEW(co,double,3**nk);
2804  if(*ithermal!=0){
2805  RENEW(t0,double,*nk);
2806  RENEW(t1old,double,*nk);
2807  RENEW(t1,double,*nk);
2808  if(*nam>0) RENEW(iamt1,ITG,*nk);
2809  }
2810  RENEW(nactdof,ITG,mt**nk);
2811  if(*ntrans>0) RENEW(inotr,ITG,2**nk);
2812  RENEW(kon,ITG,*nkon);
2813  RENEW(ipkon,ITG,*ne);
2814  RENEW(lakon,char,8**ne);
2815  RENEW(ielmat,ITG,mi[2]**ne);
2816  if(*norien>0) RENEW(ielorien,ITG,mi[2]**ne);
2817  RENEW(nodeboun,ITG,*nboun);
2818  RENEW(ndirboun,ITG,*nboun);
2819  if(*nam>0) RENEW(iamboun,ITG,*nboun);
2820  RENEW(xboun,double,*nboun);
2821  RENEW(xbounold,double,*nboun);
2822  RENEW(ikboun,ITG,*nboun);
2823  RENEW(ilboun,ITG,*nboun);
2824 
2825  /* recovering the original multiple point constraints */
2826 
2827  RENEW(ipompc,ITG,*nmpc);
2828  RENEW(nodempc,ITG,3**mpcend);
2829  RENEW(coefmpc,double,*mpcend);
2830  RENEW(labmpc,char,20**nmpc+1);
2831  RENEW(ikmpc,ITG,*nmpc);
2832  RENEW(ilmpc,ITG,*nmpc);
2833  RENEW(fmpc,double,*nmpc);
2834 
2835  *nmpc=nmpcold;
2836  *mpcend=mpcendold;
2837  for(i=0;i<*nmpc;i++){ipompc[i]=ipompcold[i];}
2838  for(i=0;i<3**mpcend;i++){nodempc[i]=nodempcold[i];}
2839  for(i=0;i<*mpcend;i++){coefmpc[i]=coefmpcold[i];}
2840  for(i=0;i<20**nmpc;i++){labmpc[i]=labmpcold[i];}
2841  for(i=0;i<*nmpc;i++){ikmpc[i]=ikmpcold[i];}
2842  for(i=0;i<*nmpc;i++){ilmpc[i]=ilmpcold[i];}
2843  SFREE(ipompcold);SFREE(nodempcold);SFREE(coefmpcold);
2844  SFREE(labmpcold);SFREE(ikmpcold);SFREE(ilmpcold);
2845 
2846  RENEW(vold,double,mt**nk);
2847  RENEW(veold,double,mt**nk);
2848  RENEW(eme,double,6*mi[0]**ne);
2849  if(*nener==1)RENEW(ener,double,mi[0]**ne);
2850 
2851 /* distributed loads */
2852 
2853  for(i=0;i<*nload;i++){
2854  if(nelemload[2*i]<=*ne*nsectors){
2855  nelemload[2*i]-=*ne*nelemload[2*i+1];
2856  }else{
2857  nelemload[2*i]-=*ne*(nsectors+nelemload[2*i+1]-1);
2858  }
2859  }
2860 
2861  /* sorting the elements with distributed loads */
2862 
2863  if(*nload>0){
2864  if(*nam>0){
2865  FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
2866  }else{
2867  FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
2868  }
2869  }
2870 
2871 /* point loads */
2872 
2873  for(i=0;i<*nforc;i++){
2874  if(nodeforc[2*i]<=*nk*nsectors){
2875  nodeforc[2*i]-=*nk*nodeforc[2*i+1];
2876  }else{
2877  nodeforc[2*i]-=*nk*(nsectors+nodeforc[2*i+1]-1);
2878  }
2879  }
2880  }
2881 
2882  SFREE(xstiff);SFREE(fric);
2883 
2884  if(dashpot){SFREE(cc);}
2885 
2886  if(nherm!=1){SFREE(xmr);SFREE(xmi);}
2887 
2888  SFREE(imddof);SFREE(imdnode);SFREE(imdboun);SFREE(imdmpc);SFREE(imdelem);
2889 
2890  *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
2891  *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
2892  *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;*veoldp=veold;
2893  *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
2894  *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
2895  *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
2896  *fmpcp=fmpc;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
2897 
2898 // (*ttime)+=(*tper);
2899 
2900  return;
2901 }
#define ITGFORMAT
Definition: CalculiX.h:52
void spooles_solve(double *b, ITG *neq)
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 drfftf(N, R, WSAVE, isave)
Definition: drfftf.f:466
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))
void sgi_solve(double *b, ITG token)
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
void pardiso_cleanup(ITG *neq, ITG *symmetryflag)
#define DMEMSET(a, b, c, d)
Definition: CalculiX.h:45
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 dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:57
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 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 writepf(d, bjr, bji, freq, nev, mode, nherm)
Definition: writepf.f:19
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 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 nident(x, px, n, id)
Definition: nident.f:25
subroutine isortiiddc(ix1, ix2, dy1, dy2, cy, n, kflag)
Definition: isortiiddc.f:5
subroutine drffti(N, WSAVE, isave)
Definition: drfftf.f:533
void tau_solve(double *b, ITG *neq)
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()
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:5
#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 pardiso_solve(double *b, ITG *neq, ITG *symmetryflag)
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 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)