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

Go to the source code of this file.

Functions

void complexfreq (double **cop, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp, ITG *ne, ITG **nodebounp, ITG **ndirbounp, double **xbounp, ITG *nboun, ITG **ipompcp, ITG **nodempcp, double **coefmpcp, char **labmpcp, ITG *nmpc, ITG *nodeforc, ITG *ndirforc, double *xforc, ITG *nforc, ITG *nelemload, char *sideload, double *xload, ITG *nload, ITG **nactdofp, ITG *neq, ITG *nzl, ITG *icol, ITG *irow, ITG *nmethod, ITG **ikmpcp, ITG **ilmpcp, ITG **ikbounp, ITG **ilbounp, double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon, double *cocon, ITG *ncocon, double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp, ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_, double **t0p, double **t1p, ITG *ithermal, double *prestr, ITG *iprestr, double **voldp, ITG *iperturb, double **stip, ITG *nzs, double *timepar, double *xmodal, double **veoldp, char *amname, double *amta, ITG *namta, ITG *nam, ITG *iamforc, ITG *iamload, ITG **iamt1p, ITG *jout, ITG *kode, char *filab, double **emep, double *xforcold, double *xloadold, double **t1oldp, ITG **iambounp, double **xbounoldp, ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon, ITG *nplkcon, double *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 **ialsetp, ITG *nprint, char *prlab, char *prset, ITG *nener, double *trab, ITG **inotrp, ITG *ntrans, double **fmpcp, char *cbody, ITG *ibody, double *xbody, ITG *nbody, double *xbodyold, ITG *istep, ITG *isolver, ITG *jq, char *output, ITG *mcs, ITG *nkon, ITG *mpcend, ITG *ics, double *cs, ITG *ntie, char *tieset, ITG *idrct, ITG *jmax, double *ctrl, ITG *itpamp, double *tietol, ITG *nalset, ITG *ikforc, ITG *ilforc, double *thicke, char *jobnamef, ITG *mei, ITG *nmat, ITG *ielprop, double *prop)
 

Function Documentation

void complexfreq ( double **  cop,
ITG nk,
ITG **  konp,
ITG **  ipkonp,
char **  lakonp,
ITG ne,
ITG **  nodebounp,
ITG **  ndirbounp,
double **  xbounp,
ITG nboun,
ITG **  ipompcp,
ITG **  nodempcp,
double **  coefmpcp,
char **  labmpcp,
ITG nmpc,
ITG nodeforc,
ITG ndirforc,
double *  xforc,
ITG nforc,
ITG nelemload,
char *  sideload,
double *  xload,
ITG nload,
ITG **  nactdofp,
ITG neq,
ITG nzl,
ITG icol,
ITG irow,
ITG nmethod,
ITG **  ikmpcp,
ITG **  ilmpcp,
ITG **  ikbounp,
ITG **  ilbounp,
double *  elcon,
ITG nelcon,
double *  rhcon,
ITG nrhcon,
double *  cocon,
ITG ncocon,
double *  alcon,
ITG nalcon,
double *  alzero,
ITG **  ielmatp,
ITG **  ielorienp,
ITG norien,
double *  orab,
ITG ntmat_,
double **  t0p,
double **  t1p,
ITG ithermal,
double *  prestr,
ITG iprestr,
double **  voldp,
ITG iperturb,
double **  stip,
ITG nzs,
double *  timepar,
double *  xmodal,
double **  veoldp,
char *  amname,
double *  amta,
ITG namta,
ITG nam,
ITG iamforc,
ITG iamload,
ITG **  iamt1p,
ITG jout,
ITG kode,
char *  filab,
double **  emep,
double *  xforcold,
double *  xloadold,
double **  t1oldp,
ITG **  iambounp,
double **  xbounoldp,
ITG iexpl,
double *  plicon,
ITG nplicon,
double *  plkcon,
ITG nplkcon,
double *  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 **  ialsetp,
ITG nprint,
char *  prlab,
char *  prset,
ITG nener,
double *  trab,
ITG **  inotrp,
ITG ntrans,
double **  fmpcp,
char *  cbody,
ITG ibody,
double *  xbody,
ITG nbody,
double *  xbodyold,
ITG istep,
ITG isolver,
ITG jq,
char *  output,
ITG mcs,
ITG nkon,
ITG mpcend,
ITG ics,
double *  cs,
ITG ntie,
char *  tieset,
ITG idrct,
ITG jmax,
double *  ctrl,
ITG itpamp,
double *  tietol,
ITG nalset,
ITG ikforc,
ITG ilforc,
double *  thicke,
char *  jobnamef,
ITG mei,
ITG nmat,
ITG ielprop,
double *  prop 
)

Definition at line 36 of file complexfreq.c.

