1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 Guido Dhondt */
3
4 /* This program is free software; you can redistribute it and/or */
5 /* modify it under the terms of the GNU General Public License as */
6 /* published by the Free Software Foundation(version 2); */
7 /* */
8
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program; if not, write to the Free Software */
16 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
17
18 #ifdef ARPACK
19
20 #include <stdio.h>
21 #include <math.h>
22 #include <stdlib.h>
23 #include <string.h>
24 #include "CalculiX.h"
25 #ifdef SPOOLES
26 #include "spooles.h"
27 #endif
28 #ifdef SGI
29 #include "sgi.h"
30 #endif
31 #ifdef TAUCS
32 #include "tau.h"
33 #endif
34 #ifdef MATRIXSTORAGE
35 #include "matrixstorage.h"
36 #endif
37 #ifdef PARDISO
38 #include "pardiso.h"
39 #endif
40 #ifdef PASTIX
41 #include "pastix.h"
42 #endif
43
arpackcs(double * co,ITG * nk,ITG ** konp,ITG ** ipkonp,char ** lakonp,ITG * ne,ITG * nodeboun,ITG * ndirboun,double * xboun,ITG * nboun,ITG * ipompc,ITG * nodempc,double * coefmpc,char * labmpc,ITG * nmpc,ITG * nodeforc,ITG * ndirforc,double * xforc,ITG * nforc,ITG * nelemload,char * sideload,double * xload,ITG * nload,ITG * nactdof,ITG * icol,ITG * jq,ITG ** irowp,ITG * neq,ITG * nzl,ITG * nmethod,ITG * ikmpc,ITG * ilmpc,ITG * ikboun,ITG * ilboun,double * elcon,ITG * nelcon,double * rhcon,ITG * nrhcon,double * alcon,ITG * nalcon,double * alzero,ITG ** ielmatp,ITG ** ielorienp,ITG * norien,double * orab,ITG * ntmat_,double * t0,double * t1,double * t1old,ITG * ithermal,double * prestr,ITG * iprestr,double * vold,ITG * iperturb,double * sti,ITG * nzs,ITG * kode,ITG * mei,double * fei,char * filab,ITG * iexpl,double * plicon,ITG * nplicon,double * plkcon,ITG * nplkcon,double ** xstatep,ITG * npmat_,char * matname,ITG * mi,ITG * ics,double * cs,ITG * mpcend,ITG * ncmat_,ITG * nstate_,ITG * mcs,ITG * nkon,char * jobnamec,char * output,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,ITG * nprint,char * prlab,char * prset,ITG * nener,ITG * isolver,double * trab,ITG * inotr,ITG * ntrans,double * ttime,double * fmpc,ITG * ipobody,ITG * ibody,double * xbody,ITG * nbody,ITG * nevtot,double * thicke,ITG * nslavs,double * tietol,ITG * mpcinfo,ITG * ntie,ITG * istep,char * tieset,ITG * nintpoint,ITG * mortar,ITG * ifacecount,ITG ** islavsurfp,double ** pslavsurfp,double ** clearinip,ITG * nmat,char * typeboun,ITG * ielprop,double * prop,char * orname,ITG * inewton,double * t0g,double * t1g)44 void arpackcs(double *co, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp,
45 ITG *ne,
46 ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun,
47 ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
48 ITG *nmpc,
49 ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc,
50 ITG *nelemload, char *sideload, double *xload,
51 ITG *nload, ITG *nactdof,
52 ITG *icol, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl,
53 ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun,
54 ITG *ilboun,
55 double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
56 double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
57 ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
58 double *t0, double *t1, double *t1old,
59 ITG *ithermal,double *prestr, ITG *iprestr,
60 double *vold,ITG *iperturb, double *sti, ITG *nzs,
61 ITG *kode, ITG *mei, double *fei,
62 char *filab,
63 ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon,
64 ITG *nplkcon,
65 double **xstatep, ITG *npmat_, char *matname, ITG *mi,
66 ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_,
67 ITG *nstate_, ITG *mcs, ITG *nkon,
68 char *jobnamec, char *output, char *set, ITG *nset,
69 ITG *istartset,
70 ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
71 char *prset, ITG *nener, ITG *isolver, double *trab,
72 ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
73 ITG *ipobody, ITG *ibody, double *xbody, ITG *nbody,
74 ITG *nevtot, double *thicke, ITG *nslavs, double *tietol,
75 ITG *mpcinfo,ITG *ntie,ITG *istep,
76 char *tieset,ITG *nintpoint,ITG *mortar,ITG *ifacecount,
77 ITG **islavsurfp,double **pslavsurfp,double **clearinip,
78 ITG *nmat,char *typeboun,ITG *ielprop,double *prop,
79 char *orname,ITG *inewton,double *t0g,double *t1g){
80
81 /* calls the Arnoldi Package (ARPACK) for cyclic symmetry calculations */
82
83 char bmat[2]="G", which[3]="LM", howmny[2]="A",*lakont=NULL,
84 description[13]=" ",fneig[132]="",filabcp[9]=" ",
85 lakonl[2]=" \0",*lakon=NULL,jobnamef[396]="",*turdir=NULL,
86 *labmpc2=NULL;
87
88 ITG *inum=NULL,k,ido,ldz,iparam[11],ipntr[14],lworkl,idir,nherm=1,
89 info,rvec=1,*select=NULL,lfin,j,lint,iout=1,nm,index,inode,id,i,idof,
90 ielas=0,icmd=0,kk,l,nkt,icntrl,*kont=NULL,*ipkont=NULL,*inumt=NULL,
91 *ielmatt=NULL,net,imag,icomplex,kkv,kk6,iinc=1,nev,ncv,kscale=1,
92 mxiter,lprev,ilength,ij,i1,i2,iel,ielset,node,indexe,nope,ml1,
93 *inocs=NULL,*ielcs=NULL,jj,l1,l2,ngraph,is,jrow,
94 *inotrt=NULL,symmetryflag=0,inputformat=0,ifreebody,
95 mass=1, stiffness=1, buckling=0, rhsi=0, intscheme=0,*ncocon=NULL,
96 coriolis=0,iworsttime,l3,iray,mt,kkx,im,ne0,*integerglob=NULL,
97 *nshcon=NULL,one=1,ncont=0,*itietri=NULL,neq2,
98 *koncont=NULL,ismallsliding=0,*itiefac=NULL,*islavsurf=NULL,
99 *islavnode=NULL,*imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,
100 *imastop=NULL,*iponoels=NULL,*inoels=NULL,*ipe=NULL,*ime=NULL,
101 mpcfree,memmpc_,icascade,maxlenmpc,nkon0,iit=-1,*irow=NULL,nasym=0,
102 kmax1,kmax2,icfd=0,*inomat=NULL,*ipkon=NULL,*kon=NULL,*ielmat=NULL,
103 *ielorien=NULL,*islavact=NULL,*islavsurfold=NULL,nslavs_prev_step,
104 maxprevcontel,iflagact=0,*nmc=NULL,icutb=0,ialeatoric=0,
105 *iponoel=NULL,*inoel=NULL,network=0,ioffr,nrhs=1,
106 ioffrl,igreen=0,mscalmethod=0,kref,*jqw=NULL,*iroww=NULL,nzsw,
107 *islavelinv=NULL,*irowtloc=NULL,*jqtloc=NULL,nboun2,
108 *ndirboun2=NULL,*nodeboun2=NULL,nmpc2,*ipompc2=NULL,*nodempc2=NULL,
109 *ikboun2=NULL,*ilboun2=NULL,*ikmpc2=NULL,*ilmpc2=NULL,mortartrafoflag=0;
110
111 double *stn=NULL,*v=NULL,*resid=NULL,*z=NULL,*workd=NULL,*vr=NULL,
112 *workl=NULL,*d=NULL,sigma,*temp_array=NULL,*vini=NULL,
113 *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[4],*fext=NULL,*emn=NULL,
114 *epn=NULL,*stiini=NULL,*fnr=NULL,*fni=NULL,fnreal,fnimag,*emeini=NULL,
115 *xstateini=NULL,theta=0,pi,*coefmpcnew=NULL,*xstiff=NULL,*vi=NULL,
116 *vt=NULL,*fnt=NULL,*stnt=NULL,*eent=NULL,*cot=NULL,t[3],ctl,stl,
117 *t1t=NULL,freq,*stx=NULL,*enern=NULL,*enernt=NULL,*xstaten=NULL,
118 *eei=NULL,*enerini=NULL,*cocon=NULL,*qfx=NULL,*qfn=NULL,*qfnt=NULL,
119 tol,fmin,fmax,xreal,ximag,*cgr=NULL,*xloadold=NULL,reltime=1.,constant,
120 vreal,vimag,*stnr=NULL,*stni=NULL,stnreal,stnimag,*vmax=NULL,
121 *stnmax=NULL,vl[4],stnl[6],dd,v1,v2,v3,bb,cc,al[3],cm,cn,tt,
122 worstpsmax,vray[3],worstumax,p1[3],p2[3],q[3],tan[3],*springarea=NULL,
123 *stxt=NULL,*eenmax=NULL,eenl[6],*emnt=NULL,*clearini=NULL,
124 *doubleglob=NULL,*shcon=NULL,*cg=NULL,*straight=NULL,*cdn=NULL,
125 *xmastnor=NULL,*areaslav=NULL,*dc=NULL,*di=NULL,*xnoels=NULL,
126 *workev=NULL,*temp_array2=NULL,*ener=NULL,*xstate=NULL,cdnimag,
127 sigmai=0,amp,ampmax,*zstorage=NULL,*au=NULL,*ad=NULL,cdnreal,
128 *b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,*pmastsurf=NULL,
129 *cdnt=NULL,*cdnr=NULL,*cdni=NULL,*eme=NULL,alea=0.1,sum,
130 *pslavsurfold=NULL,*energyini=NULL,*energy=NULL,xn[3],e1[3],e2[3],
131 *smscale=NULL,*auw=NULL,*autloc=NULL,*xboun2=NULL,*coefmpc2=NULL;
132
133 FILE *f1;
134
135 /* dummy arguments for the results call */
136
137 double *veold=NULL,*accold=NULL,bet,gam,dtime,time=1.;
138
139 ITG *ipneigh=NULL,*neigh=NULL;
140
141 #ifdef SGI
142 ITG token;
143 #endif
144
145 irow=*irowp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
146 kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;
147
148 islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
149
150 if(*nmethod==13){
151 *nmethod=2;
152 igreen=1;
153 }
154
155 if(*nener==1){NNEW(ener,double,mi[0]**ne);}
156
157 for(k=0;k<3;k++){
158 strcpy1(&jobnamef[k*132],&jobnamec[k*132],132);
159 }
160
161 if(*mortar!=1){
162 maxprevcontel=*nslavs;
163 }else if(*mortar==1){
164 maxprevcontel=*nintpoint;
165 if(*nstate_!=0){
166 if(maxprevcontel!=0){
167 NNEW(islavsurfold,ITG,2**ifacecount+2);
168 NNEW(pslavsurfold,double,3**nintpoint);
169 memcpy(&islavsurfold[0],&islavsurf[0],
170 sizeof(ITG)*(2**ifacecount+2));
171 memcpy(&pslavsurfold[0],&pslavsurf[0],
172 sizeof(double)*(3**nintpoint));
173 }
174 }
175 }
176 nslavs_prev_step=*nslavs;
177
178 mt=mi[1]+1;
179 pi=4.*atan(1.);
180 constant=180./pi;
181
182 /* copying the frequency parameters */
183
184 if(igreen==0){
185 nev=mei[0];
186 ncv=mei[1];
187 mxiter=mei[2];
188 tol=fei[0];
189 fmin=2*pi*fei[1];
190 fmax=2*pi*fei[2];
191 }else{
192
193 /* Green functions: number of "eigenmodes" is the number of
194 unit forces */
195
196 nev=*nforc;
197 ncv=*nforc;
198 }
199
200 /* assigning the body forces to the elements */
201
202 if(*nbody>0){
203 /* ifreebody=*ne+1;
204 NNEW(ipobody,ITG,2*ifreebody**nbody);
205 for(k=1;k<=*nbody;k++){
206 FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
207 iendset,ialset,&inewton,nset,&ifreebody,&k));
208 RENEW(ipobody,ITG,2*(*ne+ifreebody));
209 }
210 RENEW(ipobody,ITG,2*(ifreebody-1));*/
211 if(*inewton==1){
212 printf("*ERROR in arpackcs: generalized gravity loading is not allowed in frequency calculations");
213 FORTRAN(stop,());
214 }
215 }
216
217 ne0=*ne;nkon0=*nkon;
218
219 /* contact conditions */
220
221 if(*iperturb!=0){
222
223 memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
224 maxlenmpc=mpcinfo[3];
225
226 inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
227 lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
228 &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
229 mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
230 iperturb,ikboun,nboun,co,istep,&xnoels);
231
232 if(ncont!=0){
233
234 NNEW(cg,double,3*ncont);
235 NNEW(straight,double,16*ncont);
236
237 /* 11 instead of 10: last position is reserved for the
238 local contact spring element number; needed as
239 pointer into springarea */
240
241 if(*mortar==0){
242 RENEW(kon,ITG,*nkon+11**nslavs);
243 NNEW(springarea,double,2**nslavs);
244 if(*nener==1){
245 RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
246 }
247 RENEW(ipkon,ITG,*ne+*nslavs);
248 RENEW(lakon,char,8*(*ne+*nslavs));
249
250 if(*norien>0){
251 RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
252 for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
253 }
254
255 RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
256 for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
257
258 if((maxprevcontel==0)&&(*nslavs!=0)){
259 RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
260 for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
261 xstate[k]=0.;
262 }
263 }
264 maxprevcontel=*nslavs;
265
266 NNEW(areaslav,double,*ifacecount);
267 NNEW(xmastnor,double,3*nmastnode[*ntie]);
268 }else if(*mortar==1){
269 NNEW(islavact,ITG,nslavnode[*ntie]);
270 DMEMSET(islavact,0,nslavnode[*ntie],1);
271 if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
272 NNEW(xmastnor,double,3*nmastnode[*ntie]);
273
274 *nintpoint=0;
275
276 precontact(&ncont,ntie,tieset,nset,set,istartset,
277 iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
278 cg,straight,co,vold,istep,&iinc,&iit,itiefac,
279 islavsurf,islavnode,imastnode,nslavnode,nmastnode,
280 imastop,mi,ipe,ime,tietol,&iflagact,
281 nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
282 nslavs);
283
284 /* changing the dimension of element-related fields */
285
286 RENEW(kon,ITG,*nkon+22**nintpoint);
287 RENEW(springarea,double,2**nintpoint);
288 RENEW(pmastsurf,double,6**nintpoint);
289
290 if(*nener==1){
291 RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
292 }
293 RENEW(ipkon,ITG,*ne+*nintpoint);
294 RENEW(lakon,char,8*(*ne+*nintpoint));
295
296 if(*norien>0){
297 RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
298 for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
299 }
300 RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
301 for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
302
303 /* interpolating the state variables */
304
305 if(*nstate_!=0){
306 if(maxprevcontel!=0){
307 RENEW(xstateini,double,
308 *nstate_*mi[0]*(ne0+maxprevcontel));
309 for(k=*nstate_*mi[0]*ne0;
310 k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
311 xstateini[k]=xstate[k];
312 }
313 }
314
315 RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
316 for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
317 xstate[k]=0.;
318 }
319
320 if((*nintpoint>0)&&(maxprevcontel>0)){
321
322 /* interpolation of xstate */
323
324 FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
325 &ne0,mi,xstate,pslavsurf,nstate_,
326 xstateini,islavsurf,islavsurfold,
327 pslavsurfold,tieset,ntie,itiefac));
328
329 }
330
331 if(maxprevcontel!=0){
332 SFREE(islavsurfold);SFREE(pslavsurfold);
333 }
334
335 maxprevcontel=*nintpoint;
336
337 RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
338 for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
339 xstateini[k]=xstate[k];
340 }
341 }
342 }
343
344 /* generating contact spring elements */
345
346 contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
347 ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
348 co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
349 &ne0,vini,nmethod,
350 iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
351 itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
352 imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
353 xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta,
354 xstateini,xstate,nstate_,&icutb,&ialeatoric,jobnamef,
355 &alea,auw,jqw,iroww,&nzsw);
356
357 printf("number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
358
359 /* determining the structure of the stiffness/mass matrix */
360
361 remastructar(ipompc,&coefmpc,&nodempc,nmpc,
362 &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
363 labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
364 kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
365 neq,nzs,nmethod,ithermal,iperturb,&mass,mi,ics,cs,
366 mcs,mortar,typeboun,&iit,&network,iexpl);
367 }
368 }
369
370 /* field for initial values of state variables (needed if
371 previous static step was viscoplastic and for contact */
372
373 if((*nstate_!=0)&&((*mortar==0)||(ncont==0))){
374 NNEW(xstateini,double,*nstate_*mi[0]*(ne0+*nslavs));
375 for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
376 xstateini[k]=xstate[k];
377 }
378 }
379
380 /* determining the internal forces and the stiffness coefficients */
381
382 NNEW(f,double,neq[1]);
383
384 /* allocating a field for the stiffness matrix */
385
386 NNEW(xstiff,double,(long long)27*mi[0]**ne);
387
388 iout=-1;
389 NNEW(v,double,mt**nk);
390 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
391 NNEW(fn,double,mt**nk);
392 NNEW(stx,double,6*mi[0]**ne);
393 NNEW(eme,double,6*mi[0]**ne);
394 NNEW(inum,ITG,*nk);
395 if(*iperturb==0){
396 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
397 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
398 ielorien,norien,orab,ntmat_,t0,t0,ithermal,
399 prestr,iprestr,filab,eme,emn,een,iperturb,
400 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
401 ndirboun,xboun,nboun,ipompc,
402 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
403 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
404 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
405 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
406 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
407 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
408 fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
409 &reltime,&ne0,thicke,shcon,nshcon,
410 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
411 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
412 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
413 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
414 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
415 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
416 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
417 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
418 &intscheme);
419 }else{
420 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
421 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
422 ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
423 prestr,iprestr,filab,eme,emn,een,iperturb,
424 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
425 ndirboun,xboun,nboun,ipompc,
426 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
427 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
428 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
429 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
430 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
431 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
432 fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
433 &reltime,&ne0,thicke,shcon,nshcon,
434 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
435 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
436 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
437 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
438 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
439 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
440 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
441 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
442 &intscheme);
443 }
444 SFREE(f);SFREE(v);SFREE(fn);SFREE(stx);SFREE(eme);SFREE(inum);
445 iout=1;
446
447 /* for the frequency analysis linear strain and elastic properties
448 are used */
449
450 iperturb[1]=0;ielas=1;
451
452 /* determining the maximum number of sectors to be plotted */
453
454 ngraph=1;
455 for(j=0;j<*mcs;j++){
456 if(cs[17*j+4]>ngraph) ngraph=cs[17*j+4];
457 }
458
459 /* assigning nodes and elements to sectors */
460
461 NNEW(inocs,ITG,*nk);
462 NNEW(ielcs,ITG,*ne);
463 ielset=cs[12];
464 if((*mcs!=1)||(ielset!=0)){
465 for(i=0;i<*nk;i++) inocs[i]=-1;
466 for(i=0;i<*ne;i++) ielcs[i]=-1;
467 }
468
469 for(i=0;i<*mcs;i++){
470 is=cs[17*i+4];
471 if((is==1)&&(*mcs==1)) continue;
472 ielset=cs[17*i+12];
473 if(ielset==0) continue;
474 for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
475 if(ialset[i1]>0){
476 iel=ialset[i1]-1;
477 if(ipkon[iel]<0) continue;
478 ielcs[iel]=i;
479 indexe=ipkon[iel];
480 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
481 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
482 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
483 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
484 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
485 else if (strcmp1(&lakon[8*iel+3],"6")==0)nope=6;
486 else if (strcmp1(&lakon[8*iel],"ES")==0){
487 lakonl[0]=lakon[8*iel+7];
488 nope=atoi(lakonl)+1;}
489 else continue;
490
491 for(i2=0;i2<nope;++i2){
492 node=kon[indexe+i2]-1;
493 inocs[node]=i;
494 }
495 }
496 else{
497 iel=ialset[i1-2]-1;
498 do{
499 iel=iel-ialset[i1];
500 if(iel>=ialset[i1-1]-1) break;
501 if(ipkon[iel]<0) continue;
502 ielcs[iel]=i;
503 indexe=ipkon[iel];
504 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
505 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
506 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
507 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
508 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
509 else {nope=6;}
510 for(i2=0;i2<nope;++i2){
511 node=kon[indexe+i2]-1;
512 inocs[node]=i;
513 }
514 }while(1);
515 }
516 }
517 }
518
519 /* loop over the nodal diameters */
520
521 for(nm=cs[1];nm<=cs[2];++nm){
522
523 NNEW(ad,double,neq[1]);
524 NNEW(au,double,nzs[1]);
525
526 NNEW(adb,double,neq[1]);
527 NNEW(aub,double,nzs[1]);
528
529 NNEW(fext,double,neq[1]);
530 if(*iperturb==0){
531 FORTRAN(mafillsmcs,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun
532 ,nboun,
533 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
534 nforc,nelemload,sideload,xload,nload,xbody,ipobody,
535 nbody,cgr,
536 ad,au,fext,nactdof,icol,jq,irow,&neq[1],nzl,nmethod,
537 ikmpc,ilmpc,ikboun,ilboun,
538 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
539 ielorien,norien,orab,ntmat_,
540 t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
541 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
542 xstiff,npmat_,&dtime,matname,mi,
543 ics,cs,&nm,ncmat_,labmpc,&mass,&stiffness,&buckling,
544 &rhsi,
545 &intscheme,mcs,&coriolis,ibody,xloadold,&reltime,
546 ielcs,veold,
547 springarea,thicke,integerglob,doubleglob,
548 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
549 pmastsurf,mortar,clearini,ielprop,prop,&ne0,&kscale,
550 xstateini,xstate,nstate_,set,nset));
551 }
552 else{
553 FORTRAN(mafillsmcs,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,
554 nboun,
555 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
556 nforc,nelemload,sideload,xload,nload,xbody,ipobody,
557 nbody,cgr,
558 ad,au,fext,nactdof,icol,jq,irow,&neq[1],nzl,nmethod,
559 ikmpc,ilmpc,ikboun,ilboun,
560 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
561 ielorien,norien,orab,ntmat_,
562 t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
563 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
564 xstiff,npmat_,&dtime,matname,mi,
565 ics,cs,&nm,ncmat_,labmpc,&mass,&stiffness,&buckling,
566 &rhsi,
567 &intscheme,mcs,&coriolis,ibody,xloadold,&reltime,
568 ielcs,veold,
569 springarea,thicke,integerglob,doubleglob,
570 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
571 pmastsurf,mortar,clearini,ielprop,prop,&ne0,&kscale,
572 xstateini,xstate,nstate_,set,nset));
573
574 if(nasym==1){
575 RENEW(au,double,nzs[2]+nzs[1]);
576 RENEW(aub,double,nzs[2]+nzs[1]);
577 symmetryflag=2;
578 inputformat=1;
579
580 FORTRAN(mafillsmcsas,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,
581 nboun,ipompc,nodempc,coefmpc,nmpc,nodeforc,
582 ndirforc,xforc,nforc,nelemload,sideload,xload,
583 nload,xbody,ipobody,nbody,cgr,ad,au,fext,nactdof,
584 icol,jq,irow,&neq[1],nzl,nmethod,ikmpc,ilmpc,
585 ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,alcon,
586 nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
587 t0,t1,ithermal,prestr,iprestr,vold,iperturb,sti,
588 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,
589 nplkcon,xstiff,npmat_,&dtime,matname,mi,ics,cs,
590 &nm,ncmat_,labmpc,&mass,&stiffness,&buckling,
591 &rhsi,&intscheme,mcs,&coriolis,ibody,xloadold,
592 &reltime,ielcs,veold,springarea,thicke,
593 integerglob,doubleglob,tieset,istartset,iendset,
594 ialset,ntie,&nasym,nstate_,xstateini,xstate,
595 pslavsurf,pmastsurf,mortar,clearini,ielprop,prop,
596 &ne0,&kscale,set,nset));
597
598 }
599 }
600
601 SFREE(fext);
602
603 if(*nmethod==0){
604
605 /* error occurred in mafill: storing the geometry in frd format */
606
607 ++*kode;
608 NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
609 if(strcmp1(&filab[1044],"ZZS")==0){
610 NNEW(neigh,ITG,40**ne);
611 NNEW(ipneigh,ITG,*nk);
612 }
613
614 frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
615 kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
616 nstate_,istep,&iinc,ithermal,qfn,&j,&nm,trab,inotr,
617 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
618 mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
619 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
620 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,
621 ielprop,prop,sti);
622
623 if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
624 SFREE(inum);FORTRAN(stop,());
625
626 }
627
628 /* LU decomposition of the left hand matrix */
629
630 if(igreen==0){
631 if(nasym==1){sigma=0.;}else{sigma=1.;}
632 }else{
633 if(*nforc>0){sigma=xforc[0];}
634 }
635
636 if(*isolver==0){
637 #ifdef SPOOLES
638 spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,&symmetryflag,
639 &inputformat,&nzs[2]);
640 #else
641 printf("*ERROR in arpackcs: the SPOOLES library is not linked\n\n");
642 FORTRAN(stop,());
643 #endif
644 }
645 else if(*isolver==4){
646 #ifdef SGI
647 token=1;
648 sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,token);
649 #else
650 printf("*ERROR in arpackcs: the SGI library is not linked\n\n");
651 FORTRAN(stop,());
652 #endif
653 }
654 else if(*isolver==5){
655 #ifdef TAUCS
656 tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],nzs);
657 #else
658 printf("*ERROR in arpackcs: the TAUCS library is not linked\n\n");
659 FORTRAN(stop,());
660 #endif
661 }
662 else if(*isolver==6){
663 #ifdef MATRIXSTORAGE
664 matrixstorage(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1],
665 ntrans,inotr,trab,co,nk,nactdof,jobnamec,mi,ipkon,
666 lakon,kon,ne,mei,nboun,nmpc,cs,mcs,ithermal,nmethod);
667 #else
668 printf("*ERROR in arpack: the MATRIXSTORAGE library is not linked\n\n");
669 FORTRAN(stop,());
670 #endif
671 }
672 else if(*isolver==7){
673 #ifdef PARDISO
674 pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,
675 &symmetryflag,&inputformat,jq,&nzs[2]);
676 #else
677 printf("*ERROR in arpack: the PARDISO library is not linked\n\n");
678 FORTRAN(stop,());
679 #endif
680 }
681 else if(*isolver==8){
682 #ifdef PASTIX
683 pastix_factor_main(ad,au,adb,aub,&sigma,icol,irow,&neq[1],nzs,
684 &symmetryflag,&inputformat,jq,&nzs[2]);
685 #else
686 printf("*ERROR in arpack: the PASTIX library is not linked\n\n");
687 FORTRAN(stop,());
688 #endif
689 }
690
691 // SFREE(au);SFREE(ad);
692
693 /* calculating the eigenvalues and eigenmodes */
694
695 if(igreen==0){
696
697 printf(" Calculating the eigenvalues and the eigenmodes\n");
698
699 ido=0;
700 ldz=neq[1];
701 for(k=0;k<11;k++)iparam[k]=0;
702 iparam[0]=1;
703 iparam[2]=mxiter;
704 iparam[3]=1;
705 iparam[6]=3;
706
707 info=0;
708
709 NNEW(resid,double,neq[1]);
710 NNEW(z,double,(long long)ncv*neq[1]);
711 NNEW(workd,double,3*neq[1]);
712
713 if(nasym==1){
714 lworkl=3*ncv*(2+ncv);
715 NNEW(workl,double,lworkl);
716 FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
717 workl,&lworkl,&info));
718 }else{
719 lworkl=ncv*(8+ncv);
720 NNEW(workl,double,lworkl);
721 FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
722 iparam,ipntr,workd,workl,&lworkl,&info));
723 }
724
725 NNEW(temp_array,double,neq[1]);
726
727 while((ido==-1)||(ido==1)||(ido==2)){
728 if(ido==-1){
729 if(nasym==1){
730 FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow,nzs));
731 }else{
732 FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,
733 jq,irow));
734 }
735 }
736
737 /* solve the linear equation system */
738
739 if((ido==-1)||(ido==1)){
740
741 if(ido==-1){
742 if(*isolver==0){
743 #ifdef SPOOLES
744 spooles_solve(temp_array,&neq[1]);
745 #endif
746 }
747 else if(*isolver==4){
748 #ifdef SGI
749 sgi_solve(temp_array,token);
750 #endif
751 }
752 else if(*isolver==5){
753 #ifdef TAUCS
754 tau_solve(temp_array,&neq[1]);
755 #endif
756 }
757 else if(*isolver==7){
758 #ifdef PARDISO
759 pardiso_solve(temp_array,&neq[1],&symmetryflag,&inputformat,&nrhs);
760 #endif
761 }
762 else if(*isolver==8){
763 #ifdef PASTIX
764 if(pastix_solve(temp_array,&neq[1],&symmetryflag,&nrhs)==-1)
765 printf("*WARNING in arpackcs: solving step did not converge! Continuing anyway!\n");
766 #endif
767 }
768 for(jrow=0;jrow<neq[1];jrow++){
769 workd[ipntr[1]-1+jrow]=temp_array[jrow];
770 }
771 }
772 else if(ido==1){
773 if(*isolver==0){
774 #ifdef SPOOLES
775 spooles_solve(&workd[ipntr[2]-1],&neq[1]);
776 #endif
777 }
778 else if(*isolver==4){
779 #ifdef SGI
780 sgi_solve(&workd[ipntr[2]-1],token);
781 #endif
782 }
783 else if(*isolver==5){
784 #ifdef TAUCS
785 tau_solve(&workd[ipntr[2]-1],&neq[1]);
786 #endif
787 }
788 else if(*isolver==7){
789 #ifdef PARDISO
790 pardiso_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag,&inputformat,&nrhs);
791 #endif
792 }
793 else if(*isolver==8){
794 #ifdef PASTIX
795 if(pastix_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag,&nrhs)==-1)
796 printf("*WARNING in arpackcs: solving step did not converge! Continuing anyway!\n");
797 #endif
798 }
799 for(jrow=0;jrow<neq[1];jrow++){
800 workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
801 }
802
803 }
804 }
805
806 if(ido==2){
807 if(nasym==1){
808 FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
809 adb,aub,jq,irow,nzs));
810 }else{
811 FORTRAN(op,(neq,&workd[ipntr[0]-1],&workd[ipntr[1]-1],
812 adb,aub,jq,irow));
813 }
814 }
815
816 if(nasym==1){
817 FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
818 iparam,ipntr,workd,workl,&lworkl,&info));
819 }else{
820 FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,
821 z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
822 }
823 }
824
825 if(info!=0){
826 printf("*ERROR in d[n,s]aupd: info=%" ITGFORMAT "\n",info);
827 printf(" # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
828 }
829
830 NNEW(select,ITG,ncv);
831
832 if(nasym==1){
833 NNEW(d,double,nev+1);
834 NNEW(di,double,nev+1);
835 NNEW(workev,double,3*ncv);
836 FORTRAN(dneupd,(&rvec,howmny,select,d,di,z,&ldz,&sigma,&sigmai,
837 workev,bmat,&neq[1],which,&nev,&tol,resid,
838 &ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
839 SFREE(workev);
840 NNEW(dc,double,2*nev);
841 NNEW(nmc,ITG,nev);
842
843 /* storing as complex number and taking the square root */
844
845 for(j=0;j<nev;j++){
846 dc[2*j]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])+d[j])/sqrt(2.);
847 dc[2*j+1]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])-d[j])/sqrt(2.);
848 if(di[j]<0.) dc[2*j+1]=-dc[2*j+1];
849 nmc[j]=nm;
850 }
851 FORTRAN(writeevcscomplex,(dc,&nev,nmc,&fmin,&fmax));
852 SFREE(di);SFREE(dc);SFREE(nmc);
853 }else{
854 NNEW(d,double,nev);
855 FORTRAN(dseupd,(&rvec,howmny,select,d,z,&ldz,&sigma,bmat,&neq[1],which,&nev,
856 &tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
857 FORTRAN(writeevcs,(d,&nev,&nm,&fmin,&fmax));
858 }
859 SFREE(select);SFREE(resid);SFREE(workd);SFREE(workl);
860
861 if(info!=0){
862 printf("*ERROR in d[n,s]eupd: info=%" ITGFORMAT "\n",info);
863 }
864
865 /* check the normalization of the eigenmodes */
866
867 for(j=0;j<nev;j++){
868 kref=j*neq[1];
869 if(nasym==1){
870 FORTRAN(opas,(&neq[1],&z[kref],temp_array,
871 adb,aub,jq,irow,nzs));
872 }else{
873 FORTRAN(op,(neq,&z[kref],temp_array,
874 adb,aub,jq,irow));
875 }
876 sum=0;
877 for(k=0;k<neq[1];k++){sum+=z[kref+k]*temp_array[k];}
878 printf("U^T*M*U=%f for eigenmode %d\n",sum,j+1);
879
880 /* normalizing the eigenmode (if not yet normalized by ARPACK */
881
882 if(fabs(sum-1.)>1.e-5){
883 printf("normalizing the eigenmode \n");
884 for(k=0;k<neq[1];k++){z[kref+k]/=sqrt(sum);}
885 }
886
887 }
888
889 SFREE(temp_array);
890
891 /* for double eigenmodes: the eigenmode for which the largest
892 amplitude has the lowest dof comes first */
893
894 neq2=neq[1]/2;
895 for (j=0;j<nev;j+=2){
896
897 ampmax=0.;kmax1=0;
898 for(k=0;k<neq2;k++){
899 amp=z[2*j*neq[1]+k]*z[2*j*neq[1]+k]+
900 z[2*j*neq[1]+neq2+k]*z[2*j*neq[1]+neq2+k];
901 if(amp>ampmax){ampmax=amp;kmax1=k;}
902 }
903
904 ampmax=0.;kmax2=0;
905 for(k=0;k<neq2;k++){
906 amp=z[(2*j+1)*neq[1]+k]*z[(2*j+1)*neq[1]+k]+
907 z[(2*j+1)*neq[1]+neq2+k]*z[(2*j+1)*neq[1]+neq2+k];
908 if(amp>ampmax){ampmax=amp;kmax2=k;}
909 }
910
911 if(kmax2<kmax1){
912 printf("exchange!\n");
913 NNEW(zstorage,double,neq[1]);
914 memcpy(zstorage,&z[2*j*neq[1]],sizeof(double)*neq[1]);
915 memcpy(&z[2*j*neq[1]],&z[(2*j+1)*neq[1]],sizeof(double)*neq[1]);
916 memcpy(&z[(2*j+1)*neq[1]],zstorage,sizeof(double)*neq[1]);
917 SFREE(zstorage);
918 }
919 }
920 }else{
921
922 /* Green functions */
923
924 /* storing omega_0^2 into d */
925
926 NNEW(d,double,*nforc);
927 for(j=0;j<*nforc;j++){d[j]=xforc[0];}
928 NNEW(z,double,neq[1]);
929
930 }
931
932 /* writing the eigenvalues and mass matrix to a binary file */
933
934 if(mei[3]==1){
935
936 strcpy(fneig,jobnamec);
937 strcat(fneig,".eig");
938
939 /* the first time the file is erased before writing, all subsequent
940 times the data is appended */
941
942 if(*nevtot==0){
943 if((f1=fopen(fneig,"wb"))==NULL){
944 printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
945 exit(0);
946 }
947
948 /* storing a one as indication that this was a
949 cyclic symmetry calculation */
950
951 if(fwrite(&one,sizeof(ITG),1,f1)!=1){
952 printf("*ERROR saving the cyclic symmetry flag to the eigenvalue file...");
953 exit(0);
954 }
955
956 /* Hermitian */
957
958 if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
959 printf("*ERROR saving the Hermitian flag to the eigenvalue file...");
960 exit(0);
961 }
962
963 /* perturbation parameter iperturb[0] */
964
965 if(fwrite(&iperturb[0],sizeof(ITG),1,f1)!=1){
966 printf("*ERROR saving the perturbation flag to the eigenvalue file...");
967 exit(0);
968 }
969
970 /* reference displacements */
971
972 if(iperturb[0]==1){
973 if(fwrite(vold,sizeof(double),mt**nk,f1)!=mt**nk){
974 printf("*ERROR saving the reference displacements to the eigenvalue file...");
975 exit(0);
976 }
977 }
978
979 }else{
980 if((f1=fopen(fneig,"ab"))==NULL){
981 printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
982 exit(0);
983 }
984 }
985
986 /* nodal diameter */
987
988 if(fwrite(&nm,sizeof(ITG),1,f1)!=1){
989 printf("*ERROR saving the nodal diameter to the eigenvalue file...");
990 exit(0);
991 }
992
993 /* number of eigenfrequencies */
994
995 if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
996 printf("*ERROR saving the number of eigenfrequencies to the eigenvalue file...");
997 exit(0);
998 }
999
1000 /* the eigenfrequencies are stored as radians/time */
1001
1002 if(fwrite(d,sizeof(double),nev,f1)!=nev){
1003 printf("*ERROR saving the eigenfrequencies to the eigenvalue file...");
1004 exit(0);
1005 }
1006 if(*nevtot==0){
1007
1008 /* stiffness matrix */
1009
1010 if(fwrite(ad,sizeof(double),neq[1],f1)!=neq[1]){
1011 printf("*ERROR saving the diagonal of the stiffness matrix to the eigenvalue file...");
1012 exit(0);
1013 }
1014 if(fwrite(au,sizeof(double),nzs[1],f1)!=nzs[1]){
1015 printf("*ERROR saving the off-diagonal terms of the stiffness matrix to the eigenvalue file...");
1016 exit(0);
1017 }
1018
1019 /* mass matrix */
1020
1021 if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
1022 printf("*ERROR saving the diagonal of the mass matrix to the eigenvalue file...");
1023 exit(0);
1024 }
1025 if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
1026 printf("*ERROR saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
1027 exit(0);
1028 }
1029 }
1030
1031 // SFREE(ad);SFREE(au);
1032 }
1033
1034 if(igreen==0){
1035
1036 /* calculating the participation factors and the relative effective
1037 modal mass */
1038
1039 if(*ithermal!=2){
1040 FORTRAN(effectivemodalmass,(neq,nactdof,mi,adb,aub,jq,irow,&nev,z,co,nk));
1041 }
1042 }
1043
1044 /* calculating the displacements and the stresses and storing */
1045 /* the results in frd format for each valid eigenmode */
1046
1047 /* for energy calculations in other sectors the stress and the
1048 mechanical strain have to be calculated */
1049
1050 if((ngraph>1)&&(strcmp1(&filab[522],"ENER")==0)){
1051 strcpy1(&filabcp[0],&filab[174],4);
1052 strcpy1(&filabcp[4],&filab[2697],4);
1053 strcpy1(&filab[174],"S ",4);
1054 strcpy1(&filab[2697],"ME ",4);
1055 }
1056
1057
1058 NNEW(v,double,2*mt**nk);
1059 NNEW(fn,double,2*mt**nk);
1060 NNEW(inum,ITG,*nk);
1061 NNEW(stx,double,2*6*mi[0]**ne);
1062 NNEW(eme,double,2*6*mi[0]**ne);
1063
1064 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1653],"MAXS")==0)||
1065 (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
1066 (strcmp1(&filab[1044],"ERR ")==0))
1067 NNEW(stn,double,12**nk);
1068
1069 if((strcmp1(&filab[261],"E ")==0)||(strcmp1(&filab[2523],"MAXE")==0))
1070 NNEW(een,double,12**nk);
1071 if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,2**nk);
1072 if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,12**nk);
1073 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1074 &&(*mortar==1)) NNEW(cdn,double,12**nk);
1075
1076 NNEW(temp_array,double,neq[1]);
1077 NNEW(coefmpcnew,double,*mpcend);
1078
1079 /* creating total fields for ngraph segments */
1080
1081 NNEW(cot,double,3**nk*ngraph);
1082 if(*ntrans>0){NNEW(inotrt,ITG,2**nk*ngraph);}
1083 if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0))
1084
1085 // real and imaginary part of the displacements
1086
1087 NNEW(vt,double,2*mt**nk*ngraph);
1088 if(strcmp1(&filab[87],"NT ")==0)
1089 NNEW(t1t,double,*nk*ngraph);
1090 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
1091 (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0))
1092
1093 // real and imaginary part of the stresses
1094
1095 NNEW(stnt,double,2*6**nk*ngraph);
1096 if(strcmp1(&filab[261],"E ")==0) NNEW(eent,double,2*6**nk*ngraph);
1097 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0))
1098
1099 // real and imaginary part of the forces
1100
1101 NNEW(fnt,double,2*mt**nk*ngraph);
1102 if(strcmp1(&filab[2697],"ME ")==0) NNEW(emnt,double,2*6**nk*ngraph);
1103 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1104 &&(*mortar==1)) NNEW(cdnt,double,2*6**nk*ngraph);
1105 if(strcmp1(&filab[522],"ENER")==0)
1106
1107 // real and imaginary part of the internal energy
1108
1109 NNEW(enernt,double,*nk*ngraph);
1110
1111 // stresses at the integration points for the error estimator and contact conditions
1112
1113 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
1114 ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0)))
1115 NNEW(stxt,double,2*6*mi[0]**ne*ngraph);
1116
1117 NNEW(kont,ITG,*nkon*ngraph);
1118 NNEW(ipkont,ITG,*ne*ngraph);
1119 for(l=0;l<*ne*ngraph;l++)ipkont[l]=-1;
1120 NNEW(lakont,char,8**ne*ngraph);
1121 NNEW(inumt,ITG,*nk*ngraph);
1122 NNEW(ielmatt,ITG,mi[2]**ne*ngraph);
1123
1124 nkt=ngraph**nk;
1125 net=ngraph**ne;
1126
1127 /* copying the coordinates of the first sector */
1128
1129 for(l=0;l<3**nk;l++){cot[l]=co[l];}
1130 if(*ntrans>0){for(l=0;l<*nk;l++){inotrt[2*l]=inotr[2*l];}}
1131 for(l=0;l<*nkon;l++){kont[l]=kon[l];}
1132 for(l=0;l<*ne;l++){ipkont[l]=ipkon[l];}
1133 for(l=0;l<8**ne;l++){lakont[l]=lakon[l];}
1134 for(l=0;l<*ne;l++){ielmatt[mi[2]*l]=ielmat[mi[2]*l];}
1135
1136 /* generating the coordinates for the other sectors */
1137
1138 icntrl=1;
1139
1140 FORTRAN(rectcyl,(cot,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1141
1142 for(jj=0;jj<*mcs;jj++){
1143 is=cs[17*jj+4];
1144 for(i=1;i<is;i++){
1145
1146 theta=i*2.*pi/cs[17*jj];
1147
1148 for(l=0;l<*nk;l++){
1149 if(inocs[l]==jj){
1150 cot[3*l+i*3**nk]=cot[3*l];
1151 cot[1+3*l+i*3**nk]=cot[1+3*l]+theta;
1152 cot[2+3*l+i*3**nk]=cot[2+3*l];
1153 if(*ntrans>0){inotrt[2*l+i*2**nk]=inotrt[2*l];}
1154 }
1155 }
1156 for(l=0;l<*nkon;l++){kont[l+i**nkon]=kon[l]+i**nk;}
1157 for(l=0;l<*ne;l++){
1158 if(ielcs[l]==jj){
1159 if(ipkon[l]>=0){
1160 ipkont[l+i**ne]=ipkon[l]+i**nkon;
1161 ielmatt[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
1162 for(l1=0;l1<8;l1++){
1163 l2=8*l+l1;
1164 lakont[l2+i*8**ne]=lakon[l2];
1165 }
1166 }
1167 }
1168 }
1169 }
1170 }
1171
1172 icntrl=-1;
1173
1174 FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
1175 &imag,mi,emnt));
1176
1177 /* check that the tensor fields which are extrapolated from the
1178 integration points are requested in global coordinates */
1179
1180 if(strcmp1(&filab[174],"S ")==0){
1181 if((strcmp1(&filab[179],"L")==0)&&(*norien>0)){
1182 printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculations\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1183 strcpy1(&filab[179],"G",1);
1184 }
1185 }
1186
1187 if(strcmp1(&filab[261],"E ")==0){
1188 if((strcmp1(&filab[266],"L")==0)&&(*norien>0)){
1189 printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1190 strcpy1(&filab[266],"G",1);
1191 }
1192 }
1193
1194 if(strcmp1(&filab[1479],"PHS ")==0){
1195 if((strcmp1(&filab[1484],"L")==0)&&(*norien>0)){
1196 printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1197 strcpy1(&filab[1484],"G",1);
1198 }
1199 }
1200
1201 if(strcmp1(&filab[1653],"MAXS")==0){
1202 if((strcmp1(&filab[1658],"L")==0)&&(*norien>0)){
1203 printf("\n*WARNING in arpackcs: element fields in cyclic symmetry calculation\n cannot be requested in local orientations;\n the global orientation will be used \n\n");
1204 strcpy1(&filab[1658],"G",1);
1205 }
1206 }
1207
1208 if(strcmp1(&filab[2523],"MAXE")==0){
1209 if((strcmp1(&filab[2528],"L")==0)&&(*norien>0)){
1210 printf("\n*WARNING in arpackcs: 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[1658],"G",1);
1212 }
1213 }
1214
1215 /* allocating fields for magnitude and phase information of
1216 displacements and stresses */
1217
1218 if(strcmp1(&filab[870],"PU")==0){
1219 NNEW(vr,double,mt*nkt);
1220 NNEW(vi,double,mt*nkt);
1221 }
1222
1223 if(strcmp1(&filab[1479],"PHS")==0){
1224 NNEW(stnr,double,6*nkt);
1225 NNEW(stni,double,6*nkt);
1226 }
1227
1228 if(strcmp1(&filab[2610],"PRF")==0){
1229 NNEW(fnr,double,mt*nkt);
1230 NNEW(fni,double,mt*nkt);
1231 }
1232
1233 if(strcmp1(&filab[3915],"PCON")==0){
1234 NNEW(cdnr,double,6*nkt);
1235 NNEW(cdni,double,6*nkt);
1236 }
1237
1238 if(strcmp1(&filab[1566],"MAXU")==0){
1239 NNEW(vmax,double,4*nkt);
1240 }
1241
1242 if(strcmp1(&filab[1653],"MAXS")==0){
1243 NNEW(stnmax,double,7*nkt);
1244 }
1245
1246 if(strcmp1(&filab[2523],"MAXE")==0){
1247 NNEW(eenmax,double,7*nkt);
1248 }
1249
1250 /* determining a local coordinate system based on the rotation axis */
1251
1252 FORTRAN(localaxescs,(cs,mcs,e1,e2,xn));
1253 NNEW(turdir,char,nev);
1254
1255 /* start of output calculations */
1256
1257 lfin=0;
1258 for(j=0;j<nev;++j){
1259 lint=lfin;
1260 if(igreen==0){
1261 lfin=lfin+neq[1];
1262 }else{
1263
1264 /* Green function: solve the system ([K]-w**2*[M]){X}={Y} */
1265
1266 node=nodeforc[2*j];
1267 idir=ndirforc[j];
1268
1269 /* check whether degree of freedom is active */
1270
1271 if(nactdof[mt*(node-1)+idir]==0){
1272 printf("*ERROR in linstatic: degree of freedom corresponding to node %d \n and direction %d is not active: no unit force can be applied\n",node,idir);
1273 FORTRAN(stop,());
1274 }
1275
1276 /* defining a unit force on the rhs */
1277
1278 DMEMSET(z,0,neq[1],0.);
1279 z[nactdof[mt*(node-1)+idir]-1]=1.;
1280
1281 /* solving the system */
1282
1283 if(neq[1]>0){
1284 if(*isolver==0){
1285 #ifdef SPOOLES
1286 spooles_solve(z,&neq[1]);
1287 #endif
1288 }
1289 else if(*isolver==7){
1290 #ifdef PARDISO
1291 pardiso_solve(z,&neq[1],&symmetryflag,&inputformat,&nrhs);
1292 #endif
1293
1294 }
1295 else if(*isolver==8){
1296 #ifdef PASTIX
1297 if(pastix_solve(z,&neq[1],&symmetryflag,&nrhs)==-1)
1298 printf("*WARNING in arpackcs: solving step did not converge! Continuing anyway!\n");
1299 #endif
1300
1301 }
1302 }
1303 }
1304
1305 if(mei[3]==1){
1306 if(fwrite(&z[lint],sizeof(double),neq[1],f1)!=neq[1]){
1307 printf("*ERROR saving data to the eigenvalue file...");
1308 exit(0);
1309 }
1310 }
1311
1312 if(igreen==0){
1313
1314 /* check whether the frequency belongs to the requested
1315 interval */
1316
1317 if(fmin>-0.5){
1318 if(fmin*fmin>d[j]) continue;
1319 }
1320 if(fmax>-0.5){
1321 if(fmax*fmax<d[j]) continue;
1322 }
1323 }
1324
1325 if(*nprint>0)FORTRAN(writehe,(&j));
1326
1327 NNEW(eei,double,6*mi[0]*ne0);
1328 if(*nener==1){
1329 NNEW(stiini,double,6*mi[0]*ne0);
1330 NNEW(emeini,double,6*mi[0]*ne0);
1331 NNEW(enerini,double,mi[0]*ne0);}
1332
1333 DMEMSET(v,0,2*mt**nk,0.);
1334
1335 /* calculating the strains, stresses... (real and imaginary part) for
1336 one specific eigenvalue */
1337
1338 for(k=0;k<neq[1];k+=neq[1]/2){
1339
1340 if(k==0) {kk=0;kkv=0;kk6=0;kkx=0;if(*nprint>0)FORTRAN(writere,());}
1341 else {kk=*nk;kkv=mt**nk;kk6=6**nk;kkx=6*mi[0]**ne;
1342 if(*nprint>0)FORTRAN(writeim,());}
1343
1344 /* generating the cyclic MPC's (needed for nodal diameters
1345 different from 0 */
1346
1347 for(i=0;i<*nmpc;i++){
1348 index=ipompc[i]-1;
1349 /* check whether thermal mpc */
1350 if(nodempc[3*index+1]==0) continue;
1351 coefmpcnew[index]=coefmpc[index];
1352 while(1){
1353 index=nodempc[3*index+2];
1354 if(index==0) break;
1355 index--;
1356
1357 icomplex=0;
1358 inode=nodempc[3*index];
1359 if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
1360 icomplex=atoi(&labmpc[20*i+6]);}
1361 else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
1362 for(ij=0;ij<*mcs;ij++){
1363 lprev=cs[ij*17+13];
1364 ilength=cs[ij*17+3];
1365 FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
1366 if(id!=0){
1367 if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
1368 }
1369 }
1370 }
1371
1372 if(icomplex!=0){
1373 idir=nodempc[3*index+1];
1374 idof=nactdof[mt*(inode-1)+idir]-1;
1375 if(idof<=-1){xreal=1.;ximag=1.;}
1376 else{xreal=z[lint+idof];ximag=z[lint+idof+neq[1]/2];}
1377 if(k==0) {
1378 if(fabs(xreal)<1.e-30)xreal=1.e-30;
1379 coefmpcnew[index]=coefmpc[index]*
1380 (cs[17*(icomplex-1)+14]+ximag/xreal*cs[17*(icomplex-1)+15]);}
1381 else {
1382 if(fabs(ximag)<1.e-30)ximag=1.e-30;
1383 coefmpcnew[index]=coefmpc[index]*
1384 (cs[17*(icomplex-1)+14]-xreal/ximag*cs[17*(icomplex-1)+15]);}
1385 }
1386 else{coefmpcnew[index]=coefmpc[index];}
1387 }
1388 }
1389
1390 if(*iperturb==0){
1391 results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1392 &stx[kkx],elcon,
1393 nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1394 norien,orab,ntmat_,t0,t0,ithermal,
1395 prestr,iprestr,filab,&eme[kkx],&emn[kk6],&een[kk6],iperturb,
1396 f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
1397 nodeboun,ndirboun,xboun,nboun,ipompc,
1398 nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1399 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1400 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1401 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1402 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1403 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1404 nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1405 &ne0,thicke,shcon,nshcon,
1406 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1407 mortar,islavact,&cdn[kk6],islavnode,nslavnode,ntie,clearini,
1408 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
1409 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
1410 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
1411 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
1412 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
1413 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
1414 &intscheme);}
1415 else{
1416 results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
1417 &stx[kkx],elcon,
1418 nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1419 norien,orab,ntmat_,t0,t1old,ithermal,
1420 prestr,iprestr,filab,&eme[kkx],&emn[kk6],&een[kk6],iperturb,
1421 f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
1422 nodeboun,ndirboun,xboun,nboun,ipompc,
1423 nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1424 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1425 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
1426 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
1427 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1428 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1429 nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1430 &ne0,thicke,shcon,nshcon,
1431 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1432 mortar,islavact,&cdn[kk6],islavnode,nslavnode,ntie,clearini,
1433 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
1434 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
1435 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
1436 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
1437 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
1438 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
1439 &intscheme);
1440 }
1441
1442 }
1443 SFREE(eei);
1444 if(*nener==1){SFREE(stiini);SFREE(emeini);SFREE(enerini);}
1445
1446 /* determining the turning direction */
1447
1448 FORTRAN(turningdirection,(v,e1,e2,xn,mi,nk,&turdir[j],lakon,ipkon,
1449 kon,ne,co));
1450
1451 if(strcmp1(&filab[1566],"MAXU")==0){
1452
1453 /* determining the ray vector; the components of the
1454 ray vector are the coordinates of the node in node set
1455 RAY */
1456
1457 iray=0;
1458 for(i=0;i<*nset;i++){
1459 if(strcmp1(&set[81*i],"RAYN")==0){
1460 iray=ialset[istartset[i]-1];
1461 vray[0]=co[3*iray-3];
1462 vray[1]=co[3*iray-2];
1463 vray[2]=co[3*iray-1];
1464 break;
1465 }
1466 }
1467 if(iray==0){
1468 printf("/n*ERROR in arpackcs: no light ray vector/n/n");
1469 FORTRAN(stop,());
1470 }
1471
1472 /* initialization */
1473
1474 for(l1=0;l1<4**nk;l1++){vmax[l1]=0.;}
1475
1476 /* vector p1 is a point on the rotation axis
1477 vector p2 is a unit vector along the axis */
1478
1479 for(l2=0;l2<3;l2++){p1[l2]=cs[5+l2];}
1480 for(l2=0;l2<3;l2++){p2[l2]=cs[8+l2]-p1[l2];}
1481 dd=sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]);
1482 for(l2=0;l2<3;l2++){p2[l2]/=dd;}
1483
1484 /* determine the time for the worst displacement
1485 orthogonal to a give light ray vector ; */
1486
1487 for(l1=0;l1<*nk;l1++){
1488
1489 /* determining a vector through node (l1+1) and
1490 orthogonal to the rotation axis */
1491
1492 for(l2=0;l2<3;l2++){q[l2]=co[3*l1+l2]-p1[l2];}
1493 dd=q[0]*p2[0]+q[1]*p2[1]+q[2]*p2[2];
1494 for(l2=0;l2<3;l2++){q[l2]-=dd*p2[l2];}
1495
1496 /* determining a vector tan orthogonal to vector q
1497 and the ray vector */
1498
1499 tan[0]=q[1]*vray[2]-q[2]*vray[1];
1500 tan[1]=q[2]*vray[0]-q[0]*vray[2];
1501 tan[2]=q[0]*vray[1]-q[1]*vray[0];
1502
1503 printf("tangent= %" ITGFORMAT ",%e,%e,%e\n",l1,tan[0],tan[1],tan[2]);
1504
1505 worstumax=0.;
1506 iworsttime=0;
1507 for(l3=0;l3<360;l3++){
1508 ctl=cos(l3/constant);
1509 stl=sin(l3/constant);
1510 for(l2=1;l2<4;l2++){
1511 l=mt*l1+l2;
1512 vl[l2]=ctl*v[l]-stl*v[l+mt**nk];
1513 }
1514
1515 /* displacement component along the tangent vector
1516 (no absolute value!) */
1517
1518 dd=vl[1]*tan[0]+vl[2]*tan[1]+vl[3]*tan[2];
1519 if(dd>worstumax){
1520 worstumax=dd;
1521 iworsttime=l3;
1522 }
1523 }
1524 ctl=cos(iworsttime/constant);
1525 stl=sin(iworsttime/constant);
1526 for(l2=1;l2<4;l2++){
1527 l=mt*l1+l2;
1528 vl[l2]=ctl*v[l]-stl*v[l+mt**nk];
1529 }
1530 vmax[4*l1]=1.*iworsttime;
1531 vmax[4*l1+1]=vl[1];
1532 vmax[4*l1+2]=vl[2];
1533 vmax[4*l1+3]=vl[3];
1534
1535 }
1536 }
1537
1538 /* determine the worst principal stress anywhere
1539 in the structure as a function of time;
1540 the worst principal stress is the maximum
1541 of the absolute value of the principal stresses */
1542
1543 if(strcmp1(&filab[1653],"MAXS")==0){
1544
1545 /* determining the set of nodes for the
1546 worst principal stress calculation */
1547
1548 ielset=0;
1549 for(i=0;i<*nset;i++){
1550 if(strcmp1(&set[81*i],"STRESSDOMAINN")==0){
1551 ielset=i+1;
1552 break;
1553 }
1554 }
1555 if(ielset==0){
1556 printf("\n*ERROR in arpackcs: no node set for MAXS\n");
1557 printf(" (must have the name STRESSDOMAIN)\n\n");
1558 FORTRAN(stop,());
1559 }
1560
1561 for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
1562 if(ialset[i1]>0){
1563 l1=ialset[i1]-1;
1564
1565 worstpsmax=0.;
1566 for(l3=0;l3<360;l3++){
1567 ctl=cos(l3/constant);
1568 stl=sin(l3/constant);
1569 for(l2=0;l2<6;l2++){
1570 l=6*l1+l2;
1571 stnl[l2]=ctl*stn[l]-stl*stn[l+6**nk];
1572 }
1573
1574 /* determining the eigenvalues */
1575
1576 v1=stnl[0]+stnl[1]+stnl[2];
1577 v2=stnl[1]*stnl[2]+stnl[0]*stnl[2]+stnl[0]*stnl[1]-
1578 (stnl[5]*stnl[5]+stnl[4]*stnl[4]+stnl[3]*stnl[3]);
1579 v3=stnl[0]*(stnl[1]*stnl[2]-stnl[5]*stnl[5])
1580 -stnl[3]*(stnl[3]*stnl[2]-stnl[4]*stnl[5])
1581 +stnl[4]*(stnl[3]*stnl[5]-stnl[4]*stnl[1]);
1582 bb=v2-v1*v1/3.;
1583 cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1584 if(fabs(bb)<=1.e-10){
1585 if(fabs(cc)>1.e-10){
1586 al[0]=-pow(cc,(1./3.));
1587 }else{
1588 al[0]=0.;
1589 }
1590 al[1]=al[0];
1591 al[2]=al[0];
1592 }else{
1593 cm=2.*sqrt(-bb/3.);
1594 cn=3.*cc/(cm*bb);
1595 if(fabs(cn)>1.){
1596 if(cn>1.){
1597 cn=1.;
1598 }else{
1599 cn=-1.;
1600 }
1601 }
1602 tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1603 al[0]=cm*cos(tt);
1604 al[1]=cm*cos(tt+2.*pi/3.);
1605 al[2]=cm*cos(tt+4.*pi/3.);
1606 }
1607 for(l2=0;l2<3;l2++){
1608 al[l2]+=v1/3.;
1609 }
1610 dd=fabs(al[0]);
1611 if(fabs(al[1])>dd) dd=fabs(al[1]);
1612 if(fabs(al[2])>dd) dd=fabs(al[2]);
1613 if(dd>worstpsmax){
1614 worstpsmax=dd;
1615 stnmax[7*l1]=dd;
1616 for(l2=1;l2<7;l2++){
1617 stnmax[7*l1+l2]=stnl[l2-1];
1618 }
1619 }
1620 }
1621
1622 }else{
1623 l1=ialset[i1-2]-1;
1624 do{
1625 l1=l1-ialset[i1];
1626 if(l1>=ialset[i1-1]-1) break;
1627
1628 worstpsmax=0.;
1629 for(l3=0;l3<360;l3++){
1630 ctl=cos(l3/constant);
1631 stl=sin(l3/constant);
1632 for(l2=0;l2<6;l2++){
1633 l=6*l1+l2;
1634 stnl[l2]=ctl*stn[l]-stl*stn[l+6**nk];
1635 }
1636
1637 /* determining the eigenvalues */
1638
1639 v1=stnl[0]+stnl[1]+stnl[2];
1640 v2=stnl[1]*stnl[2]+stnl[0]*stnl[2]+stnl[0]*stnl[1]-
1641 (stnl[5]*stnl[5]+stnl[4]*stnl[4]+stnl[3]*stnl[3]);
1642 v3=stnl[0]*(stnl[1]*stnl[2]-stnl[5]*stnl[5])
1643 -stnl[3]*(stnl[3]*stnl[2]-stnl[4]*stnl[5])
1644 +stnl[4]*(stnl[3]*stnl[5]-stnl[4]*stnl[1]);
1645 bb=v2-v1*v1/3.;
1646 cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1647 if(fabs(bb)<=1.e-10){
1648 if(fabs(cc)>1.e-10){
1649 al[0]=-pow(cc,(1./3.));
1650 }else{
1651 al[0]=0.;
1652 }
1653 al[1]=al[0];
1654 al[2]=al[0];
1655 }else{
1656 cm=2.*sqrt(-bb/3.);
1657 cn=3.*cc/(cm*bb);
1658 if(fabs(cn)>1.){
1659 if(cn>1.){
1660 cn=1.;
1661 }else{
1662 cn=-1.;
1663 }
1664 }
1665 tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1666 al[0]=cm*cos(tt);
1667 al[1]=cm*cos(tt+2.*pi/3.);
1668 al[2]=cm*cos(tt+4.*pi/3.);
1669 }
1670 for(l2=0;l2<3;l2++){
1671 al[l2]+=v1/3.;
1672 }
1673 dd=fabs(al[0]);
1674 if(fabs(al[1])>dd) dd=fabs(al[1]);
1675 if(fabs(al[2])>dd) dd=fabs(al[2]);
1676 if(dd>worstpsmax){
1677 worstpsmax=dd;
1678 stnmax[7*l1]=dd;
1679 for(l2=1;l2<7;l2++){
1680 stnmax[7*l1+l2]=stnl[l2-1];
1681 }
1682 }
1683 }
1684
1685 }while(1);
1686 }
1687 }
1688 }
1689
1690 /* determine the worst principal strain anywhere
1691 in the structure as a function of time;
1692 the worst principal strain is the maximum
1693 of the absolute value of the principal strains,
1694 times its original sign */
1695
1696 if(strcmp1(&filab[2523],"MAXE")==0){
1697
1698 /* determining the set of nodes for the
1699 worst principal strain calculation */
1700
1701 ielset=0;
1702 for(i=0;i<*nset;i++){
1703 if(strcmp1(&set[81*i],"STRAINDOMAINN")==0){
1704 ielset=i+1;
1705 break;
1706 }
1707 }
1708 if(ielset==0){
1709 printf("\n*ERROR in arpackcs: no node set for MAXE\n");
1710 printf(" (must have the name STRAINDOMAIN)\n\n");
1711 FORTRAN(stop,());
1712 }
1713
1714 for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
1715 if(ialset[i1]>0){
1716 l1=ialset[i1]-1;
1717
1718 worstpsmax=0.;
1719 for(l3=0;l3<360;l3++){
1720 ctl=cos(l3/constant);
1721 stl=sin(l3/constant);
1722 for(l2=0;l2<6;l2++){
1723 l=6*l1+l2;
1724 eenl[l2]=ctl*een[l]-stl*een[l+6**nk];
1725 }
1726
1727 /* determining the eigenvalues */
1728
1729 v1=eenl[0]+eenl[1]+eenl[2];
1730 v2=eenl[1]*eenl[2]+eenl[0]*eenl[2]+eenl[0]*eenl[1]-
1731 (eenl[5]*eenl[5]+eenl[4]*eenl[4]+eenl[3]*eenl[3]);
1732 v3=eenl[0]*(eenl[1]*eenl[2]-eenl[5]*eenl[5])
1733 -eenl[3]*(eenl[3]*eenl[2]-eenl[4]*eenl[5])
1734 +eenl[4]*(eenl[3]*eenl[5]-eenl[4]*eenl[1]);
1735 bb=v2-v1*v1/3.;
1736 cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1737 if(fabs(bb)<=1.e-10){
1738 if(fabs(cc)>1.e-10){
1739 al[0]=-pow(cc,(1./3.));
1740 }else{
1741 al[0]=0.;
1742 }
1743 al[1]=al[0];
1744 al[2]=al[0];
1745 }else{
1746 cm=2.*sqrt(-bb/3.);
1747 cn=3.*cc/(cm*bb);
1748 if(fabs(cn)>1.){
1749 if(cn>1.){
1750 cn=1.;
1751 }else{
1752 cn=-1.;
1753 }
1754 }
1755 tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1756 al[0]=cm*cos(tt);
1757 al[1]=cm*cos(tt+2.*pi/3.);
1758 al[2]=cm*cos(tt+4.*pi/3.);
1759 }
1760 for(l2=0;l2<3;l2++){
1761 al[l2]+=v1/3.;
1762 }
1763 dd=fabs(al[0]);
1764 if(fabs(al[1])>dd) dd=fabs(al[1]);
1765 if(fabs(al[2])>dd) dd=fabs(al[2]);
1766 if(dd>worstpsmax){
1767 worstpsmax=dd;
1768 eenmax[7*l1]=dd;
1769 for(l2=1;l2<7;l2++){
1770 eenmax[7*l1+l2]=eenl[l2-1];
1771 }
1772 }
1773 }
1774
1775 }else{
1776 l1=ialset[i1-2]-1;
1777 do{
1778 l1=l1-ialset[i1];
1779 if(l1>=ialset[i1-1]-1) break;
1780
1781 worstpsmax=0.;
1782 for(l3=0;l3<360;l3++){
1783 ctl=cos(l3/constant);
1784 stl=sin(l3/constant);
1785 for(l2=0;l2<6;l2++){
1786 l=6*l1+l2;
1787 eenl[l2]=ctl*een[l]-stl*een[l+6**nk];
1788 }
1789
1790 /* determining the eigenvalues */
1791
1792 v1=eenl[0]+eenl[1]+eenl[2];
1793 v2=eenl[1]*eenl[2]+eenl[0]*eenl[2]+eenl[0]*eenl[1]-
1794 (eenl[5]*eenl[5]+eenl[4]*eenl[4]+eenl[3]*eenl[3]);
1795 v3=eenl[0]*(eenl[1]*eenl[2]-eenl[5]*eenl[5])
1796 -eenl[3]*(eenl[3]*eenl[2]-eenl[4]*eenl[5])
1797 +eenl[4]*(eenl[3]*eenl[5]-eenl[4]*eenl[1]);
1798 bb=v2-v1*v1/3.;
1799 cc=-2.*v1*v1*v1/27.+v1*v2/3.-v3;
1800 if(fabs(bb)<=1.e-10){
1801 if(fabs(cc)>1.e-10){
1802 al[0]=-pow(cc,(1./3.));
1803 }else{
1804 al[0]=0.;
1805 }
1806 al[1]=al[0];
1807 al[2]=al[0];
1808 }else{
1809 cm=2.*sqrt(-bb/3.);
1810 cn=3.*cc/(cm*bb);
1811 if(fabs(cn)>1.){
1812 if(cn>1.){
1813 cn=1.;
1814 }else{
1815 cn=-1.;
1816 }
1817 }
1818 tt=(atan2(sqrt(1.-cn*cn),cn))/3.;
1819 al[0]=cm*cos(tt);
1820 al[1]=cm*cos(tt+2.*pi/3.);
1821 al[2]=cm*cos(tt+4.*pi/3.);
1822 }
1823 for(l2=0;l2<3;l2++){
1824 al[l2]+=v1/3.;
1825 }
1826 dd=fabs(al[0]);
1827 if(fabs(al[1])>dd) dd=fabs(al[1]);
1828 if(fabs(al[2])>dd) dd=fabs(al[2]);
1829 if(dd>worstpsmax){
1830 worstpsmax=dd;
1831 eenmax[7*l1]=dd;
1832 for(l2=1;l2<7;l2++){
1833 eenmax[7*l1+l2]=eenl[l2-1];
1834 }
1835 }
1836 }
1837
1838 }while(1);
1839 }
1840 }
1841 }
1842
1843 /* mapping the results to the other sectors */
1844
1845 for(l=0;l<*nk;l++){inumt[l]=inum[l];}
1846
1847 icntrl=2;imag=1;
1848
1849 FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
1850
1851 if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1852 for(l=0;l<mt**nk;l++){vt[l]=v[l];}
1853 for(l=0;l<mt**nk;l++){vt[l+mt**nk*ngraph]=v[l+mt**nk];}}
1854 if(strcmp1(&filab[87],"NT ")==0){
1855 if(*iperturb==0){
1856 for(l=0;l<*nk;l++){t1t[l]=t1[l];};
1857 }else{
1858 for(l=0;l<*nk;l++){t1t[l]=t1old[l];};
1859 }
1860 }
1861 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1862 for(l=0;l<6**nk;l++){stnt[l]=stn[l];}
1863 for(l=0;l<6**nk;l++){stnt[l+6**nk*ngraph]=stn[l+6**nk];}}
1864 if(strcmp1(&filab[261],"E ")==0){
1865 for(l=0;l<6**nk;l++){eent[l]=een[l];};
1866 for(l=0;l<6**nk;l++){eent[l+6**nk*ngraph]=een[l+6**nk];}}
1867 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
1868 for(l=0;l<mt**nk;l++){fnt[l]=fn[l];}
1869 for(l=0;l<mt**nk;l++){fnt[l+mt**nk*ngraph]=fn[l+mt**nk];}}
1870 if(strcmp1(&filab[522],"ENER")==0){
1871 for(l=0;l<*nk;l++){enernt[l]=enern[l];}}
1872 // for(l=0;l<*nk;l++){enernt[l+*nk*ngraph]=enern[l+*nk];}}
1873 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
1874 ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0))){
1875 for(l=0;l<6*mi[0]**ne;l++){stxt[l]=stx[l];}
1876 for(l=0;l<6*mi[0]**ne;l++){stxt[l+6*mi[0]**ne*ngraph]=stx[l+6*mi[0]**ne];}}
1877 if(strcmp1(&filab[2697],"ME ")==0){
1878 for(l=0;l<6**nk;l++){emnt[l]=emn[l];};
1879 for(l=0;l<6**nk;l++){emnt[l+6**nk*ngraph]=emn[l+6**nk];}}
1880 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
1881 &&(*mortar==1)){
1882 for(l=0;l<6**nk;l++){cdnt[l]=cdn[l];}
1883 for(l=0;l<6**nk;l++){cdnt[l+6**nk*ngraph]=cdn[l+6**nk];}}
1884
1885 for(jj=0;jj<*mcs;jj++){
1886 ilength=cs[17*jj+3];
1887 is=cs[17*jj+4];
1888 lprev=cs[17*jj+13];
1889 for(i=1;i<is;i++){
1890
1891 for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
1892
1893 theta=i*nm*2.*pi/cs[17*jj];
1894 ctl=cos(theta);
1895 stl=sin(theta);
1896
1897 if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1898 for(l1=0;l1<*nk;l1++){
1899 if(inocs[l1]==jj){
1900
1901 /* check whether node lies on axis */
1902
1903 ml1=-l1-1;
1904 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1905 if(id!=0){
1906 if(ics[lprev+id-1]==ml1){
1907 for(l2=0;l2<4;l2++){
1908 l=mt*l1+l2;
1909 vt[l+mt**nk*i]=v[l];
1910 }
1911 continue;
1912 }
1913 }
1914 for(l2=0;l2<4;l2++){
1915 l=mt*l1+l2;
1916 vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
1917 }
1918 }
1919 }
1920 }
1921
1922 /* imaginary part of the displacements in cylindrical
1923 coordinates */
1924
1925 if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)){
1926 for(l1=0;l1<*nk;l1++){
1927 if(inocs[l1]==jj){
1928
1929 /* check whether node lies on axis */
1930
1931 ml1=-l1-1;
1932 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1933 if(id!=0){
1934 if(ics[lprev+id-1]==ml1){
1935 for(l2=0;l2<4;l2++){
1936 l=mt*l1+l2;
1937 vt[l+mt**nk*(i+ngraph)]=v[l+mt**nk];
1938 }
1939 continue;
1940 }
1941 }
1942 for(l2=0;l2<4;l2++){
1943 l=mt*l1+l2;
1944 vt[l+mt**nk*(i+ngraph)]=stl*v[l]+ctl*v[l+mt**nk];
1945 }
1946 }
1947 }
1948 }
1949
1950 if(strcmp1(&filab[87],"NT ")==0){
1951 if(*iperturb==0){
1952 for(l=0;l<*nk;l++){
1953 if(inocs[l]==jj) t1t[l+*nk*i]=t1[l];
1954 }
1955 }else{
1956 for(l=0;l<*nk;l++){
1957 if(inocs[l]==jj) t1t[l+*nk*i]=t1old[l];
1958 }
1959 }
1960 }
1961
1962 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1963 for(l1=0;l1<*nk;l1++){
1964 if(inocs[l1]==jj){
1965
1966 /* check whether node lies on axis */
1967
1968 ml1=-l1-1;
1969 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1970 if(id!=0){
1971 if(ics[lprev+id-1]==ml1){
1972 for(l2=0;l2<6;l2++){
1973 l=6*l1+l2;
1974 stnt[l+6**nk*i]=stn[l];
1975 }
1976 continue;
1977 }
1978 }
1979 for(l2=0;l2<6;l2++){
1980 l=6*l1+l2;
1981 stnt[l+6**nk*i]=ctl*stn[l]-stl*stn[l+6**nk];
1982 }
1983 }
1984 }
1985 }
1986
1987 /* imaginary part of the stresses in cylindrical
1988 coordinates */
1989
1990 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)){
1991 for(l1=0;l1<*nk;l1++){
1992 if(inocs[l1]==jj){
1993
1994 /* check whether node lies on axis */
1995
1996 ml1=-l1-1;
1997 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
1998 if(id!=0){
1999 if(ics[lprev+id-1]==ml1){
2000 for(l2=0;l2<6;l2++){
2001 l=6*l1+l2;
2002 stnt[l+6**nk*(i+ngraph)]=stn[l+6**nk];
2003 }
2004 continue;
2005 }
2006 }
2007 for(l2=0;l2<6;l2++){
2008 l=6*l1+l2;
2009 stnt[l+6**nk*(i+ngraph)]=stl*stn[l]+ctl*stn[l+6**nk];
2010 }
2011 }
2012 }
2013 }
2014
2015 if(strcmp1(&filab[261],"E ")==0){
2016 for(l1=0;l1<*nk;l1++){
2017 if(inocs[l1]==jj){
2018
2019 /* check whether node lies on axis */
2020
2021 ml1=-l1-1;
2022 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2023 if(id!=0){
2024 if(ics[lprev+id-1]==ml1){
2025 for(l2=0;l2<6;l2++){
2026 l=6*l1+l2;
2027 eent[l+6**nk*i]=een[l];
2028 }
2029 continue;
2030 }
2031 }
2032 for(l2=0;l2<6;l2++){
2033 l=6*l1+l2;
2034 eent[l+6**nk*i]=ctl*een[l]-stl*een[l+6**nk];
2035 }
2036 }
2037 }
2038 }
2039
2040 /* imaginary part of the strains in cylindrical
2041 coordinates */
2042
2043 if(strcmp1(&filab[261],"E ")==0){
2044 for(l1=0;l1<*nk;l1++){
2045 if(inocs[l1]==jj){
2046
2047 /* check whether node lies on axis */
2048
2049 ml1=-l1-1;
2050 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2051 if(id!=0){
2052 if(ics[lprev+id-1]==ml1){
2053 for(l2=0;l2<6;l2++){
2054 l=6*l1+l2;
2055 eent[l+6**nk*(i+ngraph)]=een[l+6**nk];
2056 }
2057 continue;
2058 }
2059 }
2060 for(l2=0;l2<6;l2++){
2061 l=6*l1+l2;
2062 eent[l+6**nk*(i+ngraph)]=stl*een[l]+ctl*een[l+6**nk];
2063 }
2064 }
2065 }
2066 }
2067
2068 /* real part of the mechanical strains */
2069
2070 if(strcmp1(&filab[2697],"ME ")==0){
2071 for(l1=0;l1<*nk;l1++){
2072 if(inocs[l1]==jj){
2073
2074 /* check whether node lies on axis */
2075
2076 ml1=-l1-1;
2077 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2078 if(id!=0){
2079 if(ics[lprev+id-1]==ml1){
2080 for(l2=0;l2<6;l2++){
2081 l=6*l1+l2;
2082 emnt[l+6**nk*i]=emn[l];
2083 }
2084 continue;
2085 }
2086 }
2087 for(l2=0;l2<6;l2++){
2088 l=6*l1+l2;
2089 emnt[l+6**nk*i]=ctl*emn[l]-stl*emn[l+6**nk];
2090 }
2091 }
2092 }
2093 }
2094
2095 /* imaginary part of the mechanical strains in cylindrical
2096 coordinates */
2097
2098 if(strcmp1(&filab[2697],"ME ")==0){
2099 for(l1=0;l1<*nk;l1++){
2100 if(inocs[l1]==jj){
2101
2102 /* check whether node lies on axis */
2103
2104 ml1=-l1-1;
2105 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2106 if(id!=0){
2107 if(ics[lprev+id-1]==ml1){
2108 for(l2=0;l2<6;l2++){
2109 l=6*l1+l2;
2110 emnt[l+6**nk*(i+ngraph)]=emn[l+6**nk];
2111 }
2112 continue;
2113 }
2114 }
2115 for(l2=0;l2<6;l2++){
2116 l=6*l1+l2;
2117 emnt[l+6**nk*(i+ngraph)]=stl*emn[l]+ctl*emn[l+6**nk];
2118 }
2119 }
2120 }
2121 }
2122
2123 /* real part of the contact stresses */
2124
2125 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2126 &&(*mortar==1)){
2127 for(l1=0;l1<*nk;l1++){
2128 if(inocs[l1]==jj){
2129
2130 /* check whether node lies on axis */
2131
2132 ml1=-l1-1;
2133 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2134 if(id!=0){
2135 if(ics[lprev+id-1]==ml1){
2136 for(l2=0;l2<6;l2++){
2137 l=6*l1+l2;
2138 cdnt[l+6**nk*i]=cdn[l];
2139 }
2140 continue;
2141 }
2142 }
2143 for(l2=0;l2<6;l2++){
2144 l=6*l1+l2;
2145 cdnt[l+6**nk*i]=ctl*cdn[l]-stl*cdn[l+6**nk];
2146 }
2147 }
2148 }
2149 }
2150
2151 /* imaginary part of the contact stresses */
2152
2153 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2154 &&(*mortar==1)){
2155 for(l1=0;l1<*nk;l1++){
2156 if(inocs[l1]==jj){
2157
2158 /* check whether node lies on axis */
2159
2160 ml1=-l1-1;
2161 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2162 if(id!=0){
2163 if(ics[lprev+id-1]==ml1){
2164 for(l2=0;l2<6;l2++){
2165 l=6*l1+l2;
2166 cdnt[l+6**nk*(i+ngraph)]=cdn[l+6**nk];
2167 }
2168 continue;
2169 }
2170 }
2171 for(l2=0;l2<6;l2++){
2172 l=6*l1+l2;
2173 cdnt[l+6**nk*(i+ngraph)]=stl*cdn[l]+ctl*cdn[l+6**nk];
2174 }
2175 }
2176 }
2177 }
2178
2179 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
2180 for(l1=0;l1<*nk;l1++){
2181 if(inocs[l1]==jj){
2182
2183 /* check whether node lies on axis */
2184
2185 ml1=-l1-1;
2186 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2187 if(id!=0){
2188 if(ics[lprev+id-1]==ml1){
2189 for(l2=0;l2<4;l2++){
2190 l=mt*l1+l2;
2191 fnt[l+mt**nk*i]=fn[l];
2192 }
2193 continue;
2194 }
2195 }
2196 for(l2=0;l2<4;l2++){
2197 l=mt*l1+l2;
2198 fnt[l+mt**nk*i]=ctl*fn[l]-stl*fn[l+mt**nk];
2199 }
2200 }
2201 }
2202 }
2203
2204 /* imaginary part of the forces in cylindrical
2205 coordinates */
2206
2207 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)){
2208 for(l1=0;l1<*nk;l1++){
2209 if(inocs[l1]==jj){
2210
2211 /* check whether node lies on axis */
2212
2213 ml1=-l1-1;
2214 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
2215 if(id!=0){
2216 if(ics[lprev+id-1]==ml1){
2217 for(l2=0;l2<4;l2++){
2218 l=mt*l1+l2;
2219 fnt[l+mt**nk*(i+ngraph)]=fn[l+mt**nk];
2220 }
2221 continue;
2222 }
2223 }
2224 for(l2=0;l2<4;l2++){
2225 l=mt*l1+l2;
2226 fnt[l+mt**nk*(i+ngraph)]=stl*fn[l]+ctl*fn[l+mt**nk];
2227 }
2228 }
2229 }
2230 }
2231
2232 }
2233 }
2234
2235 icntrl=-2;imag=0;
2236
2237 FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
2238 &imag,mi,emnt));
2239
2240 FORTRAN(rectcylvi,(cot,&vt[mt**nk*ngraph],&fnt[mt**nk*ngraph],
2241 &stnt[6**nk*ngraph],qfnt,&eent[6**nk*ngraph],cs,&nkt,&icntrl,
2242 t,filab,&imag,mi,&emnt[6**nk*ngraph]));
2243
2244 /* internal energy calculation */
2245
2246 for(jj=0;jj<*mcs;jj++){
2247 ilength=cs[17*jj+3];
2248 is=cs[17*jj+4];
2249 lprev=cs[17*jj+13];
2250 for(i=1;i<is;i++){
2251
2252 for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
2253
2254 theta=i*nm*2.*pi/cs[17*jj];
2255 ctl=cos(theta);
2256 stl=sin(theta);
2257
2258 if(strcmp1(&filab[522],"ENER")==0){
2259 ioffr=6**nk*i;
2260 for(l1=0;l1<*nk;l1++){
2261 ioffrl=6*l1+ioffr;
2262 enernt[l1+*nk*i]=(stnt[ioffrl]*emnt[ioffrl]+
2263 stnt[1+ioffrl]*emnt[1+ioffrl]+
2264 stnt[2+ioffrl]*emnt[2+ioffrl])/2.+
2265 stnt[3+ioffrl]*emnt[3+ioffrl]+
2266 stnt[4+ioffrl]*emnt[4+ioffrl]+
2267 stnt[5+ioffrl]*emnt[5+ioffrl];
2268 }
2269 }
2270
2271 }
2272 }
2273
2274 /* determining magnitude and phase angle for the displacements */
2275
2276 if(strcmp1(&filab[870],"PU")==0){
2277 for(l1=0;l1<nkt;l1++){
2278 for(l2=0;l2<4;l2++){
2279 l=mt*l1+l2;
2280 vreal=vt[l];
2281 vimag=vt[l+mt**nk*ngraph];
2282 vr[l]=sqrt(vreal*vreal+vimag*vimag);
2283 if(vr[l]<1.e-10){
2284 vi[l]=0.;
2285 }else if(fabs(vreal)<1.e-10){
2286 if(vimag>0){vi[l]=90.;}
2287 else{vi[l]=-90.;}
2288 }
2289 else{
2290 vi[l]=atan(vimag/vreal)*constant;
2291 if(vreal<0) vi[l]+=180.;
2292 }
2293 }
2294 }
2295 }
2296
2297 /* determining magnitude and phase for the stress */
2298
2299 if(strcmp1(&filab[1479],"PHS")==0){
2300 for(l1=0;l1<nkt;l1++){
2301 for(l2=0;l2<6;l2++){
2302 l=6*l1+l2;
2303 stnreal=stnt[l];
2304 stnimag=stnt[l+6**nk*ngraph];
2305 stnr[l]=sqrt(stnreal*stnreal+stnimag*stnimag);
2306 if(stnr[l]<1.e-10){
2307 stni[l]=0.;
2308 }else if(fabs(stnreal)<1.e-10){
2309 if(stnimag>0){stni[l]=90.;}
2310 else{stni[l]=-90.;}
2311 }
2312 else{
2313 stni[l]=atan(stnimag/stnreal)*constant;
2314 if(stnreal<0) stni[l]+=180.;
2315 }
2316 }
2317 }
2318 }
2319
2320 /* determining magnitude and phase angle for the forces */
2321
2322 if(strcmp1(&filab[2610],"PRF")==0){
2323 for(l1=0;l1<nkt;l1++){
2324 for(l2=0;l2<4;l2++){
2325 l=mt*l1+l2;
2326 fnreal=fnt[l];
2327 fnimag=fnt[l+mt**nk*ngraph];
2328 fnr[l]=sqrt(fnreal*fnreal+fnimag*fnimag);
2329 if(fnr[l]<1.e-10){
2330 fni[l]=0.;
2331 }else if(fabs(fnreal)<1.e-10){
2332 if(fnimag>0){fni[l]=90.;}
2333 else{fni[l]=-90.;}
2334 }
2335 else{
2336 fni[l]=atan(fnimag/fnreal)*constant;
2337 if(fnreal<0) fni[l]+=180.;
2338 }
2339 }
2340 }
2341 }
2342
2343 /* determining magnitude and phase for the contact stress */
2344
2345 if((strcmp1(&filab[3915],"PCON")==0)&&(*mortar==1)){
2346 for(l1=0;l1<nkt;l1++){
2347 for(l2=0;l2<6;l2++){
2348 l=6*l1+l2;
2349 cdnreal=cdnt[l];
2350 cdnimag=cdnt[l+6**nk*ngraph];
2351 cdnr[l]=sqrt(cdnreal*cdnreal+cdnimag*cdnimag);
2352 if(cdnr[l]<1.e-10){
2353 cdni[l]=0.;
2354 }else if(fabs(cdnreal)<1.e-10){
2355 if(cdnimag>0){cdni[l]=90.;}
2356 else{cdni[l]=-90.;}
2357 }
2358 else{
2359 cdni[l]=atan(cdnimag/cdnreal)*constant;
2360 if(cdnreal<0) cdni[l]+=180.;
2361 }
2362 }
2363 }
2364 }
2365
2366 ++*kode;
2367 if(d[j]>=0.){
2368 freq=sqrt(d[j])/6.283185308;
2369 }else{
2370 freq=0.;
2371 }
2372
2373 if(strcmp1(&filab[1044],"ZZS")==0){
2374 NNEW(neigh,ITG,40*net);
2375 NNEW(ipneigh,ITG,nkt);
2376 }
2377
2378 /* deactivating S and ME request if only needed for ENER */
2379
2380 if((ngraph>1)&&(strcmp1(&filab[522],"ENER")==0)){
2381 strcpy1(&filab[174],&filabcp[0],4);
2382 strcpy1(&filab[2697],&filabcp[4],4);
2383 }
2384
2385 if(igreen==1) *nmethod=13;
2386
2387 /* change on 20210510: if the nodal diameter exeeds half the
2388 number of segments change the sign of the axis */
2389
2390 // if(nm>cs[0]/2){
2391 if(nm-floor(nm/cs[0])*cs[0]>cs[0]/2){
2392 for(k=5;k<11;k++){cs[k]=-cs[k];}
2393 }
2394 frd(cot,&nkt,kont,ipkont,lakont,&net,vt,stnt,inumt,nmethod,
2395 kode,filab,eent,t1t,fnt,&freq,epn,ielmatt,matname,enernt,xstaten,
2396 nstate_,istep,&iinc,ithermal,qfn,&j,&nm,trab,inotrt,
2397 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
2398 mi,stxt,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&net,
2399 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emnt,
2400 thicke,jobnamec,output,qfx,cdnt,mortar,cdnr,cdni,nmat,
2401 ielprop,prop,sti);
2402 if(nm>cs[0]/2){
2403 for(k=5;k<11;k++){cs[k]=-cs[k];}
2404 }
2405
2406 if(igreen==1) *nmethod=2;
2407
2408 /* reactivating S and ME request if only needed for ENER */
2409
2410 if((ngraph>1)&&(strcmp1(&filab[522],"ENER")==0)){
2411 strcpy1(&filabcp[0],&filab[174],4);
2412 strcpy1(&filabcp[4],&filab[2697],4);
2413 strcpy1(&filab[174],"S ",4);
2414 strcpy1(&filab[2697],"ME ",4);
2415 }
2416
2417 if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
2418
2419 if(*iperturb!=0){
2420 for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
2421 xstate[k]=xstateini[k];
2422 }
2423 }
2424
2425 /* end loop over the eigenvalues */
2426
2427 }
2428
2429 /*--------------------------------------------------------------------*/
2430 /*
2431 -----------
2432 free memory
2433 -----------
2434 */
2435 if(*isolver==0){
2436 #ifdef SPOOLES
2437 spooles_cleanup();
2438 #endif
2439 }
2440 else if(*isolver==4){
2441 #ifdef SGI
2442 sgi_cleanup(token);
2443 #endif
2444 }
2445 else if(*isolver==5){
2446 #ifdef TAUCS
2447 tau_cleanup();
2448 #endif
2449 }
2450 else if(*isolver==7){
2451 #ifdef PARDISO
2452 pardiso_cleanup(&neq[1],&symmetryflag,&inputformat);
2453 #endif
2454 }
2455 else if(*isolver==8){
2456 #ifdef PASTIX
2457 #endif
2458 }
2459
2460 /* output of the turning direction */
2461
2462 FORTRAN(writeturdircs,(xn,turdir,&nev,&nm));
2463 SFREE(turdir);
2464
2465 if(igreen==0){
2466 if((fmax>-0.5)&&(fmax*fmax>d[nev-1])){
2467 printf("\n*WARNING: not all frequencies in the requested interval might be found;\nincrease the number of requested frequencies\n");
2468 }
2469 }
2470
2471 SFREE(adb);SFREE(aub);SFREE(temp_array);SFREE(coefmpcnew);SFREE(ad);SFREE(au);
2472
2473 /* deactivating S and ME request if only needed for ENER */
2474
2475 if((ngraph>1)&&(strcmp1(&filab[522],"ENER")==0)){
2476 strcpy1(&filab[174],&filabcp[0],4);
2477 strcpy1(&filab[2697],&filabcp[4],4);
2478 if(strcmp1(&filab[174],"S ")!=0){SFREE(stn);SFREE(stnt);}
2479 if(strcmp1(&filab[2697],"ME ")!=0){SFREE(emn);SFREE(emnt);}
2480 }
2481
2482 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1653],"MAXS")==0)||
2483 (strcmp1(&filab[1479],"PHS ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
2484 (strcmp1(&filab[1044],"ERR ")==0))
2485 SFREE(stn);
2486
2487 SFREE(v);SFREE(fn);SFREE(inum);SFREE(stx);SFREE(eme);SFREE(z);SFREE(d);
2488
2489 if((strcmp1(&filab[261],"E ")==0)||(strcmp1(&filab[2523],"MAXE")==0)) SFREE(een);
2490 if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
2491 if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
2492
2493 if((strcmp1(&filab[0],"U ")==0)||(strcmp1(&filab[870],"PU ")==0)) SFREE(vt);
2494 if(strcmp1(&filab[87],"NT ")==0) SFREE(t1t);
2495 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1479],"PHS ")==0)||
2496 (strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)) SFREE(stnt);
2497 if(strcmp1(&filab[261],"E ")==0) SFREE(eent);
2498 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[2610],"PRF ")==0)) SFREE(fnt);
2499 if(strcmp1(&filab[522],"ENER")==0) SFREE(enernt);
2500 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
2501 ((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==0))) SFREE(stxt);
2502 if(strcmp1(&filab[2697],"ME ")==0) SFREE(emnt);
2503 if(((strcmp1(&filab[2175],"CONT")==0)||(strcmp1(&filab[3915],"PCON")==0))
2504 &&(*mortar==1)) SFREE(cdnt);
2505
2506 SFREE(cot);SFREE(kont);SFREE(ipkont);SFREE(lakont);SFREE(inumt);SFREE(ielmatt);
2507 if(*ntrans>0){SFREE(inotrt);}
2508
2509 if(mei[3]==1){
2510 (*nevtot)+=nev;
2511 fclose(f1);
2512 }
2513
2514 /* end loop over the nodal diameters */
2515
2516 }
2517
2518 if(*iperturb!=0){
2519 if(ncont!=0){
2520 *ne=ne0;*nkon=nkon0;
2521 // if(*nener==1){SFREE(ener);}
2522 RENEW(ipkon,ITG,*ne);
2523 RENEW(lakon,char,8**ne);
2524 RENEW(kon,ITG,*nkon);
2525 if(*norien>0){
2526 RENEW(ielorien,ITG,mi[2]**ne);
2527 }
2528 RENEW(ielmat,ITG,mi[2]**ne);
2529 SFREE(cg);
2530 SFREE(straight);
2531 SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
2532 SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
2533 SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
2534 SFREE(springarea);SFREE(xmastnor);
2535
2536 if(*mortar==0){
2537 SFREE(areaslav);
2538 }else if(*mortar==1){
2539 SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
2540 SFREE(islavact);
2541 }
2542 }
2543 }
2544
2545 if(*nener==1){SFREE(ener);}
2546 SFREE(inocs);SFREE(ielcs);SFREE(xstiff);
2547
2548 // if(*nbody>0) SFREE(ipobody);
2549
2550 if(*nstate_!=0) SFREE(xstateini);
2551
2552 if(strcmp1(&filab[870],"PU")==0){SFREE(vr);SFREE(vi);}
2553 if(strcmp1(&filab[1479],"PHS")==0){SFREE(stnr);SFREE(stni);}
2554 if(strcmp1(&filab[3915],"PCON")==0){SFREE(cdnr);SFREE(cdni);}
2555 if(strcmp1(&filab[1566],"MAXU")==0){SFREE(vmax);}
2556 if(strcmp1(&filab[1653],"MAXS")==0){SFREE(stnmax);}
2557 if(strcmp1(&filab[2523],"MAXE")==0){SFREE(eenmax);}
2558 if(strcmp1(&filab[2610],"PRF")==0){SFREE(fnr);SFREE(fni);}
2559
2560 if(*iperturb!=0){
2561 mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
2562 mpcinfo[3]=maxlenmpc;
2563 }
2564
2565 *irowp=irow;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
2566 *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;
2567
2568 *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
2569
2570 return;
2571 }
2572
2573
2574 #endif
2575