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