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
arpack(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 * shcon,ITG * nshcon,double * cocon,ITG * ncocon,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 * ncmat_,ITG * nstate_,double ** enerp,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,double * thicke,ITG * nslavs,double * tietol,ITG * nkon,ITG * mpcinfo,ITG * ntie,ITG * istep,ITG * mcs,ITG * ics,char * tieset,double * cs,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 arpack(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 *shcon, ITG *nshcon, double *cocon, ITG *ncocon,
57 double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
58 ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
59 double *t0, double *t1, double *t1old,
60 ITG *ithermal,double *prestr, ITG *iprestr,
61 double *vold,ITG *iperturb, double *sti, ITG *nzs,
62 ITG *kode, ITG *mei, double *fei,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 *ncmat_, ITG *nstate_, double **enerp, char *jobnamec,
67 char *output, char *set, ITG *nset, ITG *istartset,
68 ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
69 char *prset, ITG *nener, ITG *isolver, double *trab,
70 ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
71 ITG *ipobody, ITG *ibody,double *xbody, ITG *nbody,
72 double *thicke, ITG *nslavs, double *tietol, ITG *nkon,
73 ITG *mpcinfo,ITG *ntie,ITG *istep,ITG *mcs,ITG *ics,
74 char *tieset,double *cs,ITG *nintpoint,ITG *mortar,
75 ITG *ifacecount,ITG **islavsurfp,double **pslavsurfp,
76 double **clearinip,ITG *nmat,char *typeboun,
77 ITG *ielprop,double *prop,char *orname,ITG *inewton,
78 double *t0g,double *t1g){
79
80 /* calls the Arnoldi Package (ARPACK) */
81
82 char bmat[2]="G", which[3]="LM", howmny[2]="A", fneig[132]="",
83 description[13]=" ",*lakon=NULL,jobnamef[396]="",*labmpc2=NULL;
84
85 ITG *inum=NULL,k,ido,ldz,iparam[11],ipntr[14],lworkl,ngraph=1,im,
86 info,rvec=1,*select=NULL,lfin,j,lint,iout,ielas=0,icmd=0,mt=mi[1]+1,
87 iinc=1,nev,ncv,mxiter,jrow,ifreebody,kmax,
88 mass[2]={1,1}, stiffness=1, buckling=0, rhsi=0, intscheme=0,noddiam=-1,
89 coriolis=0,symmetryflag=0,inputformat=0,*ipneigh=NULL,*neigh=NULL,ne0,
90 *integerglob=NULL,nasym=0,zero=0,ncont=0,*itietri=NULL,kref,
91 *koncont=NULL,ismallsliding=0,*itiefac=NULL,*islavsurf=NULL,
92 *islavnode=NULL,*imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,
93 *imastop=NULL,*iponoels=NULL,*inoels=NULL,*ipe=NULL,*ime=NULL,
94 mpcfree,memmpc_,icascade,maxlenmpc,iit=-1,*irow=NULL,nherm=1,
95 icfd=0,*inomat=NULL,*ipkon=NULL,*kon=NULL,*ielmat=NULL,*ielorien=NULL,
96 *islavact=NULL,maxprevcontel,iex,nslavs_prev_step,icutb=0,
97 iflagact=0,*islavsurfold=NULL,ialeatoric=0,kscale=1,network=0,
98 *iponoel=NULL,*inoel=NULL,nrhs=1,igreen=0,node,idir,mscalmethod=0,
99 *jqw=NULL,*iroww=NULL,nzsw,
100 *islavelinv=NULL,*irowtloc=NULL,*jqtloc=NULL,nboun2,
101 *ndirboun2=NULL,*nodeboun2=NULL,nmpc2,*ipompc2=NULL,*nodempc2=NULL,
102 *ikboun2=NULL,*ilboun2=NULL,*ikmpc2=NULL,*ilmpc2=NULL,mortartrafoflag=0;
103
104 double *stn=NULL,*v=NULL,*resid=NULL,*z=NULL,*workd=NULL,
105 *workl=NULL,*d=NULL,sigma=1,*temp_array=NULL,*pslavsurfold=NULL,
106 *een=NULL,sum,cam[5],*f=NULL,*fn=NULL,qa[4],*fext=NULL,
107 *epn=NULL,*fnr=NULL,*fni=NULL,*emn=NULL,*emeini=NULL,
108 *xstateini=NULL,*xstiff=NULL,*stiini=NULL,*vini=NULL,freq,*stx=NULL,
109 *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,*fnext=NULL,
110 *physcon=NULL,*qfx=NULL,*qfn=NULL,tol,fmin,fmax,pi,*cgr=NULL,
111 *xloadold=NULL,reltime=1.,*vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,
112 *vmax=NULL,*stnmax=NULL,*springarea=NULL,*eenmax=NULL,
113 *doubleglob=NULL,*cg=NULL,*straight=NULL,*clearini=NULL,
114 *xmastnor=NULL,*areaslav=NULL,*xnoels=NULL,theta=0.,
115 *di=NULL,sigmai=0,*workev=NULL,*ener=NULL,*xstate=NULL,*dc=NULL,
116 *au=NULL,*ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,
117 *pmastsurf=NULL,*cdn=NULL,*energyini=NULL,*energy=NULL,
118 *cdnr=NULL,*cdni=NULL,*eme=NULL,alea=0.1,*smscale=NULL,zmax,*auw=NULL,
119 *autloc=NULL,*xboun2=NULL,*coefmpc2=NULL;
120
121 FILE *f1;
122
123 /* dummy arguments for the results call */
124
125 double *veold=NULL,*accold=NULL,bet,gam,dtime,time=1.;
126
127 #ifdef SGI
128 ITG token;
129 #endif
130
131 irow=*irowp;ener=*enerp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
132 kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;
133
134 islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
135
136 if(*nmethod==13){
137 *nmethod=2;
138 igreen=1;
139 }
140
141 for(k=0;k<3;k++){
142 strcpy1(&jobnamef[k*132],&jobnamec[k*132],132);
143 }
144
145 if(*mortar!=1){
146 maxprevcontel=*nslavs;
147 }else if(*mortar==1){
148 maxprevcontel=*nintpoint;
149 if(*nstate_!=0){
150 if(maxprevcontel!=0){
151 NNEW(islavsurfold,ITG,2**ifacecount+2);
152 NNEW(pslavsurfold,double,3**nintpoint);
153 memcpy(&islavsurfold[0],&islavsurf[0],
154 sizeof(ITG)*(2**ifacecount+2));
155 memcpy(&pslavsurfold[0],&pslavsurf[0],
156 sizeof(double)*(3**nintpoint));
157 }
158 }
159 }
160 nslavs_prev_step=*nslavs;
161
162 if((strcmp1(&filab[870],"PU ")==0)||
163 (strcmp1(&filab[1479],"PHS ")==0)||
164 (strcmp1(&filab[1566],"MAXU")==0)||
165 (strcmp1(&filab[1653],"MAXS")==0)){
166 printf("*WARNING in arpack: PU, PHS, MAXU or MAX was selected in a frequency calculation without cyclic symmetry;\n this is not correct; output request is removed;\n");
167 strcpy1(&filab[870]," ",4);
168 strcpy1(&filab[1479]," ",4);
169 strcpy1(&filab[1566]," ",4);
170 strcpy1(&filab[1653]," ",4);
171 }
172
173 /* copying the frequency parameters */
174
175 pi=4.*atan(1.);
176
177 if(igreen==0){
178 nev=mei[0];
179 ncv=mei[1];
180 mxiter=mei[2];
181 tol=fei[0];
182 fmin=2*pi*fei[1];
183 fmax=2*pi*fei[2];
184 }else{
185
186 /* Green functions: number of "eigenmodes" is the number of
187 unit forces */
188
189 nev=*nforc;
190 ncv=*nforc;
191 }
192
193 /* assigning the body forces to the elements */
194
195 if(*nbody>0){
196 /* ifreebody=*ne+1;
197 NNEW(ipobody,ITG,2*ifreebody**nbody);
198 for(k=1;k<=*nbody;k++){
199 FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
200 iendset,ialset,&inewton,nset,&ifreebody,&k));
201 RENEW(ipobody,ITG,2*(*ne+ifreebody));
202 }
203 RENEW(ipobody,ITG,2*(ifreebody-1));*/
204 if(*inewton==1){
205 printf("*ERROR in arpack: generalized gravity loading is not allowed in frequency calculations");
206 FORTRAN(stop,());
207 }
208 }
209
210 ne0=*ne;
211
212 /* contact conditions */
213
214 if(*iperturb!=0){
215
216 memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
217 maxlenmpc=mpcinfo[3];
218
219 inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
220 lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
221 &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
222 mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
223 iperturb,ikboun,nboun,co,istep,&xnoels);
224
225 if(ncont!=0){
226
227 NNEW(cg,double,3*ncont);
228 NNEW(straight,double,16*ncont);
229
230 /* 11 instead of 10: last position is reserved for the
231 local contact spring element number; needed as
232 pointer into springarea */
233
234 if(*mortar==0){
235 RENEW(kon,ITG,*nkon+11**nslavs);
236 NNEW(springarea,double,2**nslavs);
237 if(*nener==1){
238 RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
239 }
240 RENEW(ipkon,ITG,*ne+*nslavs);
241 RENEW(lakon,char,8*(*ne+*nslavs));
242
243 if(*norien>0){
244 RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
245 for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
246 }
247
248 RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
249 for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
250
251 if((maxprevcontel==0)&&(*nslavs!=0)){
252 RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
253 for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
254 xstate[k]=0.;
255 }
256 }
257 maxprevcontel=*nslavs;
258
259 NNEW(areaslav,double,*ifacecount);
260 NNEW(xmastnor,double,3*nmastnode[*ntie]);
261 }else if(*mortar==1){
262 NNEW(islavact,ITG,nslavnode[*ntie]);
263 DMEMSET(islavact,0,nslavnode[*ntie],1);
264 if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
265 NNEW(xmastnor,double,3*nmastnode[*ntie]);
266
267 *nintpoint=0;
268
269 precontact(&ncont,ntie,tieset,nset,set,istartset,
270 iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
271 cg,straight,co,vold,istep,&iinc,&iit,itiefac,
272 islavsurf,islavnode,imastnode,nslavnode,nmastnode,
273 imastop,mi,ipe,ime,tietol,&iflagact,
274 nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
275 nslavs);
276
277 /* changing the dimension of element-related fields */
278
279 RENEW(kon,ITG,*nkon+22**nintpoint);
280 RENEW(springarea,double,2**nintpoint);
281 RENEW(pmastsurf,double,6**nintpoint);
282
283 if(*nener==1){
284 RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
285 }
286 RENEW(ipkon,ITG,*ne+*nintpoint);
287 RENEW(lakon,char,8*(*ne+*nintpoint));
288
289 if(*norien>0){
290 RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
291 for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
292 }
293 RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
294 for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
295
296 /* interpolating the state variables */
297
298 if(*nstate_!=0){
299 if(maxprevcontel!=0){
300 RENEW(xstateini,double,
301 *nstate_*mi[0]*(ne0+maxprevcontel));
302 for(k=*nstate_*mi[0]*ne0;
303 k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
304 xstateini[k]=xstate[k];
305 }
306 }
307
308 RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
309 for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
310 xstate[k]=0.;
311 }
312
313 if((*nintpoint>0)&&(maxprevcontel>0)){
314 iex=2;
315
316 /* interpolation of xstate */
317
318 FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
319 &ne0,mi,xstate,pslavsurf,nstate_,
320 xstateini,islavsurf,islavsurfold,
321 pslavsurfold,tieset,ntie,itiefac));
322
323 }
324
325 if(maxprevcontel!=0){
326 SFREE(islavsurfold);SFREE(pslavsurfold);
327 }
328
329 maxprevcontel=*nintpoint;
330
331 RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
332 for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
333 xstateini[k]=xstate[k];
334 }
335 }
336 }
337
338 /* generating contact spring elements */
339
340 contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
341 ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
342 co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
343 &ne0,vini,nmethod,
344 iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
345 itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
346 imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
347 xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta,
348 xstateini,xstate,nstate_,&icutb,&ialeatoric,jobnamef,&alea,
349 auw,jqw,iroww,&nzsw);
350
351 printf("number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
352
353 /* determining the structure of the stiffness/mass matrix */
354
355 remastructar(ipompc,&coefmpc,&nodempc,nmpc,
356 &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
357 labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
358 kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
359 neq,nzs,nmethod,ithermal,iperturb,mass,mi,ics,cs,
360 mcs,mortar,typeboun,&iit,&network,iexpl);
361 }
362 }
363
364 /* field for initial values of state variables (needed if
365 previous static step was viscoplastic and for contact */
366
367 if((*nstate_!=0)&&((*mortar==0)||(ncont==0))){
368 NNEW(xstateini,double,*nstate_*mi[0]*(ne0+*nslavs));
369 for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
370 xstateini[k]=xstate[k];
371 }
372 }
373
374 /* determining the internal forces and the stiffness coefficients */
375
376 NNEW(f,double,neq[1]);
377
378 /* allocating a field for the stiffness matrix */
379
380 NNEW(xstiff,double,(long long)27*mi[0]**ne);
381
382 iout=-1;
383 NNEW(v,double,mt**nk);
384 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
385 NNEW(fn,double,mt**nk);
386 NNEW(stx,double,6*mi[0]**ne);
387 NNEW(eme,double,6*mi[0]**ne);
388 if(*ithermal>1){
389 NNEW(qfx,double,3*mi[0]**ne);
390 }
391 NNEW(inum,ITG,*nk);
392 if(*iperturb==0){
393 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
394 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
395 ielorien,norien,orab,ntmat_,t0,t0,ithermal,
396 prestr,iprestr,filab,eme,emn,een,iperturb,
397 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
398 ndirboun,xboun,nboun,ipompc,
399 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
400 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
401 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
402 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
403 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
404 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
405 fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
406 &reltime,&ne0,thicke,shcon,nshcon,
407 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
408 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
409 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
410 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
411 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
412 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
413 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
414 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
415 &intscheme);
416 }else{
417 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
418 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
419 ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
420 prestr,iprestr,filab,eme,emn,een,iperturb,
421 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
422 ndirboun,xboun,nboun,ipompc,
423 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
424 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
425 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
426 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
427 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
428 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
429 fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
430 &reltime,&ne0,thicke,shcon,nshcon,
431 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
432 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
433 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
434 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
435 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
436 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
437 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
438 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
439 &intscheme);
440 }
441 SFREE(f);SFREE(v);SFREE(fn);SFREE(stx);SFREE(eme);
442 if(*ithermal>1)SFREE(qfx);SFREE(inum);
443 iout=1;
444
445 /* for the frequency analysis linear strain and elastic properties
446 are used */
447
448 iperturb[1]=0;ielas=1;
449
450 /* filling in the matrix */
451
452 NNEW(ad,double,neq[1]);
453 NNEW(au,double,nzs[2]);
454
455 NNEW(adb,double,neq[1]);
456 NNEW(aub,double,nzs[1]);
457
458 NNEW(fext,double,neq[1]);
459
460 if(*iperturb==0){
461 mafillsmmain(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
462 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
463 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
464 ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
465 ikmpc,ilmpc,ikboun,ilboun,
466 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
467 ielorien,norien,orab,ntmat_,
468 t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
469 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
470 xstiff,npmat_,&dtime,matname,mi,
471 ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
472 physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
473 &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
474 xstateini,xstate,thicke,integerglob,doubleglob,
475 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
476 pmastsurf,mortar,clearini,ielprop,prop,&ne0,fnext,&kscale,
477 iponoel,inoel,&network,ntrans,inotr,trab,smscale,&mscalmethod,
478 set,nset,islavelinv,autloc,irowtloc,jqtloc,&mortartrafoflag);
479 }
480 else{
481 mafillsmmain(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
482 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
483 nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
484 ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
485 ikmpc,ilmpc,ikboun,ilboun,
486 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
487 ielorien,norien,orab,ntmat_,
488 t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
489 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
490 xstiff,npmat_,&dtime,matname,mi,
491 ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
492 physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
493 &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
494 xstateini,xstate,thicke,integerglob,doubleglob,
495 tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
496 pmastsurf,mortar,clearini,ielprop,prop,&ne0,fnext,&kscale,
497 iponoel,inoel,&network,ntrans,inotr,trab,smscale,&mscalmethod,
498 set,nset,islavelinv,autloc,irowtloc,jqtloc,&mortartrafoflag);
499
500 if(nasym==1){
501 RENEW(au,double,nzs[2]+nzs[1]);
502 RENEW(aub,double,nzs[2]+nzs[1]);
503 symmetryflag=2;
504 inputformat=1;
505
506 mafillsmasmain(co,nk,kon,ipkon,lakon,ne,nodeboun,
507 ndirboun,xboun,nboun,
508 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
509 nforc,nelemload,sideload,xload,nload,xbody,ipobody,
510 nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
511 nmethod,ikmpc,ilmpc,ikboun,ilboun,
512 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
513 ielmat,ielorien,norien,orab,ntmat_,
514 t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
515 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
516 xstiff,npmat_,&dtime,matname,mi,
517 ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
518 physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
519 &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
520 xstateini,xstate,thicke,
521 integerglob,doubleglob,tieset,istartset,iendset,
522 ialset,ntie,&nasym,pslavsurf,pmastsurf,mortar,clearini,
523 ielprop,prop,&ne0,&kscale,iponoel,inoel,&network,set,nset);
524
525 /*FORTRAN(mafillsmas,(co,nk,kon,ipkon,lakon,ne,nodeboun,
526 ndirboun,xboun,nboun,
527 ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
528 nforc,nelemload,sideload,xload,nload,xbody,ipobody,
529 nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
530 nmethod,ikmpc,ilmpc,ikboun,ilboun,
531 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
532 ielmat,ielorien,norien,orab,ntmat_,
533 t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
534 nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
535 xstiff,npmat_,&dtime,matname,mi,
536 ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
537 physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
538 &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
539 xstateini,xstate,thicke,
540 integerglob,doubleglob,tieset,istartset,iendset,
541 ialset,ntie,&nasym,pslavsurf,pmastsurf,mortar,clearini,
542 ielprop,prop,&ne0,&kscale,iponoel,inoel,&network));*/
543 }
544 }
545
546 SFREE(fext);
547
548 if(*nmethod==0){
549
550 /* error occurred in mafill: storing the geometry in frd format */
551
552 ++*kode;
553 NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
554 if(strcmp1(&filab[1044],"ZZS")==0){
555 NNEW(neigh,ITG,40**ne);
556 NNEW(ipneigh,ITG,*nk);
557 }
558
559 frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
560 kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
561 nstate_,istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
562 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
563 mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
564 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
565 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,
566 ielprop,prop,sti);
567
568 if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
569 SFREE(inum);FORTRAN(stop,());
570
571 }
572
573 /* LU decomposition of the left hand matrix */
574
575 if(igreen==0){
576 if(nasym==1){sigma=0.;}else{sigma=1.;}
577 }else{
578 if(*nforc>0){sigma=xforc[0];}
579 }
580
581 if(*isolver==0){
582 #ifdef SPOOLES
583 spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
584 &symmetryflag,&inputformat,&nzs[2]);
585 #else
586 printf("*ERROR in arpack: the SPOOLES library is not linked\n\n");
587 FORTRAN(stop,());
588 #endif
589 }
590 else if(*isolver==4){
591 #ifdef SGI
592 token=1;
593 sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
594 #else
595 printf("*ERROR in arpack: the SGI library is not linked\n\n");
596 FORTRAN(stop,());
597 #endif
598 }
599 else if(*isolver==5){
600 #ifdef TAUCS
601 tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
602 #else
603 printf("*ERROR in arpack: the TAUCS library is not linked\n\n");
604 FORTRAN(stop,());
605 #endif
606 }
607 else if(*isolver==6){
608 #ifdef MATRIXSTORAGE
609 matrixstorage(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1],
610 ntrans,inotr,trab,co,nk,nactdof,jobnamec,mi,ipkon,
611 lakon,kon,ne,mei,nboun,nmpc,cs,mcs,ithermal,nmethod);
612 #else
613 printf("*ERROR in arpack: the MATRIXSTORAGE library is not linked\n\n");
614 FORTRAN(stop,());
615 #endif
616 }
617 else if(*isolver==7){
618 #ifdef PARDISO
619 pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
620 &symmetryflag,&inputformat,jq,&nzs[2]);
621 #else
622 printf("*ERROR in arpack: the PARDISO library is not linked\n\n");
623 FORTRAN(stop,());
624 #endif
625 }
626 else if(*isolver==8){
627 #ifdef PASTIX
628 pastix_factor_main(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
629 &symmetryflag,&inputformat,jq,&nzs[2]);
630 #else
631 printf("*ERROR in arpack: the PASTIX library is not linked\n\n");
632 FORTRAN(stop,());
633 #endif
634 }
635
636 if(igreen==0){
637
638 /* calculating the eigenvalues and eigenmodes */
639
640 printf(" Calculating the eigenvalues and the eigenmodes\n\n");
641
642 ido=0;
643 ldz=neq[1];
644 iparam[0]=1;
645 iparam[2]=mxiter;
646 iparam[3]=1;
647 iparam[6]=3;
648
649 info=0;
650
651 NNEW(resid,double,neq[1]);
652 NNEW(z,double,(long long)ncv*neq[1]);
653 NNEW(workd,double,3*neq[1]);
654
655 if(nasym==1){
656 lworkl=3*ncv*(2+ncv);
657 NNEW(workl,double,lworkl);
658 FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
659 workl,&lworkl,&info));
660 }else{
661 lworkl=ncv*(8+ncv);
662 NNEW(workl,double,lworkl);
663 FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
664 workl,&lworkl,&info));
665 }
666
667 NNEW(temp_array,double,neq[1]);
668
669 while((ido==-1)||(ido==1)||(ido==2)){
670 if(ido==-1){
671 if(nasym==1){
672 FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow,nzs));
673 }else{
674 FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow));
675 }
676 }
677 if((ido==-1)||(ido==1)){
678
679 /* solve the linear equation system */
680
681 if(ido==-1){
682 if(*isolver==0){
683 #ifdef SPOOLES
684 spooles_solve(temp_array,&neq[1]);
685 #endif
686 }
687 else if(*isolver==4){
688 #ifdef SGI
689 sgi_solve(temp_array,token);
690 #endif
691 }
692 else if(*isolver==5){
693 #ifdef TAUCS
694 tau_solve(temp_array,&neq[1]);
695 #endif
696 }
697 else if(*isolver==7){
698 #ifdef PARDISO
699 pardiso_solve(temp_array,&neq[1],&symmetryflag,&inputformat,&nrhs);
700 #endif
701 }
702 else if(*isolver==8){
703 #ifdef PASTIX
704 if( pastix_solve(temp_array,&neq[1],&symmetryflag,&nrhs)==-1 )
705 printf("*WARNING in arpack: solving step didn't converge! Continuing anyway\n");
706 #endif
707 }
708 for(jrow=0;jrow<neq[1];jrow++){
709 workd[ipntr[1]-1+jrow]=temp_array[jrow];
710 }
711 }
712 else if(ido==1){
713 if(*isolver==0){
714 #ifdef SPOOLES
715 spooles_solve(&workd[ipntr[2]-1],&neq[1]);
716 #endif
717 }
718 else if(*isolver==4){
719 #ifdef SGI
720 sgi_solve(&workd[ipntr[2]-1],token);
721 #endif
722 }
723 else if(*isolver==5){
724 #ifdef TAUCS
725 tau_solve(&workd[ipntr[2]-1],&neq[1]);
726 #endif
727 }
728 else if(*isolver==7){
729 #ifdef PARDISO
730 pardiso_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag,&inputformat,&nrhs);
731 #endif
732 }
733 else if(*isolver==8){
734 #ifdef PASTIX
735 if( pastix_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag,&nrhs)==-1 )
736 printf("*WARNING in arpack: solving step didn't converge! Continuing anyway\n");
737 #endif
738 }
739 for(jrow=0;jrow<neq[1];jrow++){
740 workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
741 }
742 }
743
744 }
745
746 if(ido==2){
747 if(nasym==1){
748 FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
749 adb,aub,jq,irow,nzs));
750 }else{
751 FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
752 adb,aub,jq,irow));
753 }
754 }
755
756 if(nasym==1){
757 FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
758 iparam,ipntr,workd,workl,&lworkl,&info));
759 }else{
760 FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
761 iparam,ipntr,workd,workl,&lworkl,&info));
762 }
763 }
764
765 if(info!=0){
766 printf("*ERROR in d[n,s]aupd: info=%" ITGFORMAT "\n",info);
767 printf(" # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
768 }
769
770 NNEW(select,ITG,ncv);
771
772 if(nasym==1){
773 NNEW(d,double,nev+1);
774 NNEW(di,double,nev+1);
775 NNEW(workev,double,3*ncv);
776 FORTRAN(dneupd,(&rvec,howmny,select,d,di,z,&ldz,&sigma,&sigmai,
777 workev,bmat,&neq[1],which,&nev,&tol,resid,
778 &ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
779 SFREE(workev);
780 NNEW(dc,double,2*nev);
781
782 /* storing as complex number and taking the square root */
783
784 for(j=0;j<nev;j++){
785 dc[2*j]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])+d[j])/sqrt(2.);
786 dc[2*j+1]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])-d[j])/sqrt(2.);
787 if(di[j]<0.) dc[2*j+1]=-dc[2*j+1];
788 }
789 FORTRAN(writeevcomplex,(dc,&nev,&fmin,&fmax));
790 SFREE(di);SFREE(dc);
791 }else{
792 NNEW(d,double,nev);
793 FORTRAN(dseupd,(&rvec,howmny,select,d,z,&ldz,&sigma,bmat,&neq[1],which,
794 &nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
795 FORTRAN(writeev,(d,&nev,&fmin,&fmax));
796 }
797 SFREE(select);SFREE(workd);SFREE(workl);SFREE(resid);
798
799 if(info!=0){
800 printf("*ERROR in d[n,s]eupd: info=%" ITGFORMAT "\n",info);
801 }
802
803 /* check the normalization of the eigenmodes */
804
805 for(j=0;j<nev;j++){
806 zmax=0.;
807 kmax=-1;
808 kref=j*neq[1];
809 if(nasym==1){
810 FORTRAN(opas,(&neq[1],&z[kref],temp_array,
811 adb,aub,jq,irow,nzs));
812 }else{
813 FORTRAN(op,(neq,&z[kref],temp_array,
814 adb,aub,jq,irow));
815 }
816 sum=0;
817 for(k=0;k<neq[1];k++){sum+=z[kref+k]*temp_array[k];}
818 printf("U^T*M*U=%f for eigenmode %d\n",j+1,sum);
819
820 /* check for the maximum */
821
822 /* for(k=0;k<neq[1];k++){
823 if(fabs(z[kref+k])>zmax){
824 zmax=fabs(z[kref+k]);
825 kmax=k;
826 }
827 }
828 zmax=z[kref+k];*/
829
830 /* normalizing the eigenmode (if not yet normalized by ARPACK */
831
832 /* if((fabs(sum-1.)>1.e-5)&&(sum>0.)){
833 printf("sum=%e\n");
834 printf("normalizing the eigenmode \n");
835 if(zmax>0.){
836 for(k=0;k<neq[1];k++){z[kref+k]/=sqrt(sum);}
837 }else{
838 for(k=0;k<neq[1];k++){z[kref+k]/=-sqrt(sum);}
839 }
840 }else if(zmax<0.){
841 for(k=0;k<neq[1];k++){z[kref+k]*=-1.;}
842 }*/
843
844 if(fabs(sum-1.)>1.e-5){
845 printf("sum=%e\n",sum);
846 printf("normalizing the eigenmode \n");
847 for(k=0;k<neq[1];k++){z[kref+k]/sqrt(sum);}
848 }
849 }
850
851 SFREE(temp_array);
852
853 }else{
854
855 /* Green functions */
856
857 /* storing omega_0^2 into d */
858
859 NNEW(d,double,*nforc);
860 for(j=0;j<*nforc;j++){d[j]=xforc[0];}
861 NNEW(z,double,neq[1]);
862
863 }
864
865
866 /* writing the eigenvalues and mass matrix to a binary file */
867
868 if(mei[3]==1){
869
870 strcpy(fneig,jobnamec);
871 strcat(fneig,".eig");
872
873 if((f1=fopen(fneig,"wb"))==NULL){
874 printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
875
876 exit(0);
877 }
878
879 /* storing a zero as indication that this was not a
880 cyclic symmetry calculation */
881
882 if(fwrite(&zero,sizeof(ITG),1,f1)!=1){
883 printf("*ERROR saving the cyclic symmetry flag to the eigenvalue file...");
884 exit(0);
885 }
886
887 /* Hermitian */
888
889 if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
890 printf("*ERROR saving the Hermitian flag to the eigenvalue file...");
891 exit(0);
892 }
893
894 /* perturbation parameter iperturb[0] */
895
896 if(fwrite(&iperturb[0],sizeof(ITG),1,f1)!=1){
897 printf("*ERROR saving the perturbation flag to the eigenvalue file...");
898 exit(0);
899 }
900
901 /* reference displacements */
902
903 if(iperturb[0]==1){
904 if(fwrite(vold,sizeof(double),mt**nk,f1)!=mt**nk){
905 printf("*ERROR saving the reference displacements to the eigenvalue file...");
906 exit(0);
907 }
908 }
909
910 /* storing the number of eigenvalues */
911
912 if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
913 printf("*ERROR saving the number of eigenvalues to the eigenvalue file...");
914 exit(0);
915 }
916
917 /* the eigenfrequencies are stored as radians/time */
918
919 if(fwrite(d,sizeof(double),nev,f1)!=nev){
920 printf("*ERROR saving the eigenfrequencies to the eigenvalue file...");
921 exit(0);
922 }
923
924 /* storing the stiffness matrix */
925
926 if(fwrite(ad,sizeof(double),neq[1],f1)!=neq[1]){
927 printf("*ERROR saving the diagonal of the stiffness matrix to the eigenvalue file...");
928 exit(0);
929 }
930 if(fwrite(au,sizeof(double),nzs[2],f1)!=nzs[2]){
931 printf("*ERROR saving the off-diagonal terms of the stiffness matrix to the eigenvalue file...");
932 exit(0);
933 }
934
935 /* storing the mass matrix */
936
937 if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
938 printf("*ERROR saving the diagonal of the mass matrix to the eigenvalue file...");
939 exit(0);
940 }
941 if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
942 printf("*ERROR saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
943 exit(0);
944 }
945 }
946
947 if(igreen==0){
948
949 /* calculating the participation factors and the relative effective
950 modal mass */
951
952 if(*ithermal!=2){
953 FORTRAN(effectivemodalmass,(neq,nactdof,mi,adb,aub,jq,irow,&nev,z,co,nk));
954 }
955 }
956
957 SFREE(au);SFREE(ad);
958
959 /* calculating the displacements and the stresses and storing */
960 /* the results in frd format for each valid eigenmode */
961
962 NNEW(v,double,mt**nk);
963 NNEW(fn,double,mt**nk);
964 NNEW(stn,double,6**nk);
965 NNEW(inum,ITG,*nk);
966 NNEW(stx,double,6*mi[0]**ne);
967 NNEW(eme,double,6*mi[0]**ne);
968 if(*ithermal>1){
969 NNEW(qfn,double,3**nk);
970 NNEW(qfx,double,3*mi[0]**ne);
971 }
972
973 if(strcmp1(&filab[261],"E ")==0) NNEW(een,double,6**nk);
974 if(strcmp1(&filab[2697],"ME ")==0) NNEW(emn,double,6**nk);
975 if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
976 if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
977
978 // NNEW(temp_array,double,neq[1]);
979
980 lfin=0;
981 for(j=0;j<nev;++j){
982 lint=lfin;
983
984 if(igreen==0){
985 lfin=lfin+neq[1];
986
987 /* sum=0.;
988 for(k=0;k<neq[1];++k)
989 temp_array[k]=0.;
990 if(nasym==1){
991 FORTRAN(opas,(&neq[1],&z[lint],temp_array,adb,aub,jq,irow,nzs));
992 }else{
993 FORTRAN(op,(&neq[1],&z[lint],temp_array,adb,aub,jq,irow));
994 }
995 for(k=0;k<neq[1];++k)
996 sum+=z[lint+k]*temp_array[k];
997 for(k=0;k<neq[1];++k)
998 z[lint+k]=z[lint+k]/sqrt(sum);*/
999 }else{
1000
1001 /* Green function: solve the system ([K]-w**2*[M]){X}={Y} */
1002
1003 node=nodeforc[2*j];
1004 idir=ndirforc[j];
1005
1006 /* check whether degree of freedom is active */
1007
1008 if(nactdof[mt*(node-1)+idir]==0){
1009 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);
1010 FORTRAN(stop,());
1011 }
1012
1013 /* defining a unit force on the rhs */
1014
1015 DMEMSET(z,0,neq[1],0.);
1016 z[nactdof[mt*(node-1)+idir]-1]=1.;
1017
1018 /* solving the system */
1019
1020 if(neq[1]>0){
1021 if(*isolver==0){
1022 #ifdef SPOOLES
1023 spooles_solve(z,&neq[1]);
1024 #endif
1025 }
1026 else if(*isolver==7){
1027 #ifdef PARDISO
1028 pardiso_solve(z,&neq[1],&symmetryflag,&inputformat,&nrhs);
1029 #endif
1030 }
1031 else if(*isolver==8){
1032 #ifdef PASTIX
1033 if( pastix_solve(z,&neq[1],&symmetryflag,&nrhs)==-1 )
1034 printf("*WARNING in arpack: solving step didn't converge! Continuing anyway\n");
1035 #endif
1036
1037 }
1038 }
1039 }
1040
1041 if(mei[3]==1){
1042 if(fwrite(&z[lint],sizeof(double),neq[1],f1)!=neq[1]){
1043 printf("*ERROR saving data to the eigenvalue file...");
1044 exit(0);
1045 }
1046 }
1047
1048 if(igreen==0){
1049
1050 /* check whether the frequency belongs to the requested
1051 interval */
1052
1053 if(fmin>-0.5){
1054 if(fmin*fmin>d[j]) continue;
1055 }
1056 if(fmax>-0.5){
1057 if(fmax*fmax<d[j]) continue;
1058 }
1059 }
1060
1061 if(*nprint>0) FORTRAN(writehe,(&j));
1062
1063 NNEW(eei,double,6*mi[0]*ne0);
1064 if(*nener==1){
1065 NNEW(stiini,double,6*mi[0]*ne0);
1066 NNEW(emeini,double,6*mi[0]*ne0);
1067 NNEW(enerini,double,mi[0]*ne0);}
1068
1069 // memset(&v[0],0.,sizeof(double)*mt**nk);
1070 DMEMSET(v,0,mt**nk,0.);
1071 if(*iperturb==0){
1072 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1073 stx,elcon,
1074 nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1075 norien,orab,ntmat_,t0,t0,ithermal,
1076 prestr,iprestr,filab,eme,emn,een,iperturb,
1077 f,fn,nactdof,&iout,qa,vold,&z[lint],
1078 nodeboun,ndirboun,xboun,nboun,ipompc,
1079 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1080 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1081 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1082 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
1083 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
1084 iendset,
1085 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1086 nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
1087 &reltime,&ne0,thicke,shcon,nshcon,
1088 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1089 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1090 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
1091 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
1092 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
1093 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
1094 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
1095 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
1096 &intscheme);}
1097 else{
1098 results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
1099 stx,elcon,
1100 nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
1101 norien,orab,ntmat_,t0,t1old,ithermal,
1102 prestr,iprestr,filab,eme,emn,een,iperturb,
1103 f,fn,nactdof,&iout,qa,vold,&z[lint],
1104 nodeboun,ndirboun,xboun,nboun,ipompc,
1105 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1106 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
1107 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
1108 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1109 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1110 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1111 nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
1112 &ne0,thicke,shcon,nshcon,
1113 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
1114 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1115 islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
1116 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
1117 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
1118 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
1119 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
1120 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
1121 &intscheme);
1122 }
1123 SFREE(eei);
1124 if(*nener==1){
1125 SFREE(stiini);SFREE(emeini);SFREE(enerini);}
1126
1127 ++*kode;
1128 if(igreen==0){
1129 if(d[j]>=0.){
1130 freq=sqrt(d[j])/6.283185308;
1131 }else{
1132 freq=0.;
1133 }
1134 }else{
1135 freq=*ttime+j;
1136 }
1137
1138 if(strcmp1(&filab[1044],"ZZS")==0){
1139 NNEW(neigh,ITG,40**ne);
1140 NNEW(ipneigh,ITG,*nk);
1141 }
1142
1143 if(*iperturb==0){
1144 frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
1145 kode,filab,een,t1,fn,&freq,epn,ielmat,matname,enern,xstaten,
1146 nstate_,istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
1147 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1148 mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1149 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1150 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,
1151 ielprop,prop,sti);
1152 }else{
1153 frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
1154 kode,filab,een,t1old,fn,&freq,epn,ielmat,matname,enern,xstaten,
1155 nstate_,istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
1156 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
1157 mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
1158 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
1159 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,
1160 ielprop,prop,sti);
1161 }
1162
1163 if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
1164
1165 if(*iperturb!=0){
1166 for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
1167 xstate[k]=xstateini[k];
1168 }
1169 }
1170 }
1171
1172 /*--------------------------------------------------------------------*/
1173 /*
1174 -----------
1175 free memory
1176 -----------
1177 */
1178 if(*isolver==0){
1179 #ifdef SPOOLES
1180 spooles_cleanup();
1181 #endif
1182 }
1183 else if(*isolver==4){
1184 #ifdef SGI
1185 sgi_cleanup(token);
1186 #endif
1187 }
1188 else if(*isolver==5){
1189 #ifdef TAUCS
1190 tau_cleanup();
1191 #endif
1192 }
1193 else if(*isolver==7){
1194 #ifdef PARDISO
1195 pardiso_cleanup(&neq[1],&symmetryflag,&inputformat);
1196 #endif
1197 }
1198 else if(*isolver==8){
1199 #ifdef PASTIX
1200 #endif
1201 }
1202
1203 if(igreen==0){
1204 if((fmax>-0.5)&&(fmax*fmax>d[nev-1])){
1205 printf("\n*WARNING: not all frequencies in the requested interval might be found;\nincrease the number of requested frequencies\n");
1206 }
1207 }
1208
1209 if(mei[3]==1){
1210 fclose(f1);
1211 }
1212
1213 if(*iperturb!=0){
1214 if(ncont!=0){
1215 *ne=ne0;
1216 if(*nener==1){
1217 RENEW(ener,double,mi[0]**ne*2);
1218 }
1219 RENEW(ipkon,ITG,*ne);
1220 RENEW(lakon,char,8**ne);
1221 RENEW(kon,ITG,*nkon);
1222 if(*norien>0){
1223 RENEW(ielorien,ITG,mi[2]**ne);
1224 }
1225 RENEW(ielmat,ITG,mi[2]**ne);
1226 SFREE(cg);SFREE(straight);
1227 SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
1228 SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
1229 SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
1230 SFREE(springarea);SFREE(xmastnor);
1231
1232 if(*mortar==0){
1233 SFREE(areaslav);
1234 }else if(*mortar==1){
1235 SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
1236 SFREE(islavact);
1237 }
1238 }
1239 }
1240
1241 SFREE(adb);SFREE(aub);
1242 //SFREE(temp_array);
1243
1244 SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);SFREE(eme);
1245 SFREE(z);SFREE(d);SFREE(xstiff);
1246
1247 //SFREE(ipobody);
1248
1249 if(*ithermal>1){SFREE(qfn);SFREE(qfx);}
1250
1251 if(*nstate_!=0){SFREE(xstateini);}
1252
1253 if(strcmp1(&filab[261],"E ")==0) SFREE(een);
1254 if(strcmp1(&filab[2697],"ME ")==0) SFREE(emn);
1255 if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1256 if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
1257
1258 if(*iperturb!=0){
1259 mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
1260 mpcinfo[3]=maxlenmpc;
1261 }
1262
1263 if(igreen==1){
1264 *nmethod=13;
1265 }
1266
1267 *irowp=irow;*enerp=ener;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
1268 *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;
1269
1270 *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
1271
1272 return;
1273 }
1274
1275 #endif
1276