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
frdcyc(double * co,ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,double * v,double * stn,ITG * inum,ITG * nmethod,ITG * kode,char * filab,double * een,double * t1,double * fn,double * time,double * epn,ITG * ielmat,char * matname,double * cs,ITG * mcs,ITG * nkon,double * enern,double * xstaten,ITG * nstate_,ITG * istep,ITG * iinc,ITG * iperturb,double * ener,ITG * mi,char * output,ITG * ithermal,double * qfn,ITG * ialset,ITG * istartset,ITG * iendset,double * trab,ITG * inotr,ITG * ntrans,double * orab,ITG * ielorien,ITG * norien,double * sti,double * veold,ITG * noddiam,char * set,ITG * nset,double * emn,double * thicke,char * jobnamec,ITG * ne0,double * cdn,ITG * mortar,ITG * nmat,double * qfx,ITG * ielprop,double * prop)24 void frdcyc(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne,double *v,
25 double *stn,ITG *inum,ITG *nmethod,ITG *kode,char *filab,
26 double *een,double *t1,double *fn,double *time,double *epn,
27 ITG *ielmat,char *matname, double *cs, ITG *mcs, ITG *nkon,
28 double *enern, double *xstaten, ITG *nstate_, ITG *istep,
29 ITG *iinc, ITG *iperturb, double *ener, ITG *mi, char *output,
30 ITG *ithermal, double *qfn, ITG *ialset, ITG *istartset,
31 ITG *iendset, double *trab, ITG *inotr, ITG *ntrans,
32 double *orab, ITG *ielorien, ITG *norien, double *sti,
33 double *veold, ITG *noddiam,char *set,ITG *nset, double *emn,
34 double *thicke,char* jobnamec,ITG *ne0,double *cdn,ITG *mortar,
35 ITG *nmat,double *qfx,ITG *ielprop,double *prop){
36
37 /* duplicates fields for static cyclic symmetric calculations */
38
39 char *lakont=NULL,description[13]=" ";
40
41 ITG nkt,icntrl,*kont=NULL,*ipkont=NULL,*inumt=NULL,*ielmatt=NULL,net,i,l,
42 imag=0,mode=-1,ngraph,*inocs=NULL,*ielcs=NULL,l1,l2,is,
43 jj,node,i1,i2,nope,iel,indexe,j,ielset,*inotrt=NULL,mt=mi[1]+1,
44 *ipneigh=NULL,*neigh=NULL,net0;
45
46 double *vt=NULL,*fnt=NULL,*stnt=NULL,*eent=NULL,*cot=NULL,*t1t=NULL,
47 *epnt=NULL,*enernt=NULL,*xstatent=NULL,theta,pi,t[3],*qfnt=NULL,
48 *vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,*vmax=NULL,*stnmax=NULL,
49 *stit=NULL,*eenmax=NULL,*fnr=NULL,*fni=NULL,*emnt=NULL,
50 *cdnr=NULL,*cdni=NULL;
51
52 pi=4.*atan(1.);
53
54 /* determining the maximum number of sectors to be plotted */
55
56 ngraph=1;
57 for(j=0;j<*mcs;j++){
58 if(cs[17*j+4]>ngraph) ngraph=cs[17*j+4];
59 }
60
61 /* assigning nodes and elements to sectors */
62
63 NNEW(inocs,ITG,*nk);
64 NNEW(ielcs,ITG,*ne);
65 ielset=cs[12];
66 if((*mcs!=1)||(ielset!=0)){
67 for(i=0;i<*nk;i++) inocs[i]=-1;
68 for(i=0;i<*ne;i++) ielcs[i]=-1;
69 }
70
71 for(i=0;i<*mcs;i++){
72 is=cs[17*i+4];
73 if(is==1) continue;
74 ielset=cs[17*i+12];
75 if(ielset==0) continue;
76 for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
77 if(ialset[i1]>0){
78 iel=ialset[i1]-1;
79 if(ipkon[iel]<0) continue;
80 ielcs[iel]=i;
81 indexe=ipkon[iel];
82 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
83 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
84 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
85 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
86 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
87 else {nope=6;}
88 for(i2=0;i2<nope;++i2){
89 node=kon[indexe+i2]-1;
90 inocs[node]=i;
91 }
92 }
93 else{
94 iel=ialset[i1-2]-1;
95 do{
96 iel=iel-ialset[i1];
97 if(iel>=ialset[i1-1]-1) break;
98 if(ipkon[iel]<0) continue;
99 ielcs[iel]=i;
100 indexe=ipkon[iel];
101 if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
102 else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
103 else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
104 else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
105 else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
106 else {nope=6;}
107 for(i2=0;i2<nope;++i2){
108 node=kon[indexe+i2]-1;
109 inocs[node]=i;
110 }
111 }while(1);
112 }
113 }
114 }
115
116 NNEW(cot,double,3**nk*ngraph);
117 if(*ntrans>0)NNEW(inotrt,ITG,2**nk*ngraph);
118
119 if((strcmp1(&filab[0],"U ")==0)||
120 (strcmp1(&filab[1131],"TT ")==0)||
121 (strcmp1(&filab[1218],"MF ")==0)||
122 (strcmp1(&filab[1305],"PT ")==0)||
123 ((strcmp1(&filab[87],"NT ")==0)&&(*ithermal>=2)))
124 NNEW(vt,double,mt**nk*ngraph);
125 if((strcmp1(&filab[87],"NT ")==0)&&(*ithermal<2))
126 NNEW(t1t,double,*nk*ngraph);
127 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
128 (strcmp1(&filab[1044],"ERR ")==0))
129 NNEW(stnt,double,6**nk*ngraph);
130 if(strcmp1(&filab[261],"E ")==0)
131 NNEW(eent,double,6**nk*ngraph);
132 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[783],"RFL ")==0))
133 NNEW(fnt,double,mt**nk*ngraph);
134 if(strcmp1(&filab[435],"PEEQ")==0)
135 NNEW(epnt,double,*nk*ngraph);
136 if(strcmp1(&filab[522],"ENER")==0)
137 NNEW(enernt,double,*nk*ngraph);
138 if(strcmp1(&filab[609],"SDV ")==0)
139 NNEW(xstatent,double,*nstate_**nk*ngraph);
140 if((strcmp1(&filab[696],"HFL ")==0)||(strcmp1(&filab[2784],"HER ")==0))
141 NNEW(qfnt,double,3**nk*ngraph);
142 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
143 (strcmp1(&filab[2175],"CONT")==0))
144 NNEW(stit,double,6*mi[0]**ne*ngraph);
145 if(strcmp1(&filab[2697],"ME ")==0)
146 NNEW(emnt,double,6**nk*ngraph);
147
148 /* the topology only needs duplication the first time it is
149 stored in the frd file (*kode=1)
150 the above two lines are not true: lakon is needed for
151 contact information in frd.f */
152
153 NNEW(kont,ITG,*nkon*ngraph);
154 NNEW(ipkont,ITG,*ne*ngraph);
155 NNEW(lakont,char,8**ne*ngraph);
156 NNEW(ielmatt,ITG,mi[2]**ne*ngraph);
157 NNEW(inumt,ITG,*nk*ngraph);
158
159 nkt=ngraph**nk;
160 net0=(ngraph-1)**ne+(*ne0);
161 net=ngraph**ne;
162
163 /* copying the coordinates of the first sector */
164
165 for(l=0;l<3**nk;l++){cot[l]=co[l];}
166 if(*ntrans>0){for(l=0;l<*nk;l++){inotrt[2*l]=inotr[2*l];}}
167
168 /* copying the topology of the first sector */
169
170 for(l=0;l<*nkon;l++){kont[l]=kon[l];}
171 for(l=0;l<*ne;l++){ipkont[l]=ipkon[l];}
172 for(l=0;l<8**ne;l++){lakont[l]=lakon[l];}
173 for(l=0;l<mi[2]**ne;l++){ielmatt[l]=ielmat[l];}
174
175 /* generating the coordinates for the other sectors */
176
177 icntrl=1;
178
179 FORTRAN(rectcyl,(cot,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
180
181 for(jj=0;jj<*mcs;jj++){
182 is=cs[17*jj+4];
183 for(i=1;i<is;i++){
184
185 theta=i*2.*pi/cs[17*jj];
186
187 for(l=0;l<*nk;l++){
188 if(inocs[l]==jj){
189 cot[3*l+i*3**nk]=cot[3*l];
190 cot[1+3*l+i*3**nk]=cot[1+3*l]+theta;
191 cot[2+3*l+i*3**nk]=cot[2+3*l];
192 }
193 }
194
195 if(*ntrans>0){
196 for(l=0;l<*nk;l++){
197 if(inocs[l]==jj){
198 inotrt[2*l+i*2**nk]=inotrt[2*l];
199 }
200 }
201 }
202
203 for(l=0;l<*nkon;l++){kont[l+i**nkon]=kon[l]+i**nk;}
204 for(l=0;l<*ne;l++){
205 if(ielcs[l]==jj){
206 if(ipkon[l]>=0){
207 ipkont[l+i**ne]=ipkon[l]+i**nkon;
208 ielmatt[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
209 for(l1=0;l1<8;l1++){
210 l2=8*l+l1;
211 lakont[l2+i*8**ne]=lakon[l2];
212 }
213 }
214 else ipkont[l+i**ne]=-1;
215 }
216 }
217 }
218 }
219
220 icntrl=-1;
221
222 FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
223 &imag,mi,emnt));
224
225 /* mapping the results to the other sectors */
226
227 for(l=0;l<*nk;l++){inumt[l]=inum[l];}
228
229 icntrl=2;
230
231 FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filab,&imag,mi,emn));
232
233 if((strcmp1(&filab[0],"U ")==0)||
234 (strcmp1(&filab[1131],"TT ")==0)||
235 (strcmp1(&filab[1218],"MF ")==0)||
236 (strcmp1(&filab[1305],"PT ")==0)||
237 ((strcmp1(&filab[87],"NT ")==0)&&(*ithermal>=2)))
238 for(l=0;l<mt**nk;l++){vt[l]=v[l];};
239 if((strcmp1(&filab[87],"NT ")==0)&&(*ithermal<2))
240 for(l=0;l<*nk;l++){t1t[l]=t1[l];};
241 if(strcmp1(&filab[174],"S ")==0)
242 for(l=0;l<6**nk;l++){stnt[l]=stn[l];};
243 if(strcmp1(&filab[261],"E ")==0)
244 for(l=0;l<6**nk;l++){eent[l]=een[l];};
245 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[783],"RFL ")==0))
246 for(l=0;l<mt**nk;l++){fnt[l]=fn[l];};
247 if(strcmp1(&filab[435],"PEEQ")==0)
248 for(l=0;l<*nk;l++){epnt[l]=epn[l];};
249 if(strcmp1(&filab[522],"ENER")==0)
250 for(l=0;l<*nk;l++){enernt[l]=enern[l];};
251 if(strcmp1(&filab[609],"SDV ")==0)
252 for(l=0;l<*nstate_**nk;l++){xstatent[l]=xstaten[l];};
253 if((strcmp1(&filab[696],"HFL ")==0)||(strcmp1(&filab[2784],"HER ")==0))
254 for(l=0;l<3**nk;l++){qfnt[l]=qfn[l];};
255 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
256 (strcmp1(&filab[2175],"CONT")==0))
257 for(l=0;l<6*mi[0]**ne;l++){stit[l]=sti[l];};
258 if(strcmp1(&filab[2697],"ME ")==0)
259 for(l=0;l<6**nk;l++){emnt[l]=emn[l];};
260
261 for(jj=0;jj<*mcs;jj++){
262 is=cs[17*jj+4];
263 for(i=1;i<is;i++){
264
265 for(l=0;l<*nk;l++){inumt[l+i**nk]=inum[l];}
266
267 if((strcmp1(&filab[0],"U ")==0)||
268 (strcmp1(&filab[1131],"TT ")==0)||
269 (strcmp1(&filab[1218],"MF ")==0)||
270 (strcmp1(&filab[1305],"PT ")==0)||
271 ((strcmp1(&filab[87],"NT ")==0)&&(*ithermal>=2))){
272 for(l1=0;l1<*nk;l1++){
273 if(inocs[l1]==jj){
274 for(l2=0;l2<4;l2++){
275 l=mt*l1+l2;
276 vt[l+mt**nk*i]=v[l];
277 }
278 }
279 }
280 }
281
282 if((strcmp1(&filab[87],"NT ")==0)&&(*ithermal<2)){
283 for(l=0;l<*nk;l++){
284 if(inocs[l]==jj) t1t[l+*nk*i]=t1[l];
285 }
286 }
287
288 if(strcmp1(&filab[174],"S ")==0){
289 for(l1=0;l1<*nk;l1++){
290 if(inocs[l1]==jj){
291 for(l2=0;l2<6;l2++){
292 l=6*l1+l2;
293 stnt[l+6**nk*i]=stn[l];
294 }
295 }
296 }
297 }
298
299 if(strcmp1(&filab[261],"E ")==0){
300 for(l1=0;l1<*nk;l1++){
301 if(inocs[l1]==jj){
302 for(l2=0;l2<6;l2++){
303 l=6*l1+l2;
304 eent[l+6**nk*i]=een[l];
305 }
306 }
307 }
308 }
309
310 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[783],"RFL ")==0)){
311 for(l1=0;l1<*nk;l1++){
312 if(inocs[l1]==jj){
313 for(l2=0;l2<4;l2++){
314 l=mt*l1+l2;
315 fnt[l+mt**nk*i]=fn[l];
316 }
317 }
318 }
319 }
320
321 if(strcmp1(&filab[435],"PEEQ")==0){
322 for(l=0;l<*nk;l++){
323 if(inocs[l]==jj) epnt[l+*nk*i]=epn[l];
324 }
325 }
326
327 if(strcmp1(&filab[522],"ENER")==0){
328 for(l=0;l<*nk;l++){
329 if(inocs[l]==jj) enernt[l+*nk*i]=enern[l];
330 }
331 }
332
333 if(strcmp1(&filab[609],"SDV ")==0){
334 for(l1=0;l1<*nk;l1++){
335 if(inocs[l1]==jj){
336 for(l2=0;l2<*nstate_;l2++){
337 l=*nstate_*l1+l2;
338 xstatent[l+*nstate_**nk*i]=xstaten[l];
339 }
340 }
341 }
342 }
343
344 if((strcmp1(&filab[696],"HFL ")==0)||(strcmp1(&filab[2784],"HER ")==0)){
345 for(l1=0;l1<*nk;l1++){
346 if(inocs[l1]==jj){
347 for(l2=0;l2<3;l2++){
348 l=3*l1+l2;
349 qfnt[l+3**nk*i]=qfn[l];
350 }
351 }
352 }
353 }
354
355 if(strcmp1(&filab[2697],"ME ")==0){
356 for(l1=0;l1<*nk;l1++){
357 if(inocs[l1]==jj){
358 for(l2=0;l2<6;l2++){
359 l=6*l1+l2;
360 emnt[l+6**nk*i]=emn[l];
361 }
362 }
363 }
364 }
365 }
366 }
367
368 icntrl=-2;
369
370 FORTRAN(rectcyl,(cot,vt,fnt,stnt,qfnt,eent,cs,&nkt,&icntrl,t,filab,
371 &imag,mi,emnt));
372
373 if(strcmp1(&filab[1044],"ZZS")==0){
374 NNEW(neigh,ITG,40*net);
375 NNEW(ipneigh,ITG,nkt);
376 }
377
378 frd(cot,&nkt,kont,ipkont,lakont,&net0,vt,stnt,inumt,nmethod,
379 kode,filab,eent,t1t,fnt,time,epnt,ielmatt,matname,enernt,xstatent,
380 nstate_,istep,iinc,ithermal,qfnt,&mode,noddiam,trab,inotrt,
381 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
382 mi,stit,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,&net,
383 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emnt,
384 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,ielprop,
385 prop,sti);
386
387 if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
388
389 if((strcmp1(&filab[0],"U ")==0)||
390 (strcmp1(&filab[1131],"TT ")==0)||
391 (strcmp1(&filab[1218],"MF ")==0)||
392 (strcmp1(&filab[1305],"PT ")==0)||
393 ((strcmp1(&filab[87],"NT ")==0)&&(*ithermal>=2))) SFREE(vt);
394 if((strcmp1(&filab[87],"NT ")==0)&&(*ithermal<2)) SFREE(t1t);
395 if((strcmp1(&filab[174],"S ")==0)||(strcmp1(&filab[1044],"ZZS ")==0)||
396 (strcmp1(&filab[1044],"ERR ")==0))
397 SFREE(stnt);
398 if(strcmp1(&filab[261],"E ")==0) SFREE(eent);
399 if((strcmp1(&filab[348],"RF ")==0)||(strcmp1(&filab[783],"RFL ")==0))
400 SFREE(fnt);
401 if(strcmp1(&filab[435],"PEEQ")==0) SFREE(epnt);
402 if(strcmp1(&filab[522],"ENER")==0) SFREE(enernt);
403 if(strcmp1(&filab[609],"SDV ")==0) SFREE(xstatent);
404 if((strcmp1(&filab[696],"HFL ")==0)||(strcmp1(&filab[2784],"HER ")==0)) SFREE(qfnt);
405 if((strcmp1(&filab[1044],"ZZS ")==0)||(strcmp1(&filab[1044],"ERR ")==0)||
406 (strcmp1(&filab[2175],"CONT")==0)) SFREE(stit);
407 if(strcmp1(&filab[2697],"ME ")==0) SFREE(emnt);
408
409 SFREE(kont);SFREE(ipkont);SFREE(lakont);SFREE(ielmatt);
410 SFREE(inumt);SFREE(cot);if(*ntrans>0)SFREE(inotrt);
411 SFREE(inocs);SFREE(ielcs);
412 return;
413 }
414
415