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 <string.h>
22 #include "CalculiX.h"
23 #ifdef SPOOLES
24 #include "spooles.h"
25 #endif
26 #ifdef SGI
27 #include "sgi.h"
28 #endif
29 #ifdef TAUCS
30 #include "tau.h"
31 #endif
32
expand(double * co,ITG * nk,ITG * kon,ITG * ipkon,char * lakon,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 * neq,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 * ielmat,ITG * ielorien,ITG * norien,double * orab,ITG * ntmat_,double * t0,ITG * ithermal,double * prestr,ITG * iprestr,double * vold,ITG * iperturb,double * sti,ITG * nzs,double * adb,double * aub,char * filab,double * eme,double * plicon,ITG * nplicon,double * plkcon,ITG * nplkcon,double * xstate,ITG * npmat_,char * matname,ITG * mi,ITG * ics,double * cs,ITG * mpcend,ITG * ncmat_,ITG * nstate_,ITG * mcs,ITG * nkon,double * ener,char * jobnamec,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 * ttime,double * fmpc,ITG * nev,double ** zp,ITG * iamboun,double * xbounold,ITG * nsectors,ITG * nm,ITG * icol,ITG * irow,ITG * nzl,ITG * nam,ITG * ipompcold,ITG * nodempcold,double * coefmpcold,char * labmpcold,ITG * nmpcold,double * xloadold,ITG * iamload,double * t1old,double * t1,ITG * iamt1,double * xstiff,ITG ** icolep,ITG ** jqep,ITG ** irowep,ITG * isolver,ITG * nzse,double ** adbep,double ** aubep,ITG * iexpl,ITG * ibody,double * xbody,ITG * nbody,double * cocon,ITG * ncocon,char * tieset,ITG * ntie,ITG * imddof,ITG * nmddof,ITG * imdnode,ITG * nmdnode,ITG * imdboun,ITG * nmdboun,ITG * imdmpc,ITG * nmdmpc,ITG ** izdofp,ITG * nzdof,ITG * nherm,double * xmr,double * xmi,char * typeboun,ITG * ielprop,double * prop,char * orname,ITG * itiefac,double * t0g,double * t1g)33 void expand(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon,
34 ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun,
35 ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
36 ITG *nmpc, ITG *nodeforc, ITG *ndirforc,double *xforc,
37 ITG *nforc, ITG *nelemload, char *sideload, double *xload,
38 ITG *nload, ITG *nactdof, ITG *neq,
39 ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun,
40 double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
41 double *alcon, ITG *nalcon, double *alzero, ITG *ielmat,
42 ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_,
43 double *t0,ITG *ithermal,double *prestr, ITG *iprestr,
44 double *vold,ITG *iperturb, double *sti, ITG *nzs,
45 double *adb, double *aub,char *filab, double *eme,
46 double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
47 double *xstate, ITG *npmat_, char *matname, ITG *mi,
48 ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_,
49 ITG *nstate_, ITG *mcs, ITG *nkon, double *ener,
50 char *jobnamec, char *output, char *set, ITG *nset,ITG *istartset,
51 ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
52 char *prset, ITG *nener, double *trab,
53 ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
54 ITG *nev, double **zp, ITG *iamboun, double *xbounold,
55 ITG *nsectors, ITG *nm,ITG *icol,ITG *irow,ITG *nzl, ITG *nam,
56 ITG *ipompcold, ITG *nodempcold, double *coefmpcold,
57 char *labmpcold, ITG *nmpcold, double *xloadold, ITG *iamload,
58 double *t1old,double *t1,ITG *iamt1, double *xstiff,ITG **icolep,
59 ITG **jqep,ITG **irowep,ITG *isolver,
60 ITG *nzse,double **adbep,double **aubep,ITG *iexpl,ITG *ibody,
61 double *xbody,ITG *nbody,double *cocon,ITG *ncocon,
62 char* tieset,ITG* ntie,ITG *imddof,ITG *nmddof,
63 ITG *imdnode,ITG *nmdnode,ITG *imdboun,ITG *nmdboun,
64 ITG *imdmpc,ITG *nmdmpc, ITG **izdofp, ITG *nzdof,ITG *nherm,
65 double *xmr,double *xmi,char *typeboun,ITG *ielprop,double *prop,
66 char *orname,ITG *itiefac,double *t0g,double *t1g){
67
68 /* calls the Arnoldi Package (ARPACK) for cyclic symmetry calculations */
69
70 char *filabt,lakonl[2]=" \0",*labmpc2=NULL;
71
72 ITG *inum=NULL,k,idir,j,iout=0,index,inode,id,i,idof,im,
73 ielas,icmd,kk,l,nkt,icntrl,imag=1,icomplex,kkv,kk6,iterm,
74 lprev,ilength,ij,i1,i2,iel,ielset,node,indexe,nope,ml1,nelem,
75 *inocs=NULL,*ielcs=NULL,jj,l1,l2,is,nlabel,*nshcon=NULL,
76 nodeleft,*noderight=NULL,numnodes,ileft,kflag=2,itr,locdir,
77 neqh,j1,nodenew,mt=mi[1]+1,istep=1,iinc=1,iit=-1,ne0,
78 network=0,noderight_,*izdof=*izdofp,iload,iforc,*iznode=NULL,
79 nznode,ll,icfd=0,*inomat=NULL,mortar=0,*islavact=NULL,*ipobody=NULL,
80 *islavnode=NULL,*nslavnode=NULL,*islavsurf=NULL,idirnew,
81 *iponoel=NULL,*inoel=NULL,mscalmethod=0,intscheme=0,
82 *islavelinv=NULL,*irowtloc=NULL,*jqtloc=NULL,nboun2,
83 *ndirboun2=NULL,*nodeboun2=NULL,nmpc2,*ipompc2=NULL,*nodempc2=NULL,
84 *ikboun2=NULL,*ilboun2=NULL,*ikmpc2=NULL,*ilmpc2=NULL,mortartrafoflag=0;
85
86 long long lint;
87
88 double *stn=NULL,*v=NULL,*temp_array=NULL,*vini=NULL,*csmass=NULL,
89 *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[4],*epn=NULL,summass,
90 *stiini=NULL,*emn=NULL,*emeini=NULL,*clearini=NULL,
91 *xstateini=NULL,theta,pi,*coefmpcnew=NULL,t[3],ctl,stl,
92 *stx=NULL,*enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,
93 *qfx=NULL,*qfn=NULL,xreal,ximag,*vt=NULL,sum,*voldt=NULL,
94 *coefright=NULL,coef,a[9],ratio,reltime,
95 *shcon=NULL,*springarea=NULL,*z=*zp, *zdof=NULL, *thicke=NULL,
96 *sumi=NULL,
97 *vti=NULL,*pslavsurf=NULL,*pmastsurf=NULL,*cdn=NULL,
98 *energyini=NULL,*energy=NULL,*smscale=NULL,
99 *autloc=NULL,*xboun2=NULL,*coefmpc2=NULL;
100
101 /* dummy arguments for the results call */
102
103 double *veold=NULL,*accold=NULL,bet,gam,dtime,time;
104
105 pi=4.*atan(1.);
106 neqh=neq[1]/2;
107
108 ne0=*ne;
109
110 noderight_=10;
111 NNEW(noderight,ITG,noderight_);
112 NNEW(coefright,double,noderight_);
113
114 NNEW(v,double,2*mt**nk);
115 NNEW(vt,double,mt**nk**nsectors);
116
117 NNEW(fn,double,2*mt**nk);
118 NNEW(stn,double,12**nk);
119 NNEW(inum,ITG,*nk);
120 NNEW(stx,double,6*mi[0]**ne);
121
122 nlabel=55;
123 NNEW(filabt,char,87*nlabel);
124 for(i=1;i<87*nlabel;i++) filabt[i]=' ';
125 filabt[0]='U';
126
127 NNEW(temp_array,double,neq[1]);
128 NNEW(coefmpcnew,double,*mpcend);
129
130 nkt=*nsectors**nk;
131
132 /* assigning nodes and elements to sectors */
133
134 NNEW(inocs,ITG,*nk);
135 NNEW(ielcs,ITG,*ne);
136 ielset=cs[12];
137 if((*mcs!=1)||(ielset!=0)){
138 for(i=0;i<*nk;i++) inocs[i]=-1;
139 for(i=0;i<*ne;i++) ielcs[i]=-1;
140 }
141 NNEW(csmass,double,*mcs);
142 if(*mcs==1) csmass[0]=1.;
143
144 for(i=0;i<*mcs;i++){
145 is=cs[17*i];
146 // if(is==1) continue;
147 ielset=cs[17*i+12];
148 if(ielset==0) continue;
149 for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
150 if(ialset[i1]>0){
151 iel=ialset[i1]-1;
152 if(ipkon[iel]<0) continue;
153 ielcs[iel]=i;
154 indexe=ipkon[iel];
155 if(*mcs==1){
156 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
157 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
158 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
159 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
160 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
161 else if (strcmp1(&lakon[8*iel+3],"6")==0)nope=6;
162 else if (strcmp1(&lakon[8*iel],"ES")==0){
163 lakonl[0]=lakon[8*iel+7];
164 nope=atoi(lakonl)+1;}
165 else continue;
166 }else{
167 nelem=iel+1;
168 FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
169 ielmat,&nope,t0,rhcon,nrhcon,ntmat_,
170 ithermal,&csmass[i],ielprop,prop));
171 }
172 for(i2=0;i2<nope;++i2){
173 node=kon[indexe+i2]-1;
174 inocs[node]=i;
175 }
176 }
177 else{
178 iel=ialset[i1-2]-1;
179 do{
180 iel=iel-ialset[i1];
181 if(iel>=ialset[i1-1]-1) break;
182 if(ipkon[iel]<0) continue;
183 ielcs[iel]=i;
184 indexe=ipkon[iel];
185 if(*mcs==1){
186 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
187 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
188 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
189 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
190 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
191 else {nope=6;}
192 }else{
193 nelem=iel+1;
194 FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
195 ielmat,&nope,t0,rhcon,nrhcon,ntmat_,
196 ithermal,&csmass[i],ielprop,prop));
197 }
198 for(i2=0;i2<nope;++i2){
199 node=kon[indexe+i2]-1;
200 inocs[node]=i;
201 }
202 }while(1);
203 }
204 }
205 // printf("expand.c mass = %" ITGFORMAT ",%e\n",i,csmass[i]);
206 }
207
208 /* copying imdnode into iznode
209 iznode contains the nodes in which output is requested and
210 the nodes in which loading is applied */
211
212 NNEW(iznode,ITG,*nk);
213 for(j=0;j<*nmdnode;j++){iznode[j]=imdnode[j];}
214 nznode=*nmdnode;
215
216 /* expanding imddof, imdnode, imdboun and imdmpc */
217
218 for(i=1;i<*nsectors;i++){
219 for(j=0;j<*nmddof;j++){
220 imddof[i**nmddof+j]=imddof[j]+i*neqh;
221 }
222 for(j=0;j<*nmdnode;j++){
223 imdnode[i**nmdnode+j]=imdnode[j]+i**nk;
224 }
225 for(j=0;j<*nmdboun;j++){
226 imdboun[i**nmdboun+j]=imdboun[j]+i**nboun;
227 }
228 for(j=0;j<*nmdmpc;j++){
229 imdmpc[i**nmdmpc+j]=imdmpc[j]+i**nmpc;
230 }
231 }
232 (*nmddof)*=(*nsectors);
233 (*nmdnode)*=(*nsectors);
234 (*nmdboun)*=(*nsectors);
235 (*nmdmpc)*=(*nsectors);
236
237 /* creating a field with the degrees of freedom in which the eigenmodes
238 are needed:
239 1. all dofs in which the solution is needed (=imddof)
240 2. all dofs in which loading was applied
241 */
242
243 NNEW(izdof,ITG,neqh**nsectors);
244 for(j=0;j<*nmddof;j++){izdof[j]=imddof[j];}
245 *nzdof=*nmddof;
246
247 /* generating the coordinates for the other sectors */
248
249 icntrl=1;
250
251 FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,emn));
252
253 for(jj=0;jj<*mcs;jj++){
254 is=(ITG)(cs[17*jj]+0.5);
255 for(i=1;i<is;i++){
256
257 theta=i*2.*pi/cs[17*jj];
258
259 for(l=0;l<*nk;l++){
260 if(inocs[l]==jj){
261 co[3*l+i*3**nk]=co[3*l];
262 co[1+3*l+i*3**nk]=co[1+3*l]+theta;
263 co[2+3*l+i*3**nk]=co[2+3*l];
264 if(*ntrans>0) inotr[2*l+i*2**nk]=inotr[2*l];
265 }
266 }
267 for(l=0;l<*nkon;l++){kon[l+i**nkon]=kon[l]+i**nk;}
268 for(l=0;l<*ne;l++){
269 if(ielcs[l]==jj){
270 if(ipkon[l]>=0){
271 ipkon[l+i**ne]=ipkon[l]+i**nkon;
272 ielmat[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
273 if(*norien>0) ielorien[l+i**ne]=ielorien[l];
274 for(l1=0;l1<8;l1++){
275 l2=8*l+l1;
276 lakon[l2+i*8**ne]=lakon[l2];
277 }
278 }else{
279 ipkon[l+i**ne]=ipkon[l];
280 }
281 }
282 }
283 }
284 }
285
286 icntrl=-1;
287
288 FORTRAN(rectcyl,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,&imag,mi,emn));
289
290 /* expand nactdof */
291
292 for(i=1;i<*nsectors;i++){
293 lint=i*mt**nk;
294 for(j=0;j<mt**nk;j++){
295 if(nactdof[j]>0){
296 nactdof[lint+j]=nactdof[j]+i*neqh;
297 }else{
298 nactdof[lint+j]=0;
299 }
300 }
301 }
302
303 /* copying the boundary conditions
304 (SPC's must be defined in cylindrical coordinates) */
305
306 for(i=1;i<*nsectors;i++){
307 for(j=0;j<*nboun;j++){
308 nodeboun[i**nboun+j]=nodeboun[j]+i**nk;
309 ndirboun[i**nboun+j]=ndirboun[j];
310 xboun[i**nboun+j]=xboun[j];
311 xbounold[i**nboun+j]=xbounold[j];
312 if(*nam>0) iamboun[i**nboun+j]=iamboun[j];
313 ikboun[i**nboun+j]=ikboun[j]+8*i**nk;
314 ilboun[i**nboun+j]=ilboun[j]+i**nboun;
315 }
316 }
317
318 /* distributed loads */
319
320 for(i=0;i<*nload;i++){
321 if(nelemload[2*i+1]<*nsectors){
322 nelemload[2*i]+=*ne*nelemload[2*i+1];
323 }else{
324 nelemload[2*i]+=*ne*(nelemload[2*i+1]-(*nsectors));
325 }
326 iload=i+1;
327 FORTRAN(addizdofdload,(nelemload,sideload,ipkon,kon,lakon,
328 nactdof,izdof,nzdof,mi,&iload,iznode,&nznode,nk,
329 imdnode,nmdnode));
330 }
331
332 /* body loads */
333
334 if(*nbody>0){
335 printf("\n *WARNING in expand: body loads are not allowed for modal dynamics\n and steady state dynamics calculations in cyclic symmetric structures;\n the load is removed\n\n");
336 *nbody=0;
337 }
338
339 /* sorting the elements with distributed loads */
340
341 if(*nload>0){
342 if(*nam>0){
343 FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
344 }else{
345 FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
346 }
347 }
348
349 /* point loads */
350
351 i=0;
352 while(i<*nforc){
353 node=nodeforc[2*i];
354
355 /* checking for a cylindrical transformation;
356 comparison with the cyclic symmetry system */
357
358 /* carthesian coordinate system */
359
360 if(nodeforc[2*i+1]<*nsectors){
361 nodeforc[2*i]+=*nk*nodeforc[2*i+1];
362 }else{
363 nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
364 }
365 i++;iforc=i;
366 FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
367 nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc,
368 ntrans,inotr));
369 }
370
371 /* loop over all eigenvalues; the loop starts from the highest eigenvalue
372 so that the reuse of z is not a problem
373 z before: real and imaginary part for a segment for all eigenvalues
374 z after: real part for all segments for all eigenvalues */
375
376 if(*nherm==1){
377 NNEW(zdof,double,(long long)*nev**nzdof);
378 }else{
379 NNEW(zdof,double,(long long)2**nev**nzdof);
380 NNEW(sumi,double,*nev);
381 }
382
383 for(j=*nev-1;j>-1;--j){
384 lint=(long long)2*j*neqh;
385
386 /* calculating the cosine and sine of the phase angle */
387
388 for(jj=0;jj<*mcs;jj++){
389 theta=nm[j]*2.*pi/cs[17*jj];
390 cs[17*jj+14]=cos(theta);
391 cs[17*jj+15]=sin(theta);
392 }
393
394 /* generating the cyclic MPC's (needed for nodal diameters
395 different from 0 */
396
397 NNEW(eei,double,6*mi[0]**ne);
398
399 DMEMSET(v,0,2*mt**nk,0.);
400
401 for(k=0;k<2*neqh;k+=neqh){
402
403 for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
404
405 if(k==0) {kk=0;kkv=0;kk6=0;}
406 else {kk=*nk;kkv=mt**nk;kk6=6**nk;}
407 for(i=0;i<*nmpc;i++){
408 index=ipompc[i]-1;
409 /* check whether thermal mpc */
410 if(nodempc[3*index+1]==0) continue;
411 coefmpcnew[index]=coefmpc[index];
412 while(1){
413 index=nodempc[3*index+2];
414 if(index==0) break;
415 index--;
416
417 icomplex=0;
418 inode=nodempc[3*index];
419 if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
420 icomplex=atoi(&labmpc[20*i+6]);}
421 else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
422 for(ij=0;ij<*mcs;ij++){
423 lprev=cs[ij*17+13];
424 ilength=cs[ij*17+3];
425 FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
426 if(id!=0){
427 if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
428 }
429 }
430 }
431
432 if(icomplex!=0){
433 idir=nodempc[3*index+1];
434 idof=nactdof[mt*(inode-1)+idir]-1;
435 if(idof<=-1){xreal=1.;ximag=1.;}
436 else{xreal=z[lint+idof];ximag=z[lint+idof+neqh];}
437 if(k==0) {
438 if(fabs(xreal)<1.e-30)xreal=1.e-30;
439 coefmpcnew[index]=coefmpc[index]*
440 (cs[17*(icomplex-1)+14]+
441 ximag/xreal*cs[17*(icomplex-1)+15]);}
442 else {
443 if(fabs(ximag)<1.e-30)ximag=1.e-30;
444 coefmpcnew[index]=coefmpc[index]*
445 (cs[17*(icomplex-1)+14]-
446 xreal/ximag*cs[17*(icomplex-1)+15]);}
447 }
448 else{coefmpcnew[index]=coefmpc[index];}
449 }
450 }
451
452 results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
453 stx,elcon,
454 nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
455 norien,orab,ntmat_,t0,t0,ithermal,
456 prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
457 f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
458 nodeboun,ndirboun,xboun,nboun,ipompc,
459 nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neqh,veold,accold,
460 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
461 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
462 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
463 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
464 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
465 nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
466 &ne0,thicke,shcon,nshcon,
467 sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
468 &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
469 islavsurf,ielprop,prop,energyini,energy,&iit,iponoel,
470 inoel,nener,orname,&network,ipobody,xbody,ibody,typeboun,
471 itiefac,tieset,smscale,&mscalmethod,nbody,t0g,t1g,
472 islavelinv,autloc,irowtloc,jqtloc,&nboun2,
473 ndirboun2,nodeboun2,xboun2,&nmpc2,ipompc2,nodempc2,coefmpc2,
474 labmpc2,ikboun2,ilboun2,ikmpc2,ilmpc2,&mortartrafoflag,
475 &intscheme);
476
477 }
478 SFREE(eei);
479
480 /* mapping the results to the other sectors */
481
482 if(*nherm!=1)NNEW(vti,double,mt**nk**nsectors);
483
484 icntrl=2;imag=1;
485
486 FORTRAN(rectcylexp,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,
487 iznode,&nznode,nsectors,nk,emn));
488
489 /* basis sector */
490
491 for(ll=0;ll<nznode;ll++){
492 l1=iznode[ll]-1;
493 for(l2=0;l2<mt;l2++){
494 l=mt*l1+l2;
495 vt[l]=v[l];
496 if(*nherm!=1)vti[l]=v[l+mt**nk];
497 }
498 }
499
500 /* other sectors */
501
502 for(jj=0;jj<*mcs;jj++){
503 ilength=cs[17*jj+3];
504 lprev=cs[17*jj+13];
505 for(i=1;i<*nsectors;i++){
506
507 theta=i*nm[j]*2.*pi/cs[17*jj];
508 ctl=cos(theta);
509 stl=sin(theta);
510
511 for(ll=0;ll<nznode;ll++){
512 l1=iznode[ll]-1;
513 if(inocs[l1]==jj){
514
515 /* check whether node lies on axis */
516
517 ml1=-l1-1;
518 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
519 if(id!=0){
520 if(ics[lprev+id-1]==ml1){
521 for(l2=0;l2<mt;l2++){
522 l=mt*l1+l2;
523 vt[l+mt**nk*i]=v[l];
524 if(*nherm!=1)vti[l+mt**nk*i]=v[l+mt**nk];
525 }
526 continue;
527 }
528 }
529 for(l2=0;l2<mt;l2++){
530 l=mt*l1+l2;
531 vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
532 if(*nherm!=1)vti[l+mt**nk*i]=stl*v[l]+ctl*v[l+mt**nk];
533 }
534 }
535 }
536 }
537 }
538
539 icntrl=-2;imag=0;
540
541 FORTRAN(rectcylexp,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,
542 &imag,mi,iznode,&nznode,nsectors,nk,emn));
543
544 /* storing the displacements into the expanded eigenvectors */
545
546 for(ll=0;ll<nznode;ll++){
547 i=iznode[ll]-1;
548 for(j1=0;j1<mt;j1++){
549
550 for(k=0;k<*nsectors;k++){
551 /* C-convention */
552 idof=nactdof[mt*(k**nk+i)+j1]-1;
553 if(idof>-1){
554 FORTRAN(nident,(izdof,&idof,nzdof,&id));
555 if(id!=0){
556 if(izdof[id-1]==idof){
557 if(*nherm==1){
558 zdof[(long long)j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
559 }else{
560 zdof[(long long)2*j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
561 zdof[(long long)(2*j+1)**nzdof+id-1]=vti[k*mt**nk+mt*i+j1];
562 }
563 }
564 }
565 }
566 }
567 }
568 }
569
570 if(*nherm!=1) SFREE(vti);
571
572 /* normalizing the eigenvectors with the mass */
573
574 sum=0.;
575 summass=0.;
576 for(i=0;i<*mcs;i++){
577 if (nm[j]==0||(nm[j]==(ITG)((cs[17*i]/2))&&(fmod(cs[17*i],2.)==0.))){
578 sum+=cs[17*i]*csmass[i];
579 }else{
580 sum+=cs[17*i]*csmass[i]/2.;
581 }
582 summass+=csmass[i];
583 }
584 if(fabs(summass)>1.e-20){
585 sum=sqrt(sum/summass);
586 }else{
587 printf("*ERROR in expand.c: total mass of structure is zero\n");
588 printf(" maybe no element sets were specified for the\n");
589 printf(" cyclic symmetry ties\n");
590 FORTRAN(stop,());
591 }
592
593 if(*nherm==1){
594 for(i=0;i<*nzdof;i++){zdof[(long long)j**nzdof+i]/=sum;}
595 }else{
596 for(i=0;i<*nzdof;i++){zdof[(long long)(2*j)**nzdof+i]/=sum;}
597 for(i=0;i<*nzdof;i++){zdof[(long long)(2*j+1)**nzdof+i]/=sum;}
598 sumi[j]=sqrt(sum);
599 }
600 }
601
602 /* mapping the reference displacements to the other sectors
603 expanding vold */
604
605 if(iperturb[0]==1){
606 NNEW(voldt,double,mt**nk**nsectors);
607
608 icntrl=2;
609
610 FORTRAN(rectcylvold,(co,vold,cs,&icntrl,mi,
611 iznode,&nznode,nsectors,nk));
612
613 /* basis sector */
614
615 for(ll=0;ll<nznode;ll++){
616 l1=iznode[ll]-1;
617 for(l2=0;l2<mt;l2++){
618 l=mt*l1+l2;
619 voldt[l]=vold[l];
620 }
621 }
622
623 /* other sectors */
624
625 for(jj=0;jj<*mcs;jj++){
626 ilength=cs[17*jj+3];
627 lprev=cs[17*jj+13];
628 for(i=1;i<*nsectors;i++){
629
630 for(ll=0;ll<nznode;ll++){
631 l1=iznode[ll]-1;
632 if(inocs[l1]==jj){
633
634 /* check whether node lies on axis */
635
636 ml1=-l1-1;
637 FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
638 if(id!=0){
639 if(ics[lprev+id-1]==ml1){
640 for(l2=0;l2<mt;l2++){
641 l=mt*l1+l2;
642 voldt[l+mt**nk*i]=vold[l];
643 }
644 continue;
645 }
646 }
647 for(l2=0;l2<mt;l2++){
648 l=mt*l1+l2;
649 voldt[l+mt**nk*i]=vold[l];
650 }
651 }
652 }
653 }
654 }
655
656 icntrl=-2;
657
658 FORTRAN(rectcylvold,(co,voldt,cs,&icntrl,mi,iznode,&nznode,nsectors,nk));
659
660 memcpy(&vold[0],&voldt[0],sizeof(double)*mt*nkt);
661 SFREE(voldt);
662
663 }
664
665 /* copying zdof into z */
666
667 if(*nherm==1){
668 RENEW(z,double,(long long)*nev**nzdof);
669 memcpy(&z[0],&zdof[0],(long long)sizeof(double)**nev**nzdof);
670 }else{
671 RENEW(z,double,(long long)2**nev**nzdof);
672 memcpy(&z[0],&zdof[0],(long long)sizeof(double)*2**nev**nzdof);
673 for(i=0;i<*nev;i++){
674 for(j=0;j<*nev;j++){
675 xmr[i**nev+j]/=(sumi[i]*sumi[j]);
676 xmi[i**nev+j]/=(sumi[i]*sumi[j]);
677 }
678 }
679 SFREE(sumi);
680 }
681 SFREE(zdof);
682
683 /* copying the multiple point constraints */
684
685 *nmpc=0;
686 *mpcend=0;
687 for(i=0;i<*nsectors;i++){
688 if(i==0){
689 ileft=*nsectors-1;
690 }else{
691 ileft=i-1;
692 }
693 for(j=0;j<*nmpcold;j++){
694 if(noderight_>10){
695 noderight_=10;
696 RENEW(noderight,ITG,noderight_);
697 RENEW(coefright,double,noderight_);
698 }
699 ipompc[*nmpc]=*mpcend+1;
700 ikmpc[*nmpc]=ikmpc[j]+8*i**nk;
701 ilmpc[*nmpc]=ilmpc[j]+i**nmpcold;
702 strcpy1(&labmpc[20**nmpc],&labmpcold[20*j],20);
703 if(strcmp1(&labmpcold[20*j],"CYCLIC")==0){
704
705 /* identifying the nodes on the master side
706 corresponding to a given slave node */
707
708 index=ipompcold[j]-1;
709 nodeleft=nodempcold[3*index];
710
711 /* slave node term direction */
712
713 idir=nodempcold[3*index+1];
714 index=nodempcold[3*index+2]-1;
715 numnodes=0;
716 do{
717 node=nodempcold[3*index];
718 if(nodempcold[3*index+1]==idir){
719 noderight[numnodes]=node;
720 coefright[numnodes]=coefmpcold[index];
721 numnodes++;
722 if(numnodes>=noderight_){
723 noderight_=(ITG)(1.5*noderight_);
724 RENEW(noderight,ITG,noderight_);
725 RENEW(coefright,double,noderight_);
726 }
727 }else{
728
729 /* master node term direction */
730
731 idirnew=nodempcold[3*index+1];
732 }
733 index=nodempcold[3*index+2]-1;
734 if(index==-1) break;
735 }while(1);
736
737 /* if not successful:
738 The cyclic MPC consists of one term corresponding
739 to the slave node and one or more corresponding to
740 the master node; the direction of the first master
741 node term can correspond to the direction of the
742 slave node term, but does not have to. Any master
743 node term direction will actually do */
744
745 if(numnodes==0){
746
747 /* identifying the nodes on the master side
748 corresponding to a given slave node */
749
750 index=ipompcold[j]-1;
751 nodeleft=nodempcold[3*index];
752 idir=idirnew;
753 index=nodempcold[3*index+2]-1;
754 numnodes=0;
755 do{
756 node=nodempcold[3*index];
757 if(nodempcold[3*index+1]==idir){
758 noderight[numnodes]=node;
759 coefright[numnodes]=coefmpcold[index];
760 numnodes++;
761 if(numnodes>=noderight_){
762 noderight_=(ITG)(1.5*noderight_);
763 RENEW(noderight,ITG,noderight_);
764 RENEW(coefright,double,noderight_);
765 }
766 }
767 index=nodempcold[3*index+2]-1;
768 if(index==-1) break;
769 }while(1);
770
771 }
772
773 if(numnodes>0){
774 sum=0.;
775 for(k=0;k<numnodes;k++){
776 sum+=coefright[k];
777 }
778 for(k=0;k<numnodes;k++){
779 coefright[k]/=sum;
780 }
781 }else{coefright[0]=1.;}
782 nodempc[3**mpcend]=nodeleft+i**nk;
783 nodempc[3**mpcend+1]=idir;
784 nodempc[3**mpcend+2]=*mpcend+2;
785 coefmpc[*mpcend]=1.;
786 for(k=0;k<numnodes;k++){
787 (*mpcend)++;
788 nodempc[3**mpcend]=noderight[k]+ileft**nk;
789 nodempc[3**mpcend+1]=idir;
790 nodempc[3**mpcend+2]=*mpcend+2;
791 coefmpc[*mpcend]=-coefright[k];
792 }
793 nodempc[3**mpcend+2]=0;
794 (*mpcend)++;
795 }else{
796 index=ipompcold[j]-1;
797 iterm=0;
798 do{
799 iterm++;
800 node=nodempcold[3*index];
801 idir=nodempcold[3*index+1];
802 coef=coefmpcold[index];
803
804 /* check whether SUBCYCLIC MPC: if the current node
805 is an independent node of a CYCLIC MPC, the
806 node in the new MPC should be the cylic previous
807 one */
808
809 nodenew=node+i**nk;
810 if(strcmp1(&labmpcold[20*j],"SUBCYCLIC")==0){
811 for(ij=0;ij<*mcs;ij++){
812 lprev=cs[ij*17+13];
813 ilength=cs[ij*17+3];
814 FORTRAN(nident,(&ics[lprev],&node,&ilength,&id));
815 if(id!=0){
816 if(ics[lprev+id-1]==node){
817 nodenew=node+ileft**nk;
818 break;
819 }
820 }
821 }
822 }
823
824 /* modification of the MPC coefficients if
825 cylindrical coordinate system is active
826 it is assumed that all terms in the MPC are
827 either in the radial, the circumferential
828 or axial direction */
829
830 if(*ntrans<=0){itr=0;}
831 else if(inotr[2*node-2]==0){itr=0;}
832 else{itr=inotr[2*node-2];}
833
834 if(iterm==1) locdir=-1;
835
836 if((itr!=0)&&(idir!=0)){
837 if(trab[7*itr-1]<0){
838 FORTRAN(transformatrix,(&trab[7*itr-7],
839 &co[3*node-3],a));
840 if(iterm==1){
841 for(k=0;k<3;k++){
842 if(fabs(a[3*k+idir-1]-coef)<1.e-10){
843 FORTRAN(transformatrix,(&trab[7*itr-7],
844 &co[3*nodenew-3],a));
845 coef=a[3*k+idir-1];
846 locdir=k;
847 break;
848 }
849 if(fabs(a[3*k+idir-1]+coef)<1.e-10){
850 FORTRAN(transformatrix,(&trab[7*itr-7],
851 &co[3*nodenew-3],a));
852 coef=-a[3*k+idir-1];
853 locdir=k;
854 break;
855 }
856 }
857 }else{
858 if(locdir!=-1){
859 if(fabs(a[3*locdir+idir-1])>1.e-10){
860 ratio=coef/a[3*locdir+idir-1];
861 }else{ratio=0.;}
862 FORTRAN(transformatrix,(&trab[7*itr-7],
863 &co[3*nodenew-3],a));
864 coef=ratio*a[3*locdir+idir-1];
865 }
866 }
867 }
868 }
869
870 nodempc[3**mpcend]=nodenew;
871 nodempc[3**mpcend+1]=idir;
872 coefmpc[*mpcend]=coef;
873 index=nodempcold[3*index+2]-1;
874 if(index==-1) break;
875 nodempc[3**mpcend+2]=*mpcend+2;
876 (*mpcend)++;
877 }while(1);
878 nodempc[3**mpcend+2]=0;
879 (*mpcend)++;
880 }
881 (*nmpc)++;
882 }
883 }
884
885 /* copying the temperatures */
886
887 if(*ithermal!=0){
888 for(i=1;i<*nsectors;i++){
889 lint=i**nk;
890 for(j=0;j<*nk;j++){
891 t0[lint+j]=t0[j];
892 t1old[lint+j]=t1old[j];
893 t1[lint+j]=t1[j];
894 }
895 }
896 if(*nam>0){
897 for(i=1;i<*nsectors;i++){
898 lint=i**nk;
899 for(j=0;j<*nk;j++){
900 iamt1[lint+j]=iamt1[j];
901 }
902 }
903 }
904 }
905
906 /* copying the contact definition */
907
908 // if(*nmethod==4){
909
910 /* first find the startposition to append the expanded contact fields*/
911
912 /* for(j=0; j<*nset; j++){
913 if(iendset[j]>tint){
914 tint=iendset[j];
915 }
916 }
917 tint++;*/
918 /* now append and expand the contact definitons*/
919 /* NNEW(tchar1,char,81);
920 NNEW(tchar2,char,81);
921 NNEW(tchar3,char,81);
922 for(i=0; i<*ntie; i++){
923 if(tieset[i*(81*3)+80]=='C'){
924 memcpy(tchar2,&tieset[i*(81*3)+81],81);
925 tchar2[80]='\0';
926 memcpy(tchar3,&tieset[i*(81*3)+81+81],81);
927 tchar3[80]='\0';*/
928 //a contact constraint was found, so append and expand the information
929 /* for(j=0; j<*nset; j++){
930 memcpy(tchar1,&set[j*81],81);
931 tchar1[80]='\0';
932 if(strcmp(tchar1,tchar2)==0){*/
933 /* dependent nodal surface was found,copy the original information first */
934 /* tnstart=tint;
935 for(k=0; k<iendset[j]-istartset[j]+1; k++){
936 ialset[tint-1]=ialset[istartset[j]-1+k];
937 tint++;
938 }*/
939 /* now append the expanded information */
940 /* for(l=1; l<*nsectors; l++){
941 for(k=0; k<iendset[j]-istartset[j]+1; k++){
942 ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?ialset[istartset[j]-1+k]+*nk*l:-1;
943 tint++;
944 }
945 }
946 tnend=tint-1;*/
947 /* now replace the information in istartset and iendset*/
948 /* istartset[j]=tnstart;
949 iendset[j]=tnend;
950 }
951 else if(strcmp(tchar1,tchar3)==0){*/
952 /* independent element face surface was found */
953 /* tnstart=tint;
954 for(k=0; k<iendset[j]-istartset[j]+1; k++){
955 ialset[tint-1]=ialset[istartset[j]-1+k];
956 tint++;
957 }*/
958 /* now append the expanded information*/
959 /* for(l=1; l<*nsectors; l++){
960 for(k=0; k<iendset[j]-istartset[j]+1; k++){
961 tint2=((ITG)(ialset[istartset[j]-1+k]))/10;
962 ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?(tint2+*ne*l)*10+(ialset[istartset[j]-1+k]-(tint2*10)):-1;
963 tint++;
964 }
965 }
966 tnend=tint-1;*/
967 /* now replace the information in istartset and iendset*/
968 /* istartset[j]=tnstart;
969 iendset[j]=tnend;
970 }
971 }
972 }
973 }
974 SFREE(tchar1);
975 SFREE(tchar2);
976 SFREE(tchar3);
977 } */
978
979 *nk=nkt;
980 (*ne)*=(*nsectors);
981 (*nkon)*=(*nsectors);
982 (*nboun)*=(*nsectors);
983 neq[1]=neqh**nsectors;
984
985 *zp=z;*izdofp=izdof;
986
987 SFREE(temp_array);SFREE(coefmpcnew);SFREE(noderight);SFREE(coefright);
988 SFREE(v);SFREE(vt);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
989 SFREE(inocs);SFREE(ielcs);SFREE(filabt);SFREE(iznode);SFREE(csmass);
990
991 return;
992 }
993