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