72  {
73 
74  char fneig[132]="",description[13]=" ",*lakon=NULL,*labmpc=NULL,
75  *lakont=NULL;
76 
77  ITG nev,i,j,k,idof,*inum=NULL,*ipobody=NULL,inewton=0,id,
78  iinc=0,l,iout=1,ielas,icmd,ifreebody,mode,m,nherm,
79  *kon=NULL,*ipkon=NULL,*ielmat=NULL,*ielorien=NULL,*islavact=NULL,
80  *inotr=NULL,*nodeboun=NULL,*ndirboun=NULL,*iamboun=NULL,*ikboun=NULL,
81  *ilboun=NULL,*nactdof=NULL,*ipompc=NULL,*nodempc=NULL,*ikmpc=NULL,
82  *ilmpc=NULL,nsectors,nmd,nevd,*nm=NULL,*iamt1=NULL,*islavnode=NULL,
83  ngraph=1,nkg,neg,ne0,ij,lprev,nope,indexe,ilength,*nslavnode=NULL,
84  *ipneigh=NULL,*neigh=NULL,index,im,cyclicsymmetry,inode,
85  *ialset=*ialsetp,mt=mi[1]+1,kmin,kmax,i1,
86  *iter=NULL,lint,lfin,kk,kkv,kk6,kkx,icomplex,igeneralizedforce,
87  idir,*inumt=NULL,icntrl,imag,jj,is,l1,*inocs=NULL,ml1,l2,nkt,net,
88  *ipkont=NULL,*ielmatt=NULL,*inotrt=NULL,*kont=NULL,node,iel,*ielcs=NULL,
89  ielset,*istartnmd=NULL,*iendnmd=NULL,inmd,neqact,*nshcon=NULL,
90  *ipev=NULL,icfd=0,*inomat=NULL,mortar=0,*islavsurf=NULL;
91 
92  long long i2;
93 
94  double *d=NULL, *z=NULL,*stiini=NULL,*cc=NULL,*v=NULL,*zz=NULL,*emn=NULL,
95  *stn=NULL, *stx=NULL, *een=NULL, *adb=NULL,*xstiff=NULL,*cdn=NULL,
96  *aub=NULL,*f=NULL, *fn=NULL,*epn=NULL,*xstateini=NULL,
97  *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*qfn=NULL,
98  *qfx=NULL, *cgr=NULL, *au=NULL,dtime,reltime,*t0=NULL,*t1=NULL,*t1old=NULL,
99  sum,qa[3],cam[5],accold[1],bet,gam,*ad=NULL,alpham,betam,
100  *co=NULL,*xboun=NULL,*xbounold=NULL,*vold=NULL,*emeini=NULL,
101  *eme=NULL,*ener=NULL,*coefmpc=NULL,*fmpc=NULL,*veold=NULL,
102  *adc=NULL,*auc=NULL,*zc=NULL,*fnr=NULL,*fni=NULL,setnull,deltmx,dd,
103  theta,*vini=NULL,*vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,*vmax=NULL,
104  *stnmax=NULL,*cstr=NULL,*sti=*stip,time0=0.0,time=0.0,zero=0.0,
105  *springarea=NULL,*eenmax=NULL,*aa=NULL,*bb=NULL,*xx=NULL,
106  *eiga=NULL,*eigb=NULL,*eigxx=NULL,*temp=NULL,*coefmpcnew=NULL,xreal,
107  ximag,t[3],*vt=NULL,*t1t=NULL,*stnt=NULL,*eent=NULL,*fnt=NULL,*enernt=NULL,
108  *stxt=NULL,pi,ctl,stl,*cot=NULL,*qfnt=NULL,vreal,vimag,constant,stnreal,
109  stnimag,freq,*emnt=NULL,*shcon=NULL,*eig=NULL,*clearini=NULL,
110  *eigxr=NULL,*eigxi=NULL,*xmac=NULL,*bett=NULL,*betm=NULL,*xmaccpx=NULL,
111  fmin=0.,fmax=1.e30,*xmr=NULL,*xmi=NULL,*zi=NULL,*eigx=NULL,
112  *pslavsurf=NULL,*pmastsurf=NULL,*cdnr=NULL,*cdni=NULL,*tinc,*tper,*tmin,*tmax;
113 
114  FILE *f1;
115 
116 #ifdef SGI
117  ITG token;
118 #endif
119 
120  pi=4.*atan(1.);
121  constant=180./pi;
122 
123  co=*cop;kon=*konp;ipkon=*ipkonp;lakon=*lakonp;ielmat=*ielmatp;
124  ielorien=*ielorienp;inotr=*inotrp;nodeboun=*nodebounp;
125  ndirboun=*ndirbounp;iamboun=*iambounp;xboun=*xbounp;
126  xbounold=*xbounoldp;ikboun=*ikbounp;ilboun=*ilbounp;nactdof=*nactdofp;
127  vold=*voldp;eme=*emep;ener=*enerp;ipompc=*ipompcp;nodempc=*nodempcp;
128  coefmpc=*coefmpcp;labmpc=*labmpcp;ikmpc=*ikmpcp;ilmpc=*ilmpcp;
129  fmpc=*fmpcp;veold=*veoldp;iamt1=*iamt1p;t0=*t0p;t1=*t1p;t1old=*t1oldp;
130 
131  tinc=&timepar[0];
132  tper=&timepar[1];
133  tmin=&timepar[2];
134  tmax=&timepar[3];
135 
136  if(ithermal[0]<=1){
137  kmin=1;kmax=3;
138  }else if(ithermal[0]==2){
139  kmin=0;kmax=mi[1];if(kmax>2)kmax=2;
140  }else{
141  kmin=0;kmax=3;
142  }
143 
144  NNEW(xstiff,double,(long long)27*mi[0]**ne);
145 
146  dtime=*tinc;
147 
148  alpham=xmodal[0];
149  betam=xmodal[1];
150 
151  dd=ctrl[16];deltmx=ctrl[26];
152 
153  /* determining nzl */
154 
155  *nzl=0;
156  for(i=neq[1];i>0;i--){
157  if(icol[i-1]>0){
158  *nzl=i;
159  break;
160  }
161  }
162 
163  /* opening the eigenvalue file and checking for cyclic symmetry */
164 
165  strcpy(fneig,jobnamec);
166  strcat(fneig,".eig");
167 
168  if((f1=fopen(fneig,"rb"))==NULL){
169  printf("*ERROR in complexfreq: cannot open eigenvalue file for reading");
170  exit(0);
171  }
172 
173  printf(" *INFO in complexfreq: if there are problems reading the .eig file this may be due to:\n");
174  printf(" 1) the nonexistence of the .eig file\n");
175  printf(" 2) other boundary conditions than in the input deck\n");
176  printf(" which created the .eig file\n\n");
177 
178  if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
179  printf("*ERROR in complexfreq reading the cyclic symmetry flag in the eigenvalue file");
180  exit(0);
181  }
182 
183  if(fread(&nherm,sizeof(ITG),1,f1)!=1){
184  printf("*ERROR in complexfreq reading the Hermitian flag in the eigenvalue file");
185  exit(0);
186  }
187 
188  if(nherm!=1){
189  printf("*ERROR in complexfreq: the eigenvectors in the .eig-file result\n");
190  printf(" from a non-Hermitian eigenvalue problem. The complex\n");
191  printf(" frequency procedure cannot handle that yet\n\n");
192  FORTRAN(stop,());
193  }
194 
195  nsectors=1;
196 
197  if(!cyclicsymmetry){
198 
199  nkg=*nk;
200  neg=*ne;
201 
202  if(fread(&nev,sizeof(ITG),1,f1)!=1){
203  printf("*ERROR in complexfreq reading the number of eigenvalues in the eigenvalue file...");
204  exit(0);
205  }
206 
207 
208  if(nherm==1){
209  NNEW(d,double,nev);
210  if(fread(d,sizeof(double),nev,f1)!=nev){
211  printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
212  exit(0);
213  }
214  }else{
215  NNEW(d,double,2*nev);
216  if(fread(d,sizeof(double),2*nev,f1)!=2*nev){
217  printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
218  exit(0);
219  }
220  }
221 
222  NNEW(ad,double,neq[1]);
223  NNEW(adb,double,neq[1]);
224  NNEW(au,double,nzs[2]);
225  NNEW(aub,double,nzs[1]);
226 
227  /* reading the stiffness matrix */
228 
229  if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
230  printf("*ERROR in complexfreq reading the diagonal of the stiffness matrix in the eigenvalue file...");
231  exit(0);
232  }
233 
234  if(fread(au,sizeof(double),nzs[2],f1)!=nzs[2]){
235  printf("*ERROR in complexfreq reading the off-diagonal terms of the stiffness matrix in the eigenvalue file...");
236  exit(0);
237  }
238 
239  /* reading the mass matrix */
240 
241  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
242  printf("*ERROR in complexfreq reading the diagonal of the mass matrix in eigenvalue file...");
243  exit(0);
244  }
245 
246  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
247  printf("*ERROR in complexfreq reading the off-diagonals of the mass matrix in the eigenvalue file...");
248  exit(0);
249  }
250 
251  /* reading the eigenvectors */
252 
253  if(nherm==1){
254  NNEW(z,double,neq[1]*nev);
255  if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
256  printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
257  exit(0);
258  }
259  }else{
260  NNEW(z,double,2*neq[1]*nev);
261  if(fread(z,sizeof(double),2*neq[1]*nev,f1)!=2*neq[1]*nev){
262  printf("*ERROR in complexfreq reading the eigenvectors in the eigenvalue file...");
263  exit(0);
264  }
265  }
266 
267  NNEW(nm,ITG,nev);
268  for(i=0;i<nev;i++){nm[i]=-1;}
269  }
270  else{
271 
272  if(*nmethod==6){
273  printf("*ERROR in complexfreq: Coriolis forces cannot\n");
274  printf(" be combined with cyclic symmetry\n\n");
275  FORTRAN(stop,());
276  }
277 
278  nev=0;
279  do{
280  if(fread(&nmd,sizeof(ITG),1,f1)!=1){
281  break;
282  }
283  if(fread(&nevd,sizeof(ITG),1,f1)!=1){
284  printf("*ERROR in complexfreq reading the number of eigenvalues in the eigenvalue file...");
285  exit(0);
286  }
287  if(nev==0){
288  if(nherm==1){NNEW(d,double,nevd);
289  }else{NNEW(d,double,2*nevd);}
290  NNEW(nm,ITG,nevd);
291  }else{
292  printf("*ERROR in complexfreq: flutter forces cannot\n");
293  printf(" be combined with multiple modal diameters\n");
294  printf(" in cyclic symmetry calculations\n\n");
295  FORTRAN(stop,());
296  }
297 
298  if(nherm==1){
299  if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
300  printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
301  exit(0);
302  }
303  }else{
304  if(fread(&d[nev],sizeof(double),2*nevd,f1)!=2*nevd){
305  printf("*ERROR in complexfreq reading the eigenvalues in the eigenvalue file...");
306  exit(0);
307  }
308  }
309  for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
310 
311  if(nev==0){
312  NNEW(adb,double,neq[1]);
313  NNEW(aub,double,nzs[1]);
314 
315  if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
316  printf("*ERROR in complexfreq reading the diagonal of the mass matrix in the eigenvalue file...");
317  exit(0);
318  }
319 
320  if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
321  printf("*ERROR in complexfreq reading the off-diagonals of the mass matrix in the eigenvalue file...");
322  exit(0);
323  }
324  }
325 
326  if(nev==0){
327  NNEW(z,double,neq[1]*nevd);
328  }else{
329  RENEW(z,double,(long long)neq[1]*(nev+nevd));
330  }
331 
332  if(fread(&z[neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
333  printf("*ERROR in complexfreq reading eigenvectors in the eigenvalue file...");
334  exit(0);
335  }
336  nev+=nevd;
337  }while(1);
338 
339  /* removing double eigenmodes */
340 
341  j=-1;
342  for(i=0;i<nev;i++){
343  if((i/2)*2==i){
344  j++;
345  if(nherm==1){
346  d[j]=d[i];
347  }else{
348  d[2*j]=d[2*i];d[2*j+1]=d[2*i+1];
349  }
350  nm[j]=nm[i];
351  for(k=0;k<neq[1];k++){
352  z[j*neq[1]+k]=z[i*neq[1]+k];
353  }
354  }
355  }
356  nev=j+1;
357  if(nherm==1){RENEW(d,double,nev);}else{RENEW(d,double,2*nev);}
358  RENEW(nm,ITG,nev);
359  RENEW(z,double,neq[1]*nev);
360 
361  /* determining the maximum amount of segments */
362 
363  for(i=0;i<*mcs;i++){
364  if(cs[17*i]>nsectors) nsectors=(ITG)(cs[17*i]+0.5);
365  }
366 
367  /* determining the maximum number of sectors to be plotted */
368 
369  for(j=0;j<*mcs;j++){
370  if(cs[17*j+4]>ngraph) ngraph=(ITG)cs[17*j+4];
371  }
372  nkg=*nk*ngraph;
373  neg=*ne*ngraph;
374 
375  }
376 
377  fclose(f1);
378 
379  /* assigning nodes and elements to sectors */
380 
381  if(cyclicsymmetry){
382  NNEW(inocs,ITG,*nk);
383  NNEW(ielcs,ITG,*ne);
384  ielset=cs[12];
385  if((*mcs!=1)||(ielset!=0)){
386  for(i=0;i<*nk;i++) inocs[i]=-1;
387  for(i=0;i<*ne;i++) ielcs[i]=-1;
388  }
389 
390  for(i=0;i<*mcs;i++){
391  is=cs[17*i+4];
392  if((is==1)&&(*mcs==1)) continue;
393  ielset=cs[17*i+12];
394  if(ielset==0) continue;
395  for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
396  if(ialset[i1]>0){
397  iel=ialset[i1]-1;
398  if(ipkon[iel]<0) continue;
399  ielcs[iel]=i;
400  indexe=ipkon[iel];
401  if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
402  else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
403  else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
404  else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
405  else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
406  else {nope=6;}
407  for(i2=0;i2<nope;++i2){
408  node=kon[indexe+i2]-1;
409  inocs[node]=i;
410  }
411  }
412  else{
413  iel=ialset[i1-2]-1;
414  do{
415  iel=iel-ialset[i1];
416  if(iel>=ialset[i1-1]-1) break;
417  if(ipkon[iel]<0) continue;
418  ielcs[iel]=i;
419  indexe=ipkon[iel];
420  if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
421  else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
422  else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
423  else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
424  else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
425  else {nope=6;}
426  for(i2=0;i2<nope;++i2){
427  node=kon[indexe+i2]-1;
428  inocs[node]=i;
429  }
430  }while(1);
431  }
432  }
433  }
434  }
435 
436  /* check for rigid body modes
437  if there is a jump of 1.e4 in two subsequent eigenvalues
438  all eigenvalues preceding the jump are considered to
439  be rigid body modes and their frequency is set to zero */
440 
441  if(nherm==1){
442  setnull=1.;
443  for(i=nev-2;i>-1;i--){
444  if(fabs(d[i])<0.0001*fabs(d[i+1])) setnull=0.;
445  d[i]*=setnull;
446  }
447  }else{
448  setnull=1.;
449  for(i=nev-2;i>-1;i--){
450  if(sqrt(d[2*i]*d[2*i]+d[2*i+1]*d[2*i+1])<
451  0.0001*sqrt(d[2*i+2]*d[2*i+2]+d[2*i+3]*d[2*i+3])) setnull=0.;
452  d[2*i]*=setnull;
453  d[2*i+1]*=setnull;
454  }
455  }
456 
457  /* determining the frequency ranges corresponding to one
458  and the same nodal diameter */
459 
460  if(cyclicsymmetry){
461  NNEW(istartnmd,ITG,nev);
462  NNEW(iendnmd,ITG,nev);
463  nmd=0;
464  inmd=nm[0];
465  istartnmd[0]=1;
466  for(i=1;i<nev;i++){
467  if(nm[i]==inmd) continue;
468  iendnmd[nmd]=i;
469  nmd++;
470  istartnmd[nmd]=i+1;
471  inmd=nm[i];
472  }
473  iendnmd[nmd]=nev;
474  nmd++;
475  RENEW(istartnmd,ITG,nmd);
476  RENEW(iendnmd,ITG,nmd);
477  }
478 
479  if(*nmethod==6){
480 
481  /* Coriolis */
482 
483  neqact=neq[1];
484 
485  /* assigning the body forces to the elements */
486 
487  ifreebody=*ne+1;
488  NNEW(ipobody,ITG,2**ne);
489  for(k=1;k<=*nbody;k++){
490  FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
491  iendset,ialset,&inewton,nset,&ifreebody,&k));
492  RENEW(ipobody,ITG,2*(*ne+ifreebody));
493  }
494  RENEW(ipobody,ITG,2*(ifreebody-1));
495 
496  if(cyclicsymmetry){
497  printf("*ERROR in complexfreq: dashpots are not allowed in combination with cyclic symmetry\n");
498  FORTRAN(stop,());
499  }
500 
501  NNEW(adc,double,neq[1]);
502  NNEW(auc,double,nzs[1]);
503  FORTRAN(mafillcorio,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
504  xboun,nboun,
505  ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
506  nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
507  adc,auc,nactdof,icol,jq,irow,neq,nzl,nmethod,
508  ikmpc,ilmpc,ikboun,ilboun,
509  elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
510  ielorien,norien,orab,ntmat_,
511  t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
512  nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
513  xstiff,npmat_,&dtime,matname,mi,ncmat_,
514  ttime,&time0,istep,&iinc,ibody));
515 
516  /* zc = damping matrix * eigenmodes */
517 
518  NNEW(zc,double,neq[1]*nev);
519  for(i=0;i<nev;i++){
520  FORTRAN(op_corio,(&neq[1],&z[i*neq[1]],&zc[i*neq[1]],adc,auc,
521  jq,irow));
522  }
523  SFREE(adc);SFREE(auc);
524 
525  /* cc is the reduced damping matrix (damping matrix mapped onto
526  space spanned by eigenmodes) */
527 
528  NNEW(cc,double,nev*nev);
529  for(i=0;i<nev;i++){
530  for(j=0;j<=i;j++){
531  for(k=0;k<neq[1];k++){
532  cc[i*nev+j]+=z[j*neq[1]+k]*zc[i*neq[1]+k];
533  }
534  }
535  }
536 
537  /* symmetric part of cc matrix */
538 
539  for(i=0;i<nev;i++){
540  for(j=i+1;j<nev;j++){
541  cc[i*nev+j]=-cc[j*nev+i];
542  }
543  }
544  SFREE(zc);
545 
546  /* solving for the complex eigenvalues */
547 
548  NNEW(aa,double,4*nev*nev);
549  NNEW(bb,double,4*nev*nev);
550  NNEW(xx,double,4*nev*nev);
551  NNEW(temp,double,4*nev*nev);
552  NNEW(eiga,double,2*nev);
553  NNEW(eigb,double,2*nev);
554  NNEW(eigxx,double,2*nev);
555  NNEW(iter,ITG,nev);
556 
557  FORTRAN(coriolissolve,(cc,&nev,aa,bb,xx,eiga,eigb,eigxx,
558  iter,d,temp));
559 
560  SFREE(aa);SFREE(bb);SFREE(temp);SFREE(eiga);SFREE(eigb);SFREE(iter);SFREE(cc);
561 
562  }else{
563 
564  /* flutter */
565 
566  /* complex force is being read (e.g. due to fluid flow) */
567 
568  if(cyclicsymmetry){
569  neqact=neq[1]/2;
570  }else{
571  neqact=neq[1];
572  }
573  NNEW(zc,double,2*neqact*nev);
574  NNEW(aa,double,4*nev*nev);
575 
576  FORTRAN(readforce,(zc,&neqact,&nev,nactdof,ikmpc,nmpc,
577  ipompc,nodempc,mi,coefmpc,jobnamef,
578  aa,&igeneralizedforce));
579 
580  NNEW(bb,double,4*nev*nev);
581  NNEW(xx,double,4*nev*nev);
582  NNEW(eiga,double,2*nev);
583  NNEW(eigb,double,2*nev);
584  NNEW(eigxx,double,2*nev);
585  NNEW(iter,ITG,nev);
586  FORTRAN(forcesolve,(zc,&nev,aa,bb,xx,eiga,eigb,eigxx,iter,d,
587  &neq[1],z,istartnmd,iendnmd,&nmd,&cyclicsymmetry,
588  &neqact,&igeneralizedforce));
589  SFREE(aa);SFREE(bb);SFREE(eiga);SFREE(eigb);SFREE(iter);SFREE(zc);
590 
591  }
592 
593 /* sorting the eigenvalues and eigenmodes according to the size of the
594  eigenvalues */
595 
596  NNEW(ipev,ITG,nev);
597  NNEW(eigxr,double,nev);
598  NNEW(aa,double,2*nev);
599  NNEW(bb,double,4*nev*nev);
600 
601  FORTRAN(sortev,(&nev,&nmd,eigxx,&cyclicsymmetry,xx,
602  eigxr,ipev,istartnmd,iendnmd,aa,bb));
603 
604  SFREE(ipev);SFREE(eigxr);SFREE(aa);SFREE(bb);
605 
606  /* storing the eigenvalues in the .dat file */
607 
608  if(cyclicsymmetry){
609  FORTRAN(writeevcscomplex,(eigxx,&nev,nm,&fmin,&fmax));
610  }else{
611  FORTRAN(writeevcomplex,(eigxx,&nev,&fmin,&fmax));
612  }
613 
614  /* storing the participation factors */
615 
616  NNEW(eigxr,double,nev);
617  NNEW(eigxi,double,nev);
618  if(nherm==1){
619  NNEW(eig,double,nev);
620  for(l=0;l<nev;l++){
621  if(d[l]<0.){
622  eig[l]=0.;
623  }else{
624  eig[l]=sqrt(d[l]);
625  }
626  }
627  }else{
628 
629  NNEW(eig,double,2*nev);
630  for(l=0;l<nev;l++){
631  eig[2*l]=sqrt(sqrt(d[2*l]*d[2*l]+d[2*l+1]*d[2*l+1])+d[2*l])/sqrt(2.);
632  eig[2*l+1]=sqrt(sqrt(d[2*l]*d[2*l]+d[2*l+1]*d[2*l+1])-d[2*l])/sqrt(2.);
633  if(d[2*l+1]<0.) eig[2*l+1]=-eig[2*l+1];
634  }
635  }
636  for(l=0;l<nev;l++){
637  mode=l+1;
638  for(k=0;k<nev;k++){
639  eigxr[k]=xx[2*l*nev+2*k];
640  eigxi[k]=xx[2*l*nev+2*k+1];
641  }
642  FORTRAN(writepf,(eig,eigxr,eigxi,&zero,&nev,&mode,&nherm));
643  }
644  SFREE(eigxr);SFREE(eigxi);SFREE(eig);SFREE(d);
645 
646  if(cyclicsymmetry){
647 
648  /* assembling the new eigenmodes */
649 
650  /* storage in zz: per eigenmode first the complete real part of
651  the eigenvector, then the complete imaginary part */
652 
653  NNEW(zz,double,(long long)2*nev*neqact);
654  for(l=0;l<nev;l++){
655  for(i=0;i<neqact;i++){
656  for(k=0;k<nev;k++){
657 
658  /* real part */
659 
660  zz[(long long)2*l*neqact+i]+=
661  (xx[2*l*nev+2*k]*z[(long long)2*k*neqact+i]
662  -xx[2*l*nev+2*k+1]*z[(long long)(2*k+1)*neqact+i]);
663 
664  /* imaginary part */
665 
666  zz[(long long)(2*l+1)*neqact+i]+=
667  (xx[2*l*nev+2*k]*z[(long long)(2*k+1)*neqact+i]
668  +xx[2*l*nev+2*k+1]*z[(long long)2*k*neqact+i]);
669  }
670  }
671  }
672 
673  /* calculating the scalar product of all old eigenmodes with
674  all new eigenmodes => nev x nev matrix */
675 
676  NNEW(xmac,double,nev*nev);
677  NNEW(xmaccpx,double,4*nev*nev);
678  NNEW(bett,double,nev);
679  NNEW(betm,double,nev);
680  FORTRAN(calcmac,(&neq[1],z,zz,&nev,xmac,xmaccpx,istartnmd,
681  iendnmd,&nmd,&cyclicsymmetry,&neqact,bett,betm));
682  FORTRAN(writemaccs,(xmac,&nev,nm));
683 
684  SFREE(xmac);SFREE(bett);SFREE(betm);SFREE(xmaccpx);
685  SFREE(z);
686 
687  /* normalizing the eigenmodes */
688 
689  NNEW(z,double,neq[1]);
690  for(l=0;l<nev;l++){
691  sum=0.;
692  DMEMSET(z,0,neq[1],0.);
693  FORTRAN(op,(&neq[1],&zz[l*neq[1]],z,adb,aub,jq,irow));
694  for(k=0;k<neq[1];k++){
695  sum+=zz[l*neq[1]+k]*z[k];
696  }
697 
698  sum=sqrt(sum);
699  for(k=0;k<neq[1];k++){
700  zz[l*neq[1]+k]/=sum;
701  }
702  }
703  SFREE(z);
704 
705  /* calculating the mass-weighted internal products (eigenvectors are not
706  necessarily orthogonal, since the matrix of the eigenvalue problem is
707  not necessarily Hermitian)
708  = orthogonality matrices */
709 
710  if(mei[3]==1){
711 
712  NNEW(xmr,double,nev*nev);
713  NNEW(xmi,double,nev*nev);
714  NNEW(z,double,neq[1]);
715  NNEW(zi,double,neq[1]);
716 
717  for(l=0;l<nev;l++){
718  DMEMSET(z,0,neq[1],0.);
719  FORTRAN(op,(&neq[1],&zz[l*neq[1]],z,adb,aub,jq,irow));
720  for(m=l;m<nev;m++){
721  for(k=0;k<neq[1];k++){
722  xmr[l*nev+m]+=zz[m*neq[1]+k]*z[k];
723  }
724  }
725 
726  memcpy(&zi[0],&zz[(2*l+1)*neqact],sizeof(double)*neqact);
727  for(k=0;k<neqact;k++){zi[neqact+k]=-zz[2*l*neqact+k];}
728  DMEMSET(z,0,neq[1],0.);
729  FORTRAN(op,(&neq[1],zi,z,adb,aub,jq,irow));
730  for(m=l;m<nev;m++){
731  for(k=0;k<neq[1];k++){
732  xmi[l*nev+m]+=zz[m*neq[1]+k]*z[k];
733  }
734  }
735  }
736 
737  /* Hermitian part of the matrix */
738 
739  for(l=0;l<nev;l++){
740  for(m=0;m<l;m++){
741  xmr[l*nev+m]=xmr[m*nev+l];
742  xmi[l*nev+m]=-xmi[m*nev+l];
743  }
744  }
745 
746  for(l=0;l<nev;l++){
747  for(m=0;m<nev;m++){
748  printf(" %f",xmr[m*nev+l]);
749  }
750  printf("\n");
751  }
752  printf("\n");
753  for(l=0;l<nev;l++){
754  for(m=0;m<nev;m++){
755  printf(" %f",xmi[m*nev+l]);
756  }
757  printf("\n");
758  }
759  SFREE(z);SFREE(zi);
760  }
761 
762  }else{
763 
764  /* no cyclic symmmetry */
765 
766  /* assembling the new eigenmodes */
767 
768  NNEW(zz,double,2*nev*neq[1]);
769  for(l=0;l<nev;l++){
770  for(j=0;j<2;j++){
771  for(i=0;i<neq[1];i++){
772  for(k=0;k<nev;k++){
773  zz[(2*l+j)*neq[1]+i]+=xx[2*l*nev+2*k+j]*
774  z[(long long)k*neq[1]+i];
775  }
776  }
777  }
778  }
779 
780  /* calculating the scalar product of all old eigenmodes with
781  all new eigenmodes => nev x nev matrix */
782 
783  NNEW(xmac,double,nev*nev);
784  NNEW(xmaccpx,double,4*nev*nev);
785  NNEW(bett,double,nev);
786  NNEW(betm,double,nev);
787  FORTRAN(calcmac,(&neq[1],z,zz,&nev,xmac,xmaccpx,istartnmd,
788  iendnmd,&nmd,&cyclicsymmetry,&neqact,bett,betm));
789  FORTRAN(writemac,(xmac,&nev));
790  SFREE(xmac);SFREE(bett);SFREE(betm);SFREE(xmaccpx);
791 
792  SFREE(z);
793 
794  /* normalizing the eigenmodes */
795 
796  NNEW(z,double,neq[1]);
797  for(l=0;l<nev;l++){
798  sum=0.;
799 
800  /* Ureal^T*M*Ureal */
801 
802  DMEMSET(z,0,neq[1],0.);
803  FORTRAN(op,(&neq[1],&zz[2*l*neq[1]],z,adb,aub,jq,irow));
804  for(k=0;k<neq[1];k++){
805  sum+=zz[2*l*neq[1]+k]*z[k];
806  }
807 
808  /* Uimag^T*M*Uimag */
809 
810  DMEMSET(z,0,neq[1],0.);
811  FORTRAN(op,(&neq[1],&zz[(2*l+1)*neq[1]],z,adb,aub,jq,irow));
812  for(k=0;k<neq[1];k++){
813  sum+=zz[(2*l+1)*neq[1]+k]*z[k];
814  }
815 
816  sum=sqrt(sum);
817  for(k=0;k<2*neq[1];k++){
818  zz[2*l*neq[1]+k]/=sum;
819  }
820  }
821  SFREE(z);
822 
823  /* calculating the mass-weighted internal products (eigenvectors are not
824  necessarily orthogonal, since the matrix of the eigenvalue problem is
825  not necessarily symmetric)
826  = orthogonality matrices */
827 
828  if(mei[3]==1){
829 
830  NNEW(xmr,double,nev*nev);
831  NNEW(xmi,double,nev*nev);
832  NNEW(z,double,neq[1]);
833 
834  for(l=0;l<nev;l++){
835  sum=0.;
836 
837  /* M*Ureal */
838 
839  DMEMSET(z,0,neq[1],0.);
840  FORTRAN(op,(&neq[1],&zz[2*l*neq[1]],z,adb,aub,jq,irow));
841 
842  /* Ureal^T*M*Ureal and Uimag^T*M*Ureal */
843 
844  for(m=l;m<nev;m++){
845  for(k=0;k<neq[1];k++){
846  xmr[l*nev+m]+=zz[2*m*neq[1]+k]*z[k];
847  }
848  for(k=0;k<neq[1];k++){
849  xmi[l*nev+m]-=zz[(2*m+1)*neq[1]+k]*z[k];
850  }
851  }
852 
853  /* M*Uimag */
854 
855  DMEMSET(z,0,neq[1],0.);
856  FORTRAN(op,(&neq[1],&zz[(2*l+1)*neq[1]],z,adb,aub,jq,irow));
857 
858  /* Ureal^T*M*Uimag and Uimag^T*M*Uimag */
859 
860  for(m=l;m<nev;m++){
861  for(k=0;k<neq[1];k++){
862  xmr[l*nev+m]+=zz[(2*m+1)*neq[1]+k]*z[k];
863  }
864  for(k=0;k<neq[1];k++){
865  xmi[l*nev+m]+=zz[2*m*neq[1]+k]*z[k];
866  }
867  }
868  }
869 
870  /* Hermitian part of the matrix */
871 
872  for(l=0;l<nev;l++){
873  for(m=0;m<l;m++){
874  xmr[l*nev+m]=xmr[m*nev+l];
875  xmi[l*nev+m]=-xmi[m*nev+l];
876  }
877  }
878 
879  for(l=0;l<nev;l++){
880  for(m=0;m<nev;m++){
881  printf(" %f",xmr[m*nev+l]);
882  }
883  printf("\n");
884  }
885  printf("\n");
886  for(l=0;l<nev;l++){
887  for(m=0;m<nev;m++){
888  printf(" %f",xmi[m*nev+l]);
889  }
890  printf("\n");
891  }
892  SFREE(z);
893  }
894 
895  }
896 
897  /*storing new eigenmodes and eigenvalues to *.eig-file for later use in
898  steady states dynamic analysis*/
899 
900  if(mei[3]==1){
901 
902  nherm=0;
903 
904  if(!cyclicsymmetry){
905  if((f1=fopen(fneig,"wb"))==NULL){
906  printf("*ERROR in complexfreq: cannot open eigenvalue file for writing...");
907  exit(0);
908  }
909 
910  /* storing a zero as indication that this was not a
911  cyclic symmetry calculation */
912 
913  if(fwrite(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
914  printf("*ERROR in complexfreq saving the cyclic symmetry flag to the eigenvalue file...");
915  exit(0);
916  }
917 
918  /* not Hermitian */
919 
920  if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
921  printf("*ERROR in complexfreq saving the Hermitian flag to the eigenvalue file...");
922  exit(0);
923  }
924 
925  /* storing the number of eigenvalues */
926 
927  if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
928  printf("*ERROR in complexfreq saving the number of eigenfrequencies to the eigenvalue file...");
929  exit(0);
930  }
931 
932  /* the eigenfrequencies are stored as (radians/time)**2
933  squaring the complexe eigenvalues first */
934 
935  NNEW(eigx,double,2*nev);
936  for(i=0;i<nev;i++){
937  eigx[2*i]=eigxx[2*i]*eigxx[2*i]-eigxx[2*i+1]*eigxx[2*i+1];
938  eigx[2*i+1]=2.*eigxx[2*i]*eigxx[2*i+1];
939  }
940 
941  if(fwrite(eigx,sizeof(double),2*nev,f1)!=2*nev){
942  printf("*ERROR in complexfreq saving the eigenfrequencies to the eigenvalue file...");
943  exit(0);
944  }
945 
946  SFREE(eigx);
947 
948  /* storing the stiffness matrix */
949 
950  if(fwrite(ad,sizeof(double),neq[1],f1)!=neq[1]){
951  printf("*ERROR in complexfreq saving the diagonal of the stiffness matrix to the eigenvalue file...");
952  exit(0);
953  }
954  if(fwrite(au,sizeof(double),nzs[2],f1)!=nzs[2]){
955  printf("*ERROR in complexfreq saving the off-diagonal entries of the stiffness matrix to the eigenvalue file...");
956  exit(0);
957  }
958 
959  /* storing the mass matrix */
960 
961  if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
962  printf("*ERROR in complexfreq saving the diagonal of the mass matrix to the eigenvalue file...");
963  exit(0);
964  }
965  if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
966  printf("*ERROR in complexfreq saving the off-diagonal entries of the mass matrix to the eigenvalue file...");
967  exit(0);
968  }
969 
970  /* storing the eigenvectors */
971 
972  lfin=0;
973  lint=0;
974  for(j=0;j<nev;++j){
975  lint=lfin;
976  lfin=lfin+2*neq[1];
977  if(fwrite(&zz[lint],sizeof(double),2*neq[1],f1)!=2*neq[1]){
978  printf("*ERROR in complexfreq saving the eigenvectors to the eigenvalue file...");
979  exit(0);
980  }
981  }
982 
983  /* storing the orthogonality matrices */
984 
985  if(fwrite(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
986  printf("*ERROR in complexfreq saving the real orthogonality matrix to the eigenvalue file...");
987  exit(0);
988  }
989 
990  if(fwrite(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
991  printf("*ERROR in complexfreq saving the imaginary orthogonality matrix to the eigenvalue file...");
992  exit(0);
993  }
994 
995  }else{
996 
997  /* opening .eig file for writing */
998 
999  if((f1=fopen(fneig,"wb"))==NULL){
1000  printf("*ERROR in complexfreq: cannot open eigenvalue file for writing...");
1001  exit(0);
1002  }
1003  /* storing a one as indication that this was a
1004  cyclic symmetry calculation */
1005 
1006  if(fwrite(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
1007  printf("*ERROR in complexfreq saving the cyclic symmetry flag to the eigenvalue file...");
1008  exit(0);
1009  }
1010 
1011  /* not Hermitian */
1012 
1013  if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
1014  printf("*ERROR in complexfreq saving the Hermitian flag to the eigenvalue file...");
1015  exit(0);
1016  }
1017 
1018  if(fwrite(&nm[0],sizeof(ITG),1,f1)!=1){
1019  printf("*ERROR in complexfreq saving the nodal diameter to the eigenvalue file...");
1020  exit(0);
1021  }
1022  if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
1023  printf("*ERROR in complexfreq saving the number of eigenvalues to the eigenvalue file...");
1024  exit(0);
1025  }
1026 
1027  /* the eigenfrequencies are stored as (radians/time)**2
1028  squaring the complexe eigenvalues first */
1029 
1030  NNEW(eigx,double,2*nev);
1031  for(i=0;i<nev;i++){
1032  eigx[2*i]=eigxx[2*i]*eigxx[2*i]-eigxx[2*i+1]*eigxx[2*i+1];
1033  eigx[2*i+1]=2.*eigxx[2*i]*eigxx[2*i+1];
1034  }
1035 
1036  if(fwrite(eigx,sizeof(double),2*nev,f1)!=2*nev){
1037  printf("*ERROR in complexfreq saving the eigenfrequencies to the eigenvalue file...");
1038  exit(0);
1039  }
1040 
1041  SFREE(eigx);
1042 
1043  /* storing the mass matrix */
1044 
1045  if(fwrite(adb,sizeof(double),*neq,f1)!=*neq){
1046  printf("*ERROR in complexfreq saving the diagonal of the mass matrix to the eigenvalue file...");
1047  exit(0);
1048  }
1049  if(fwrite(aub,sizeof(double),*nzs,f1)!=*nzs){
1050  printf("*ERROR in complexfreq saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
1051  exit(0);
1052  }
1053 
1054  /* storing the eigenvectors */
1055 
1056  lfin=0;
1057  for(j=0;j<nev;++j){
1058  lint=lfin;
1059  lfin=lfin+neq[1];
1060  if(fwrite(&zz[lint],sizeof(double),neq[1],f1)!=neq[1]){
1061  printf("*ERROR in complexfreq saving the eigenvectors to the eigenvalue file...");
1062  exit(0);
1063  }
1064  }
1065 
1066  /* storing the orthogonality matrices */
1067 
1068  if(fwrite(xmr,sizeof(double),nev*nev,f1)!=nev*nev){
1069  printf("*ERROR in complexfreq saving the real orthogonality matrix to the eigenvalue file...");
1070  exit(0);
1071  }
1072 
1073  if(fwrite(xmi,sizeof(double),nev*nev,f1)!=nev*nev){
1074  printf("*ERROR in complexfreq saving the imaginary orthogonality matrix to the eigenvalue file...");
1075  exit(0);
1076  }
1077  }
1078  SFREE(xmr);SFREE(xmi);
1079 
1080  fclose(f1);
1081  }
1082 
1083  SFREE(adb);SFREE(aub);
1084  if(!cyclicsymmetry){SFREE(ad);SFREE(au);}
1085 
1086  /* calculating the displacements and the stresses and storing */
1087  /* the results in frd format for each valid eigenmode */
1088 
1089  NNEW(v,double,2*mt**nk);
1090  NNEW(fn,double,2*mt**nk);
1091  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1653],"MAXS")==0)||
1092  (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1093  (strcmp1(&filab[1044],"ERR ")==0))
1094  NNEW(stn,double,12**nk);
1095 
1096  if((strcmp1(&filab[261],"E ")==0)||(strcmp1(&filab[2523],"MAXE")==0))
1097  NNEW(een,double,12**nk);
1098  if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,2**nk);
1099  if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,12**nk);
1100 
1101  NNEW(inum,ITG,*nk);
1102  NNEW(stx,double,2*6*mi[0]**ne);
1103 
1104  NNEW(coefmpcnew,double,*mpcend);
1105 
1106  NNEW(cot,double,3**nk*ngraph);
1107  if(*ntrans>0){NNEW(inotrt,ITG,2**nk*ngraph);}
1108  if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0))
1109 
1110 // real and imaginary part of the displacements
1111 
1112  NNEW(vt,double,2*mt**nk*ngraph);
1113  if(strcmp1(&filab[87],"NT ")==0)
1114  NNEW(t1t,double,*nk*ngraph);
1115  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1116  (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1117 
1118 // real and imaginary part of the stresses
1119 
1120  NNEW(stnt,double,2*6**nk*ngraph);
1121  if(strcmp1(&filab[261],"E ")==0)
1122  NNEW(eent,double,2*6**nk*ngraph);
1123  if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0))
1124 
1125 // real and imaginary part of the forces
1126 
1127  NNEW(fnt,double,2*mt**nk*ngraph);
1128  if(strcmp1(&filab[522],"ENER")==0)
1129  NNEW(enernt,double,*nk*ngraph);
1130  if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1131  NNEW(stxt,double,2*6*mi[0]**ne*ngraph);
1132  if(strcmp1(&filab[2697],"ME ")==0)
1133  NNEW(emnt,double,2*6**nk*ngraph);
1134 
1135  NNEW(kont,ITG,*nkon*ngraph);
1136  NNEW(ipkont,ITG,*ne*ngraph);
1137  for(l=0;l<*ne*ngraph;l++)ipkont[l]=-1;
1138  NNEW(lakont,char,8**ne*ngraph);
1139  NNEW(inumt,ITG,*nk*ngraph);
1140  NNEW(ielmatt,ITG,mi[2]**ne*ngraph);
1141 
1142  nkt=ngraph**nk;
1143  net=ngraph**ne;
1144 
1145  /* copying the coordinates of the first sector */
1146 
1147  for(l=0;l<3**nk;l++){cot[l]=co[l];}
1148  if(*ntrans>0){for(l=0;l<*nk;l++){inotrt[2*l]=inotr[2*l];}}
1149  for(l=0;l<*nkon;l++){kont[l]=kon[l];}
1150  for(l=0;l<*ne;l++){ipkont[l]=ipkon[l];}
1151  for(l=0;l<8**ne;l++){lakont[l]=lakon[l];}
1152  for(l=0;l<*ne;l++){ielmatt[mi[2]*l]=ielmat[mi[2]*l];}
1153 
1154  /* generating the coordinates for the other sectors */
1155 
1156  if(cyclicsymmetry){
1157 
1158  icntrl=1;
1159 
1160  FORTRAN(rectcyl,(cot,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1161 
1162  for(jj=0;jj<*mcs;jj++){
1163  is=cs[17*jj+4];
1164  for(i=1;i<is;i++){
1165 
1166  theta=i*2.*pi/cs[17*jj];
1167 
1168  for(l=0;l<*nk;l++){
1169  if(inocs[l]==jj){
1170  cot[3*l+i*3**nk]=cot[3*l];
1171  cot[1+3*l+i*3**nk]=cot[1+3*l]+theta;
1172  cot[2+3*l+i*3**nk]=cot[2+3*l];
1173  if(*ntrans>0){inotrt[2*l+i*2**nk]=inotrt[2*l];}
1174  }
1175  }
1176  for(l=0;l<*nkon;l++){kont[l+i**nkon]=kon[l]+i**nk;}
1177  for(l=0;l<*ne;l++){
1178  if(ielcs[l]==jj){
1179  if(ipkon[l]>=0){
1180  ipkont[l+i**ne]=ipkon[l]+i**nkon;
1181  ielmatt[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
1182  for(l1=0;l1<8;l1++){
1183  l2=8*l+l1;
1184  lakont[l2+i*8**ne]=lakon[l2];
1185  }
1186  }
1187  }
1188  }
1189  }
1190  }
1191 
1192  icntrl=-1;
1193 
1194  FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1195  &imag,mi,emnt));
1196  }
1197 
1198  /* check that the tensor fields which are extrapolated from the
1199  integration points are requested in global coordinates */
1200 
1201  if(strcmp1(&filab[174],"S ")==0){
1202  if((strcmp1(&filab[179],"L")==0)&&(*norien>0)){
1203  printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculations\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1204  strcpy1(&filab[179],"G",1);
1205  }
1206  }
1207 
1208  if(strcmp1(&filab[261],"E ")==0){
1209  if((strcmp1(&filab[266],"L")==0)&&(*norien>0)){
1210  printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1211  strcpy1(&filab[266],"G",1);
1212  }
1213  }
1214 
1215  if(strcmp1(&filab[1479],"PHS ")==0){
1216  if((strcmp1(&filab[1484],"L")==0)&&(*norien>0)){
1217  printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1218  strcpy1(&filab[1484],"G",1);
1219  }
1220  }
1221 
1222  if(strcmp1(&filab[1653],"MAXS")==0){
1223  if((strcmp1(&filab[1658],"L")==0)&&(*norien>0)){
1224  printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1225  strcpy1(&filab[1658],"G",1);
1226  }
1227  }
1228 
1229  if(strcmp1(&filab[2523],"MAXE")==0){
1230  if((strcmp1(&filab[2528],"L")==0)&&(*norien>0)){
1231  printf("\n*WARNING in complexfreq: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1232  strcpy1(&filab[1658],"G",1);
1233  }
1234  }
1235 
1236  /* allocating fields for magnitude and phase information of
1237  displacements and stresses */
1238 
1239  if(strcmp1(&filab[870],"PU")==0){
1240  NNEW(vr,double,mt*nkt);
1241  NNEW(vi,double,mt*nkt);
1242  }
1243 
1244  if(strcmp1(&filab[1479],"PHS")==0){
1245  NNEW(stnr,double,6*nkt);
1246  NNEW(stni,double,6*nkt);
1247  }
1248 
1249  if(strcmp1(&filab[1566],"MAXU")==0){
1250  NNEW(vmax,double,4*nkt);
1251  }
1252 
1253  if(strcmp1(&filab[1653],"MAXS")==0){
1254  NNEW(stnmax,double,7*nkt);
1255  }
1256 
1257  if(strcmp1(&filab[2523],"MAXE")==0){
1258  NNEW(eenmax,double,7*nkt);
1259  }
1260 
1261  /* storing the results */
1262 
1263  if(!cyclicsymmetry) (neq[1])*=2;
1264 
1265  lfin=0;
1266  for(j=0;j<nev;++j){
1267  lint=lfin;
1268  lfin=lfin+neq[1];
1269 
1270  /* calculating the cosine and sine */
1271 
1272  for(k=0;k<*mcs;k++){
1273  theta=nm[j]*2.*pi/cs[17*k];
1274  cs[17*k+14]=cos(theta);
1275  cs[17*k+15]=sin(theta);
1276  }
1277 
1278  if(*nprint>0)FORTRAN(writehe,(&j));
1279 
1280  NNEW(eei,double,6*mi[0]**ne);
1281  if(*nener==1){
1282  NNEW(stiini,double,6*mi[0]**ne);
1283  NNEW(enerini,double,mi[0]**ne);}
1284 
1285  DMEMSET(v,0,2*mt**nk,0.);
1286 
1287  for(k=0;k<neq[1];k+=neq[1]/2){
1288 
1289  for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
1290 
1291  if(k==0) {kk=0;kkv=0;kk6=0;kkx=0;if(*nprint>0)FORTRAN(writere,());}
1292  else {kk=*nk;kkv=mt**nk;kk6=6**nk;kkx=6*mi[0]**ne;
1293  if(*nprint>0)FORTRAN(writeim,());}
1294 
1295  /* generating the cyclic MPC's (needed for nodal diameters
1296  different from 0 */
1297 
1298  for(i=0;i<*nmpc;i++){
1299  index=ipompc[i]-1;
1300  /* check whether thermal mpc */
1301  if(nodempc[3*index+1]==0) continue;
1302  coefmpcnew[index]=coefmpc[index];
1303  while(1){
1304  index=nodempc[3*index+2];
1305  if(index==0) break;
1306  index--;
1307 
1308  icomplex=0;
1309  inode=nodempc[3*index];
1310  if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
1311  icomplex=atoi(&labmpc[20*i+6]);}
1312  else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
1313  for(ij=0;ij<*mcs;ij++){
1314  lprev=cs[ij*17+13];
1315  ilength=cs[ij*17+3];
1316  FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
1317  if(id!=0){
1318  if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
1319  }
1320  }
1321  }
1322 
1323  if(icomplex!=0){
1324  idir=nodempc[3*index+1];
1325  idof=nactdof[mt*(inode-1)+idir]-1;
1326  if(idof==-1){xreal=1.;ximag=1.;}
1327  else{xreal=zz[lint+idof];ximag=zz[lint+idof+neq[1]/2];}
1328  if(k==0) {
1329  if(fabs(xreal)<1.e-30)xreal=1.e-30;
1330  coefmpcnew[index]=coefmpc[index]*
1331  (cs[17*(icomplex-1)+14]+ximag/xreal*cs[17*(icomplex-1)+15]);}
1332  else {
1333  if(fabs(ximag)<1.e-30)ximag=1.e-30;
1334  coefmpcnew[index]=coefmpc[index]*
1335  (cs[17*(icomplex-1)+14]-xreal/ximag*cs[17*(icomplex-1)+15]);}
1336  }
1337  else{coefmpcnew[index]=coefmpc[index];}
1338  }
1339  }
1340 
1341  if(*iperturb==0){
1342  results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1343  &stx[kkx],elcon,
1344  nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1345  norien,orab,ntmat_,t0,t0,ithermal,
1346  prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1347  f,&fn[kkv],nactdof,&iout,qa,vold,&zz[lint+k],
1348  nodeboun,ndirboun,xboun,nboun,ipompc,
1349  nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1350  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1351  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1352  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1353  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1354  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1355  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1356  &ne0,xforc,nforc,thicke,shcon,nshcon,
1357  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1358  &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1359  islavsurf,ielprop,prop);}
1360  else{
1361  results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1362  &stx[kkx],elcon,
1363  nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1364  norien,orab,ntmat_,t0,t1old,ithermal,
1365  prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
1366  f,&fn[kkv],nactdof,&iout,qa,vold,&zz[lint+k],
1367  nodeboun,ndirboun,xboun,nboun,ipompc,
1368  nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1369  &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1370  xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1371  ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1372  xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1373  ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1374  nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1375  &ne0,xforc,nforc,thicke,shcon,nshcon,
1376  sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1377  &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1378  islavsurf,ielprop,prop);
1379  }
1380 
1381  }
1382  SFREE(eei);
1383  if(*nener==1){SFREE(stiini);SFREE(enerini);}
1384 
1385  /* changing the basic results into cylindrical coordinates
1386  (only for cyclic symmetry structures */
1387 
1388  for(l=0;l<*nk;l++){inumt[l]=inum[l];}
1389 
1390  if(cyclicsymmetry){
1391  icntrl=2;imag=1;
1392  FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1393  }
1394 
1395  /* copying the basis results (real and imaginary) into
1396  a new field */
1397 
1398  if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1399  for(l=0;l<mt**nk;l++){vt[l]=v[l];}
1400  for(l=0;l<mt**nk;l++){vt[l+mt**nk*ngraph]=v[l+mt**nk];}}
1401  if(strcmp1(&filab[87],"NT ")==0)
1402  for(l=0;l<*nk;l++){t1t[l]=t1[l];};
1403  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1404  for(l=0;l<6**nk;l++){stnt[l]=stn[l];}
1405  for(l=0;l<6**nk;l++){stnt[l+6**nk*ngraph]=stn[l+6**nk];}}
1406  if(strcmp1(&filab[261],"E ")==0){
1407  for(l=0;l<6**nk;l++){eent[l]=een[l];};
1408  for(l=0;l<6**nk;l++){eent[l+6**nk*ngraph]=een[l+6**nk];}}
1409  if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1410  for(l=0;l<mt**nk;l++){fnt[l]=fn[l];}
1411  for(l=0;l<mt**nk;l++){fnt[l+mt**nk*ngraph]=fn[l+mt**nk];}}
1412  if(strcmp1(&filab[522],"ENER")==0)
1413  for(l=0;l<*nk;l++){enernt[l]=enern[l];};
1414  if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)){
1415  for(l=0;l<6*mi[0]**ne;l++){stxt[l]=stx[l];}
1416  for(l=0;l<6*mi[0]**ne;l++){stxt[l+6*mi[0]**ne*ngraph]=stx[l+6*mi[0]**ne];}}
1417  if(strcmp1(&filab[2697],"ME ")==0){
1418  for(l=0;l<6**nk;l++){emnt[l]=emn[l];};
1419  for(l=0;l<6**nk;l++){emnt[l+6**nk*ngraph]=emn[l+6**nk];}}
1420 
1421  /* mapping the results to the other sectors
1422  (only for cyclic symmetric structures */
1423 
1424  if(cyclicsymmetry){
1425 
1426  for(jj=0;jj<*mcs;jj++){
1427  ilength=cs[17*jj+3];
1428  is=cs[17*jj+4];
1429  lprev=cs[17*jj+13];
1430  for(i=1;i<is;i++){
1431 
1432  for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
1433 
1434  theta=i*nm[j]*2.*pi/cs[17*jj];
1435  ctl=cos(theta);
1436  stl=sin(theta);
1437 
1438  if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1439  for(l1=0;l1<*nk;l1++){
1440  if(inocs[l1]==jj){
1441 
1442  /* check whether node lies on axis */
1443 
1444  ml1=-l1-1;
1445  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1446  if(id!=0){
1447  if(ics[lprev+id-1]==ml1){
1448  for(l2=0;l2<4;l2++){
1449  l=mt*l1+l2;
1450  vt[l+mt**nk*i]=v[l];
1451  }
1452  continue;
1453  }
1454  }
1455  for(l2=0;l2<4;l2++){
1456  l=mt*l1+l2;
1457  vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
1458  }
1459 
1460  }
1461  }
1462  }
1463 
1464  /* imaginary part of the displacements in cylindrical
1465  coordinates */
1466 
1467  if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1468  for(l1=0;l1<*nk;l1++){
1469  if(inocs[l1]==jj){
1470 
1471  /* check whether node lies on axis */
1472 
1473  ml1=-l1-1;
1474  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1475  if(id!=0){
1476  if(ics[lprev+id-1]==ml1){
1477  for(l2=0;l2<4;l2++){
1478  l=mt*l1+l2;
1479  vt[l+mt**nk*(i+ngraph)]=v[l+mt**nk];
1480  }
1481  continue;
1482  }
1483  }
1484  for(l2=0;l2<4;l2++){
1485  l=mt*l1+l2;
1486  vt[l+mt**nk*(i+ngraph)]=stl*v[l]+ctl*v[l+mt**nk];
1487  }
1488  }
1489  }
1490  }
1491 
1492  if(strcmp1(&filab[87],"NT ")==0){
1493  for(l=0;l<*nk;l++){
1494  if(inocs[l]==jj) t1t[l+*nk*i]=t1[l];
1495  }
1496  }
1497 
1498  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1499  for(l1=0;l1<*nk;l1++){
1500  if(inocs[l1]==jj){
1501 
1502  /* check whether node lies on axis */
1503 
1504  ml1=-l1-1;
1505  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1506  if(id!=0){
1507  if(ics[lprev+id-1]==ml1){
1508  for(l2=0;l2<6;l2++){
1509  l=6*l1+l2;
1510  stnt[l+6**nk*i]=stn[l];
1511  }
1512  continue;
1513  }
1514  }
1515  for(l2=0;l2<6;l2++){
1516  l=6*l1+l2;
1517  stnt[l+6**nk*i]=ctl*stn[l]-stl*stn[l+6**nk];
1518  }
1519  }
1520  }
1521  }
1522 
1523  /* imaginary part of the stresses in cylindrical
1524  coordinates */
1525 
1526  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1527  for(l1=0;l1<*nk;l1++){
1528  if(inocs[l1]==jj){
1529 
1530  /* check whether node lies on axis */
1531 
1532  ml1=-l1-1;
1533  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1534  if(id!=0){
1535  if(ics[lprev+id-1]==ml1){
1536  for(l2=0;l2<6;l2++){
1537  l=6*l1+l2;
1538  stnt[l+6**nk*(i+ngraph)]=stn[l+6**nk];
1539  }
1540  continue;
1541  }
1542  }
1543  for(l2=0;l2<6;l2++){
1544  l=6*l1+l2;
1545  stnt[l+6**nk*(i+ngraph)]=stl*stn[l]+ctl*stn[l+6**nk];
1546  }
1547  }
1548  }
1549  }
1550 
1551  if(strcmp1(&filab[261],"E ")==0){
1552  for(l1=0;l1<*nk;l1++){
1553  if(inocs[l1]==jj){
1554 
1555  /* check whether node lies on axis */
1556 
1557  ml1=-l1-1;
1558  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1559  if(id!=0){
1560  if(ics[lprev+id-1]==ml1){
1561  for(l2=0;l2<6;l2++){
1562  l=6*l1+l2;
1563  eent[l+6**nk*i]=een[l];
1564  }
1565  continue;
1566  }
1567  }
1568  for(l2=0;l2<6;l2++){
1569  l=6*l1+l2;
1570  eent[l+6**nk*i]=ctl*een[l]-stl*een[l+6**nk];
1571  }
1572  }
1573  }
1574  }
1575 
1576  /* imaginary part of the strains in cylindrical
1577  coordinates */
1578 
1579  if(strcmp1(&filab[261],"E ")==0){
1580  for(l1=0;l1<*nk;l1++){
1581  if(inocs[l1]==jj){
1582 
1583  /* check whether node lies on axis */
1584 
1585  ml1=-l1-1;
1586  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1587  if(id!=0){
1588  if(ics[lprev+id-1]==ml1){
1589  for(l2=0;l2<6;l2++){
1590  l=6*l1+l2;
1591  eent[l+6**nk*(i+ngraph)]=een[l+6**nk];
1592  }
1593  continue;
1594  }
1595  }
1596  for(l2=0;l2<6;l2++){
1597  l=6*l1+l2;
1598  eent[l+6**nk*(i+ngraph)]=stl*een[l]+ctl*een[l+6**nk];
1599  }
1600  }
1601  }
1602  }
1603 
1604  if(strcmp1(&filab[2697],"ME ")==0){
1605  for(l1=0;l1<*nk;l1++){
1606  if(inocs[l1]==jj){
1607 
1608  /* check whether node lies on axis */
1609 
1610  ml1=-l1-1;
1611  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1612  if(id!=0){
1613  if(ics[lprev+id-1]==ml1){
1614  for(l2=0;l2<6;l2++){
1615  l=6*l1+l2;
1616  emnt[l+6**nk*i]=emn[l];
1617  }
1618  continue;
1619  }
1620  }
1621  for(l2=0;l2<6;l2++){
1622  l=6*l1+l2;
1623  emnt[l+6**nk*i]=ctl*emn[l]-stl*emn[l+6**nk];
1624  }
1625  }
1626  }
1627  }
1628 
1629  /* imaginary part of the mechanical strains in cylindrical
1630  coordinates */
1631 
1632  if(strcmp1(&filab[2697],"ME ")==0){
1633  for(l1=0;l1<*nk;l1++){
1634  if(inocs[l1]==jj){
1635 
1636  /* check whether node lies on axis */
1637 
1638  ml1=-l1-1;
1639  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1640  if(id!=0){
1641  if(ics[lprev+id-1]==ml1){
1642  for(l2=0;l2<6;l2++){
1643  l=6*l1+l2;
1644  emnt[l+6**nk*(i+ngraph)]=emn[l+6**nk];
1645  }
1646  continue;
1647  }
1648  }
1649  for(l2=0;l2<6;l2++){
1650  l=6*l1+l2;
1651  emnt[l+6**nk*(i+ngraph)]=stl*emn[l]+ctl*emn[l+6**nk];
1652  }
1653  }
1654  }
1655  }
1656 
1657  if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1658  for(l1=0;l1<*nk;l1++){
1659  if(inocs[l1]==jj){
1660 
1661  /* check whether node lies on axis */
1662 
1663  ml1=-l1-1;
1664  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1665  if(id!=0){
1666  if(ics[lprev+id-1]==ml1){
1667  for(l2=0;l2<4;l2++){
1668  l=mt*l1+l2;
1669  fnt[l+mt**nk*i]=fn[l];
1670  }
1671  continue;
1672  }
1673  }
1674  for(l2=0;l2<4;l2++){
1675  l=mt*l1+l2;
1676  fnt[l+mt**nk*i]=ctl*fn[l]-stl*fn[l+mt**nk];
1677  }
1678  }
1679  }
1680  }
1681 
1682  /* imaginary part of the forces in cylindrical
1683  coordinates */
1684 
1685  if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1686  for(l1=0;l1<*nk;l1++){
1687  if(inocs[l1]==jj){
1688 
1689  /* check whether node lies on axis */
1690 
1691  ml1=-l1-1;
1692  FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1693  if(id!=0){
1694  if(ics[lprev+id-1]==ml1){
1695  for(l2=0;l2<4;l2++){
1696  l=mt*l1+l2;
1697  fnt[l+mt**nk*(i+ngraph)]=fn[l+mt**nk];
1698  }
1699  continue;
1700  }
1701  }
1702  for(l2=0;l2<4;l2++){
1703  l=mt*l1+l2;
1704  fnt[l+mt**nk*(i+ngraph)]=stl*fn[l]+ctl*fn[l+mt**nk];
1705  }
1706  }
1707  }
1708  }
1709 
1710  if(strcmp1(&filab[522],"ENER")==0){
1711  for(l=0;l<*nk;l++){
1712  if(inocs[l]==jj) enernt[l+*nk*i]=0.;
1713  }
1714  }
1715  }
1716  }
1717 
1718  icntrl=-2;imag=0;
1719 
1720  FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1721  &imag,mi,emnt));
1722 
1723  FORTRAN(rectcylvi,(cot,&vt[mt**nk*ngraph],&fnt[mt**nk*ngraph],
1724  &stnt[6**nk*ngraph],qfnt,&eent[6**nk*ngraph],
1725  cs,&nkt,&icntrl,t,filab,&imag,mi,&emnt[6**nk*ngraph]));
1726 
1727  }
1728 
1729  /* determining magnitude and phase angle for the displacements */
1730 
1731  if(strcmp1(&filab[870],"PU")==0){
1732  for(l1=0;l1<nkt;l1++){
1733  for(l2=0;l2<4;l2++){
1734  l=mt*l1+l2;
1735  vreal=vt[l];
1736  vimag=vt[l+mt**nk*ngraph];
1737  vr[l]=sqrt(vreal*vreal+vimag*vimag);
1738  if(fabs(vreal)<1.e-10){
1739  if(vimag>0){vi[l]=90.;}
1740  else{vi[l]=-90.;}
1741  }
1742  else{
1743  vi[l]=atan(vimag/vreal)*constant;
1744  if(vreal<0) vi[l]+=180.;
1745  }
1746  }
1747  }
1748  }
1749 
1750  /* determining magnitude and phase for the stress */
1751 
1752  if(strcmp1(&filab[1479],"PHS")==0){
1753  for(l1=0;l1<nkt;l1++){
1754  for(l2=0;l2<6;l2++){
1755  l=6*l1+l2;
1756  stnreal=stnt[l];
1757  stnimag=stnt[l+6**nk*ngraph];
1758  stnr[l]=sqrt(stnreal*stnreal+stnimag*stnimag);
1759  if(fabs(stnreal)<1.e-10){
1760  if(stnimag>0){stni[l]=90.;}
1761  else{stni[l]=-90.;}
1762  }
1763  else{
1764  stni[l]=atan(stnimag/stnreal)*constant;
1765  if(stnreal<0) stni[l]+=180.;
1766  }
1767  }
1768  }
1769  }
1770 
1771  ++*kode;
1772 
1773  /* storing the real part of the eigenfrequencies in freq */
1774 
1775  freq=eigxx[2*j]/6.283185308;
1776  if(strcmp1(&filab[1044],"ZZS")==0){
1777  NNEW(neigh,ITG,40*net);
1778  NNEW(ipneigh,ITG,nkt);
1779  }
1780  frd(cot,&nkt,kont,ipkont,lakont,&net,vt,stnt,inumt,nmethod,
1781  kode,filab,eent,t1t,fnt,&freq,epn,ielmatt,matname,enernt,xstaten,
1782  nstate_,istep,&iinc,ithermal,qfn,&j,&nm[j],trab,inotrt,
1783  ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1784  mi,stxt,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&net,
1785  cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emnt,
1786  thicke,jobnamec,output,qfx,cdn,&mortar,cdnr,cdni,nmat);
1787  if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1788  }
1789 
1790  SFREE(xstiff);if(*nbody>0) SFREE(ipobody);
1791  SFREE(cstr);SFREE(zz);SFREE(eigxx);SFREE(xx);
1792 
1793  if(cyclicsymmetry){
1794  SFREE(istartnmd);SFREE(iendnmd);
1795  }else{
1796  (neq[1])/=2;
1797  }
1798 
1799  SFREE(nm);SFREE(coefmpcnew);
1800 
1801  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1653],"MAXS")==0)||
1802  (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1803  (strcmp1(&filab[1044],"ERR ")==0))
1804  SFREE(stn);
1805 
1806  SFREE(v);SFREE(fn);SFREE(inum);SFREE(stx);
1807 
1808  if((strcmp1(&filab[261],"E ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) SFREE(een);
1809  if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1810  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
1811 
1812  if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)) SFREE(vt);
1813  if(strcmp1(&filab[87],"NT ")==0) SFREE(t1t);
1814  if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1815  (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stnt);
1816  if(strcmp1(&filab[261],"E ")==0) SFREE(eent);
1817  if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)) SFREE(fnt);
1818  if(strcmp1(&filab[522],"ENER")==0) SFREE(enernt);
1819  if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stxt);
1820  if(strcmp1(&filab[2697],"ME ")==0) SFREE(emnt);
1821 
1822  SFREE(cot);SFREE(kont);SFREE(ipkont);SFREE(lakont);SFREE(inumt);SFREE(ielmatt);
1823  if(*ntrans>0){SFREE(inotrt);}
1824 
1825  *ialsetp=ialset;
1826  *cop=co;*konp=kon;*ipkonp=ipkon;*lakonp=lakon;*ielmatp=ielmat;
1827  *ielorienp=ielorien;*inotrp=inotr;*nodebounp=nodeboun;
1828  *ndirbounp=ndirboun;*iambounp=iamboun;*xbounp=xboun;
1829  *xbounoldp=xbounold;*ikbounp=ikboun;*ilbounp=ilboun;*nactdofp=nactdof;
1830  *voldp=vold;*emep=eme;*enerp=ener;*ipompcp=ipompc;*nodempcp=nodempc;
1831  *coefmpcp=coefmpc;*labmpcp=labmpc;*ikmpcp=ikmpc;*ilmpcp=ilmpc;
1832  *fmpcp=fmpc;*veoldp=veold;*iamt1p=iamt1;*t0p=t0;*t1oldp=t1old;*t1p=t1;
1833  *stip=sti;
1834 
1835  return;
1836 }
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 writemac(mac, nev)
Definition: writemac.f:19
subroutine coriolissolve(cc, nev, a, b, x, eiga, eigb, eigcorio, iter, d, temp)
Definition: coriolissolve.f:19
subroutine mafillcorio(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)
Definition: mafillcorio.f:19
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
ITG strcpy1(char *s1, const char *s2, ITG length)
Definition: strcpy1.c:24
subroutine op_corio(n, x, y, ad, au, jq, irow)
Definition: op_corio.f:26
subroutine op(n, x, y, ad, au, jq, irow)
Definition: op.f:25
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
subroutine writehe(j)
Definition: writehe.f:19
#define DMEMSET(a, b, c, d)
Definition: CalculiX.h:45
subroutine stop()
Definition: stop.f:19
subroutine rectcylvi(co, v, fn, stn, qfn, een, cs, n, icntrl, t, filab, imag, mi, emn)
Definition: rectcylvi.f:19
subroutine writemaccs(mac, nev, nm)
Definition: writemaccs.f:19
subroutine writepf(d, bjr, bji, freq, nev, mode, nherm)
Definition: writepf.f:19
subroutine writere()
Definition: writere.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 sortev(nev, nmd, eigxx, cyclicsymmetry, x, eigxr, ipev, istartnmd, iendnmd, a, b)
Definition: sortev.f:25
subroutine writeim()
Definition: writeim.f:19
subroutine nident(x, px, n, id)
Definition: nident.f:25
subroutine forcesolve(zc, nev, a, b, x, eiga, eigb, eigxx, iter, d, neq, z, istartnmd, iendnmd, nmd, cyclicsymmetry, neqact, igeneralizedforce)
Definition: forcesolve.f:19
subroutine readforce(zc, neq, nev, nactdof, ikmpc, nmpc, ipompc, nodempc, mi, coefmpc, jobnamef, a, igeneralizedforce)
Definition: readforce.f:19
subroutine calcmac(neq, z, zz, nev, xmac, xmaccpx, istartnmd, iendnmd, nmd, cyclicsymmetry, neqact, bett, betm)
Definition: calcmac.f:45
subroutine writeevcscomplex(x, nx, nm, fmin, fmax)
Definition: writeevcscomplex.f:19
#define ITG
Definition: CalculiX.h:51
#define NNEW(a, b, c)
Definition: CalculiX.h:39
subroutine bodyforce(cbody, ibody, ipobody, nbody, set, istartset, iendset, ialset, inewton, nset, ifreebody, k)
Definition: bodyforce.f:19
subroutine writeevcomplex(x, nx, fmin, fmax)
Definition: writeevcomplex.f:19
subroutine rectcyl(co, v, fn, stn, qfn, een, cs, n, icntrl, t, filab, imag, mi, emn)
Definition: rectcyl.f:19
Hosted by OpenAircraft.com, (Michigan UAV, LLC)