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 <stdlib.h>
19 #include <math.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include "CalculiX.h"
23
24 #define min(a,b) ((a) <= (b) ? (a) : (b))
25 #define max(a,b) ((a) >= (b) ? (a) : (b))
26
frd(double * co,ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne0,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 * enern,double * xstaten,ITG * nstate_,ITG * istep,ITG * iinc,ITG * ithermal,double * qfn,ITG * mode,ITG * noddiam,double * trab,ITG * inotr,ITG * ntrans,double * orab,ITG * ielorien,ITG * norien,char * description,ITG * ipneigh,ITG * neigh,ITG * mi,double * stx,double * vr,double * vi,double * stnr,double * stni,double * vmax,double * stnmax,ITG * ngraph,double * veold,double * ener,ITG * ne,double * cs,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,double * eenmax,double * fnr,double * fni,double * emn,double * thicke,char * jobnamec,char * output,double * qfx,double * cdn,ITG * mortar,double * cdnr,double * cdni,ITG * nmat,ITG * ielprop,double * prop,double * sti)27 void frd(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne0,
28 double *v,double *stn,ITG *inum,ITG *nmethod,ITG *kode,
29 char *filab,double *een,double *t1,double *fn,double *time,
30 double *epn,ITG *ielmat,char *matname,double *enern,
31 double *xstaten,ITG *nstate_,ITG *istep,ITG *iinc,
32 ITG *ithermal,double *qfn,ITG *mode,ITG *noddiam,
33 double *trab,ITG *inotr,ITG *ntrans,double *orab,
34 ITG *ielorien,ITG *norien,char *description,ITG *ipneigh,
35 ITG *neigh,ITG *mi,double *stx,double *vr,double *vi,
36 double *stnr,double *stni,double *vmax,double *stnmax,
37 ITG *ngraph,double *veold,double *ener,ITG *ne,double *cs,
38 char *set,ITG *nset,ITG *istartset,ITG *iendset,ITG *ialset,
39 double *eenmax,double *fnr,double *fni,double *emn,
40 double *thicke,char *jobnamec,char *output,double *qfx,
41 double *cdn,ITG *mortar,double *cdnr,double *cdni,ITG *nmat,
42 ITG *ielprop,double *prop,double *sti){
43
44 /* stores the results in frd format
45
46 iselect selects which nodes are to be stored:
47 iselect=-1 means only those nodes for which inum negative
48 ist, i.e. network nodes
49 iselect=+1 means only those nodes for which inum positive
50 ist, i.e. structural nodes
51 iselect=0 means both of the above */
52
53 FILE *f1;
54
55 char c[2]="C",m1[4]=" -1",m2[4]=" -2",m3[4]=" -3",
56 p0[6]=" 0",p1[6]=" 1",p2[6]=" 2",p3[6]=" 3",p4[6]=" 4",
57 p5[6]=" 5",p6[6]=" 6",p7[6]=" 7",p8[6]=" 8",p9[6]=" 9",
58 p10[6]=" 10",p11[6]=" 11",
59 p12[6]=" 12", fneig[132]="",date[8],clock[10],newdate[21],newclock[9],
60 material[59]=" ",
61 text[2]=" ";
62
63 static ITG icounter=0,nkcoords,iaxial;
64
65 ITG null,one,i,j,k,indexe,nemax,nlayer,noutloc,iset,iselect,ncomp,nope,
66 nodes,ifield[7],nfield[2],icomp[7],ifieldstate[*nstate_],two,three,
67 icompstate[*nstate_],ip0=0,ip1=1,ip2=2,ip3=3,ip4=4,ip5=5,ip6=6,ip7=7,
68 ip8=8,ip9=9,ip10=10,ip11=11,ip12=12,imat,nelout,ioutall=0,*inumshell=NULL,
69 nterms,nout,noutplus,noutmin,mt=mi[1]+1,iflag,numfield,iorienloc;
70
71 ITG ncompscalar=1,ifieldscalar[1]={1},icompscalar[1]={0},
72 nfieldscalar[2]={1,0};
73 ITG ncompvector=3,ifieldvector[3]={1,1,1},icompvector[3]={0,1,2},
74 nfieldvector1[2]={3,0},nfieldvector0[2]={mi[1]+1,0},
75 icompvectorlast[3]={3,4,5};
76 ITG ncomptensor=6,ifieldtensor[6]={1,1,1,1,1,1},icomptensor[6]={0,1,2,3,5,4},
77 nfieldtensor[2]={6,0};
78 ITG ncompscalph=2,ifieldscalph[2]={1,2},icompscalph[2]={0,0},
79 nfieldscalph[2]={0,0};
80 ITG ncompvectph=6,ifieldvectph[6]={1,1,1,2,2,2},icompvectph[6]={1,2,3,1,2,3},
81 nfieldvectph[2]={mi[1]+1,mi[1]+1};
82 ITG ncomptensph=12,ifieldtensph[12]={1,1,1,1,1,1,2,2,2,2,2,2},
83 icomptensph[12]={0,1,2,3,5,4,0,1,2,3,5,4},nfieldtensph[2]={6,6};
84
85 int iw;
86
87 float fl;
88
89 double pi,oner,*errn=NULL,*ethn=NULL;
90
91 strcpy(fneig,jobnamec);
92 strcat(fneig,".frd");
93
94 if((f1=fopen(fneig,"ab"))==NULL){
95 printf("*EOR in frd: cannot open frd file for writing...");
96 exit(0);
97 }
98
99 /* check whether all results have to be stored (also those
100 corresponding to inactive nodes or elements) */
101
102 if(strcmp1(&output[3],"a")==0) ioutall=1;
103
104 pi=4.*atan(1.);
105 null=0;
106 one=1;two=2;three=3;
107 oner=1.;
108
109 /* determining nout, noutplus and noutmin
110 nout: number of structural and network nodes
111 noutplus: number of structural nodes
112 noutmin: number of network nodes */
113
114 if(*nmethod!=0){
115 nout=0;
116 noutplus=0;
117 noutmin=0;
118 if(ioutall==0){
119 for(i=0;i<*nk;i++){
120 if(inum[i]==0) continue;
121 nout++;
122 if(inum[i]>0) noutplus++;
123 if(inum[i]<0) noutmin++;
124 }
125 }else{
126 for(i=0;i<*nk;i++){
127 nout++;
128 if(inum[i]>0) noutplus++;
129 if(inum[i]<0) noutmin++;
130 }
131 }
132 }else{
133 nout=*nk;
134 }
135
136 /* first time something is written in the frd-file: store
137 computational metadata, the nodal coordinates and the
138 topology */
139
140 if((*kode==1)&&((*nmethod!=5)||(*mode!=0))){
141 iaxial=0.;
142
143 /* date and time */
144
145 FORTRAN(dattime,(date,clock));
146
147 for(i=0;i<20;i++) newdate[i]=' ';
148 newdate[20]='\0';
149
150 strcpy1(newdate,&date[6],2);
151 strcpy1(&newdate[2],".",1);
152 if(strcmp1(&date[4],"01")==0){
153 strcpy1(&newdate[3],"january.",8);
154 strcpy1(&newdate[11],&date[0],4);
155 }else if(strcmp1(&date[4],"02")==0){
156 strcpy1(&newdate[3],"february.",9);
157 strcpy1(&newdate[12],&date[0],4);
158 }else if(strcmp1(&date[4],"03")==0){
159 strcpy1(&newdate[3],"march.",6);
160 strcpy1(&newdate[9],&date[0],4);
161 }else if(strcmp1(&date[4],"04")==0){
162 strcpy1(&newdate[3],"april.",6);
163 strcpy1(&newdate[9],&date[0],4);
164 }else if(strcmp1(&date[4],"05")==0){
165 strcpy1(&newdate[3],"may.",4);
166 strcpy1(&newdate[7],&date[0],4);
167 }else if(strcmp1(&date[4],"06")==0){
168 strcpy1(&newdate[3],"june.",5);
169 strcpy1(&newdate[8],&date[0],4);
170 }else if(strcmp1(&date[4],"07")==0){
171 strcpy1(&newdate[3],"july.",5);
172 strcpy1(&newdate[8],&date[0],4);
173 }else if(strcmp1(&date[4],"08")==0){
174 strcpy1(&newdate[3],"august.",7);
175 strcpy1(&newdate[10],&date[0],4);
176 }else if(strcmp1(&date[4],"09")==0){
177 strcpy1(&newdate[3],"september.",10);
178 strcpy1(&newdate[13],&date[0],4);
179 }else if(strcmp1(&date[4],"10")==0){
180 strcpy1(&newdate[3],"october.",8);
181 strcpy1(&newdate[11],&date[0],4);
182 }else if(strcmp1(&date[4],"11")==0){
183 strcpy1(&newdate[3],"november.",9);
184 strcpy1(&newdate[12],&date[0],4);
185 }else if(strcmp1(&date[4],"12")==0){
186 strcpy1(&newdate[3],"december.",9);
187 strcpy1(&newdate[12],&date[0],4);
188 }
189
190 strcpy1(newclock,clock,2);
191 strcpy1(&newclock[2],":",1);
192 strcpy1(&newclock[3],&clock[2],2);
193 strcpy1(&newclock[5],":",1);
194 strcpy1(&newclock[6],&clock[4],2);
195 newclock[8]='\0';
196
197 fprintf(f1,"%5sUUSER \n",p1);
198 fprintf(f1,"%5sUDATE %20s \n",p1,newdate);
199 fprintf(f1,"%5sUTIME %8s \n",p1,newclock);
200 fprintf(f1,"%5sUHOST \n",p1);
201 fprintf(f1,"%5sUPGM CalculiX \n",p1);
202 fprintf(f1,"%5sUVERSION Version 2.18 \n",p1);
203 fprintf(f1,"%5sUCOMPILETIME Wed Sep 15 21:41:41 CEST 2021 \n",p1);
204 fprintf(f1,"%5sUDIR \n",p1);
205 fprintf(f1,"%5sUDBN \n",p1);
206
207 for(i=0;i<*nmat;i++){
208 strcpy1(material,&matname[80*i],58);
209 fprintf(f1,"%5sUMAT%5" ITGFORMAT "%58s\n",p1,i+1,material);
210 }
211
212 /* storing the header of the coordinates */
213
214 if(strcmp1(output,"asc")==0){
215 fprintf(f1,"%5s%1s %12" ITGFORMAT "%38" ITGFORMAT "\n",p2,c,nout,one);
216 }else{
217 fprintf(f1,"%5s%1s %12" ITGFORMAT "%38" ITGFORMAT "\n",p2,c,nout,three);
218 }
219
220 /* storing the coordinates themselves */
221
222 if(*nmethod!=0){
223 for(i=0;i<*nk;i++){
224 if((inum[i]==0)&&(ioutall!=1)) continue;
225 if(strcmp1(output,"asc")==0){
226 fprintf(f1,"%3s%10" ITGFORMAT "%12.5E%12.5E%12.5E\n",m1,i+1,(float)co[3*i],
227 (float)co[3*i+1],(float)co[3*i+2]);
228 }else{
229 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
230 fwrite(&co[3*i],sizeof(double),1,f1);
231 fwrite(&co[3*i+1],sizeof(double),1,f1);
232 fwrite(&co[3*i+2],sizeof(double),1,f1);
233 }
234 }
235 }else{
236 for(i=0;i<*nk;i++){
237 if(strcmp1(output,"asc")==0){
238 fprintf(f1,"%3s%10" ITGFORMAT "%12.5E%12.5E%12.5E\n",m1,i+1,(float)co[3*i],
239 (float)co[3*i+1],(float)co[3*i+2]);
240 }else{
241 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
242 fwrite(&co[3*i],sizeof(double),1,f1);
243 fwrite(&co[3*i+1],sizeof(double),1,f1);
244 fwrite(&co[3*i+2],sizeof(double),1,f1);
245 }
246 }
247 }
248
249 /* nkcoords is the number of nodes at the time when
250 the nodal coordinates are stored in the frd file */
251
252 nkcoords=*nk;
253 if(strcmp1(output,"asc")==0)fprintf(f1,"%3s\n",m3);
254
255 /* determining the number of elements */
256
257 nelout=0;
258 for(i=0;i<*ne0;i++){
259 if(((ipkon[i]<=-1)&&(ioutall==0))||(ipkon[i]==-1)){
260 continue;
261
262 /* the following elements are not stored in the .frd file: */
263
264 /* contact spring element */
265
266 }else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
267 continue;
268
269 /* film advection element */
270
271 }else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
272 continue;
273
274 /* one-noded spring element */
275
276 }else if((strcmp1(&lakon[8*i],"E")==0)&&
277 (strcmp1(&lakon[8*i+6],"1")==0)){
278 continue;
279
280 /* coupling element */
281
282 }else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
283 continue;
284
285 /* mass element */
286
287 }else if(strcmp1(&lakon[8*i],"MASS")==0){
288 continue;
289
290 /* user element */
291
292 }else if(strcmp1(&lakon[8*i],"U1")==0){
293 continue;
294 }
295 nelout++;
296 }
297
298 /* storing the topology */
299
300 if(strcmp1(output,"asc")==0){
301 fprintf(f1,"%5s%1s %12" ITGFORMAT "%38" ITGFORMAT "\n",p3,c,nelout,one);
302 }else{
303 fprintf(f1,"%5s%1s %12" ITGFORMAT "%38" ITGFORMAT "\n",p3,c,nelout,two);
304 }
305 nemax=*ne0;
306
307 for(i=0;i<*ne0;i++){
308 if(ipkon[i]<=-1){
309
310 if(ipkon[i]==-1){
311 continue;
312 }else if(ioutall!=0){
313 indexe=-2-ipkon[i];
314 }else{
315 continue;
316 }
317
318 /* next element types are not stored */
319
320 }else if(strcmp1(&lakon[8*i],"F")==0){
321 continue;
322 }else if(strcmp1(&lakon[8*i],"ESPRNGC")==0){
323 continue;
324 }else if(strcmp1(&lakon[8*i],"ESPRNGF")==0){
325 continue;
326 }else if((strcmp1(&lakon[8*i],"E")==0)&&
327 (strcmp1(&lakon[8*i+6],"1")==0)){
328 continue;
329 }else if(strcmp1(&lakon[8*i],"DCOUP3D")==0){
330 continue;
331 }else if(strcmp1(&lakon[8*i],"MASS")==0){
332 continue;
333 }else if(strcmp1(&lakon[8*i],"U1")==0){
334 continue;
335 }else{
336 indexe=ipkon[i];
337 }
338 imat=ielmat[i*mi[2]];
339 if(strcmp1(&lakon[8*i+3],"2")==0){
340
341 /* 20-node brick element */
342
343 if(((strcmp1(&lakon[8*i+6]," ")==0)||
344 (strcmp1(&filab[4],"E")==0)||
345 (strcmp1(&lakon[8*i+6],"I")==0))&&
346 (strcmp2(&lakon[8*i+6],"LC",2)!=0)){
347 if(strcmp1(output,"asc")==0){
348 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
349 m1,i+1,p4,p0,imat,m2);
350 for(j=0;j<10;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
351 fprintf(f1,"\n%3s",m2);
352 for(j=10;j<12;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
353 for(j=16;j<19;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
354 for(j=19;j<20;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
355 for(j=12;j<16;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
356 fprintf(f1,"\n");
357 }else{
358 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
359 iw=(int)ip4;fwrite(&iw,sizeof(int),1,f1);
360 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
361 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
362 for(j=0;j<10;j++){iw=(int)kon[indexe+j];
363 fwrite(&iw,sizeof(int),1,f1);}
364 for(j=10;j<12;j++){iw=(int)kon[indexe+j];
365 fwrite(&iw,sizeof(int),1,f1);}
366 for(j=16;j<19;j++){iw=(int)kon[indexe+j];
367 fwrite(&iw,sizeof(int),1,f1);}
368 for(j=19;j<20;j++){iw=(int)kon[indexe+j];
369 fwrite(&iw,sizeof(int),1,f1);}
370 for(j=12;j<16;j++){iw=(int)kon[indexe+j];
371 fwrite(&iw,sizeof(int),1,f1);}
372 }
373 }else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
374
375 /* composite material */
376
377 /* 20-node brick elements */
378
379 nlayer=0;
380 for(k=0;k<mi[2];k++){
381 if(ielmat[i*mi[2]+k]==0) break;
382 nlayer++;
383 }
384 for(k=0;k<nlayer;k++){
385 nemax++;
386 if(strcmp1(output,"asc")==0){
387 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
388 m1,nemax,p4,p0,imat,m2);
389 for(j=0;j<10;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+28+20*k+j]);
390 fprintf(f1,"\n%3s",m2);
391 for(j=10;j<12;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+28+20*k+j]);
392 for(j=16;j<19;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+28+20*k+j]);
393 for(j=19;j<20;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+28+20*k+j]);
394 for(j=12;j<16;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+28+20*k+j]);
395 fprintf(f1,"\n");
396 }else{
397 iw=(int)nemax;fwrite(&iw,sizeof(int),1,f1);
398 iw=(int)ip4;fwrite(&iw,sizeof(int),1,f1);
399 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
400 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
401 for(j=0;j<10;j++){iw=(int)kon[indexe+28+20*k+j];
402 fwrite(&iw,sizeof(int),1,f1);}
403 for(j=10;j<12;j++){iw=(int)kon[indexe+28+20*k+j];
404 fwrite(&iw,sizeof(int),1,f1);}
405 for(j=16;j<19;j++){iw=(int)kon[indexe+28+20*k+j];
406 fwrite(&iw,sizeof(int),1,f1);}
407 for(j=19;j<20;j++){iw=(int)kon[indexe+28+20*k+j];
408 fwrite(&iw,sizeof(int),1,f1);}
409 for(j=12;j<16;j++){iw=(int)kon[indexe+28+20*k+j];
410 fwrite(&iw,sizeof(int),1,f1);}
411 }
412 }
413 }else if(strcmp1(&lakon[8*i+6],"B")==0){
414
415 /* 3-node beam element */
416
417 if(strcmp1(output,"asc")==0){
418 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p12,p0,imat);
419 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe+20],
420 kon[indexe+22],kon[indexe+21]);
421 }else{
422 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
423 iw=(int)ip12;fwrite(&iw,sizeof(int),1,f1);
424 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
425 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
426 iw=(int)kon[indexe+20];fwrite(&iw,sizeof(int),1,f1);
427 iw=(int)kon[indexe+22];fwrite(&iw,sizeof(int),1,f1);
428 iw=(int)kon[indexe+21];fwrite(&iw,sizeof(int),1,f1);
429 }
430 }else{
431
432 /* 8-node 2d element */
433
434 if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
435 if(strcmp1(output,"asc")==0){
436 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
437 m1,i+1,p10,p0,imat,m2);
438 for(j=0;j<8;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+20+j]);
439 fprintf(f1,"\n");
440 }else{
441 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
442 iw=(int)ip10;fwrite(&iw,sizeof(int),1,f1);
443 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
444 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
445 for(j=0;j<8;j++){iw=(int)kon[indexe+20+j];
446 fwrite(&iw,sizeof(int),1,f1);}
447 }
448 }
449 }else if(strcmp1(&lakon[8*i+3],"8")==0){
450 if((strcmp1(&lakon[8*i+6]," ")==0)||
451 (strcmp1(&filab[4],"E")==0)){
452
453 /* 8-node brick element */
454
455 if(strcmp1(output,"asc")==0){
456 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
457 m1,i+1,p1,p0,imat,m2);
458 for(j=0;j<8;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
459 fprintf(f1,"\n");
460 }else{
461 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
462 iw=(int)ip1;fwrite(&iw,sizeof(int),1,f1);
463 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
464 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
465 for(j=0;j<8;j++){iw=(int)kon[indexe+j];
466 fwrite(&iw,sizeof(int),1,f1);}
467 }
468 }else if(strcmp1(&lakon[8*i+6],"B")==0){
469
470 /* 2-node 1d element */
471
472 if(strcmp1(&lakon[8*i+4],"R")==0){
473 if(strcmp1(output,"asc")==0){
474 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p11,p0,imat);
475 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe+8],
476 kon[indexe+9]);
477 }else{
478 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
479 iw=(int)ip11;fwrite(&iw,sizeof(int),1,f1);
480 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
481 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
482 iw=(int)kon[indexe+8];fwrite(&iw,sizeof(int),1,f1);
483 iw=(int)kon[indexe+9];fwrite(&iw,sizeof(int),1,f1);
484 }
485 }else if(strcmp1(&lakon[8*i+4],"I")==0){
486 if(strcmp1(output,"asc")==0){
487 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p11,p0,imat);
488 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe+11],
489 kon[indexe+12]);
490 }else{
491 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
492 iw=(int)ip11;fwrite(&iw,sizeof(int),1,f1);
493 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
494 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
495 iw=(int)kon[indexe+11];fwrite(&iw,sizeof(int),1,f1);
496 iw=(int)kon[indexe+12];fwrite(&iw,sizeof(int),1,f1);
497 }
498 }
499 }else{
500
501 /* 4-node 2d element */
502
503 if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
504 if((strcmp1(&lakon[8*i+4],"R")==0)||
505 (strcmp1(&lakon[8*i+4]," ")==0)){
506 if(strcmp1(output,"asc")==0){
507 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
508 m1,i+1,p9,p0,imat,m2);
509 for(j=0;j<4;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+8+j]);
510 fprintf(f1,"\n");
511 }else{
512 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
513 iw=(int)ip9;fwrite(&iw,sizeof(int),1,f1);
514 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
515 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
516 for(j=0;j<4;j++){iw=(int)kon[indexe+8+j];
517 fwrite(&iw,sizeof(int),1,f1);}
518 }
519 }else if(strcmp1(&lakon[8*i+4],"I")==0){
520 if(strcmp1(output,"asc")==0){
521 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
522 m1,i+1,p9,p0,imat,m2);
523 for(j=0;j<4;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+11+j]);
524 fprintf(f1,"\n");
525 }else{
526 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
527 iw=(int)ip9;fwrite(&iw,sizeof(int),1,f1);
528 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
529 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
530 for(j=0;j<4;j++){iw=(int)kon[indexe+11+j];
531 fwrite(&iw,sizeof(int),1,f1);}
532 }
533 }
534 }
535 }else if((strcmp1(&lakon[8*i+3],"10")==0)||
536 (strcmp1(&lakon[8*i+3],"14")==0)){
537
538 /* 10-node tetrahedral element */
539
540 if(strcmp1(output,"asc")==0){
541 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
542 m1,i+1,p6,p0,imat,m2);
543 for(j=0;j<10;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
544 fprintf(f1,"\n");
545 }else{
546 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
547 iw=(int)ip6;fwrite(&iw,sizeof(int),1,f1);
548 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
549 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
550 for(j=0;j<10;j++){iw=(int)kon[indexe+j];
551 fwrite(&iw,sizeof(int),1,f1);}
552 }
553 }else if(strcmp1(&lakon[8*i+3],"4")==0){
554
555 /* 4-node tetrahedral element */
556
557 if(strcmp1(output,"asc")==0){
558 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
559 m1,i+1,p3,p0,imat,m2);
560 for(j=0;j<4;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
561 fprintf(f1,"\n");
562 }else{
563 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
564 iw=(int)ip3;fwrite(&iw,sizeof(int),1,f1);
565 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
566 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
567 for(j=0;j<4;j++){iw=(int)kon[indexe+j];
568 fwrite(&iw,sizeof(int),1,f1);}
569 }
570 }else if(strcmp1(&lakon[8*i+3],"15")==0){
571 if(((strcmp1(&lakon[8*i+6]," ")==0)||
572 (strcmp1(&filab[4],"E")==0))&&
573 (strcmp2(&lakon[8*i+6],"LC",2)!=0)){
574
575 /* 15-node wedge element */
576
577 if(strcmp1(output,"asc")==0){
578 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
579 m1,i+1,p5,p0,imat,m2);
580 for(j=0;j<9;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
581 for(j=12;j<13;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
582 fprintf(f1,"\n%3s",m2);
583 for(j=13;j<15;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
584 for(j=9;j<12;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
585 fprintf(f1,"\n");
586 }else{
587 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
588 iw=(int)ip5;fwrite(&iw,sizeof(int),1,f1);
589 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
590 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
591 for(j=0;j<9;j++){iw=(int)kon[indexe+j];
592 fwrite(&iw,sizeof(int),1,f1);}
593 for(j=12;j<13;j++){iw=(int)kon[indexe+j];
594 fwrite(&iw,sizeof(int),1,f1);}
595 for(j=13;j<15;j++){iw=(int)kon[indexe+j];
596 fwrite(&iw,sizeof(int),1,f1);}
597 for(j=9;j<12;j++){iw=(int)kon[indexe+j];
598 fwrite(&iw,sizeof(int),1,f1);}
599 }
600
601 }else if(strcmp2(&lakon[8*i+6],"LC",2)==0){
602
603 /* composite material */
604
605 /* 15-node wedge elements */
606
607 nlayer=0;
608 for(k=0;k<mi[2];k++){
609 if(ielmat[i*mi[2]+k]==0) break;
610 nlayer++;
611 }
612 for(k=0;k<nlayer;k++){
613 nemax++;
614 if(strcmp1(output,"asc")==0){
615 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
616 m1,nemax,p5,p0,imat,m2);
617 for(j=0;j<9;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+21+15*k+j]);
618 for(j=12;j<13;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+21+15*k+j]);
619 fprintf(f1,"\n%3s",m2);
620 for(j=13;j<15;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+21+15*k+j]);
621 for(j=9;j<12;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+21+15*k+j]);
622 fprintf(f1,"\n");
623 }else{
624 iw=(int)nemax;fwrite(&iw,sizeof(int),1,f1);
625 iw=(int)ip5;fwrite(&iw,sizeof(int),1,f1);
626 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
627 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
628 for(j=0;j<9;j++){iw=(int)kon[indexe+21+15*k+j];
629 fwrite(&iw,sizeof(int),1,f1);}
630 for(j=12;j<13;j++){iw=(int)kon[indexe+21+15*k+j];
631 fwrite(&iw,sizeof(int),1,f1);}
632 for(j=13;j<15;j++){iw=(int)kon[indexe+21+15*k+j];
633 fwrite(&iw,sizeof(int),1,f1);}
634 for(j=9;j<12;j++){iw=(int)kon[indexe+21+15*k+j];
635 fwrite(&iw,sizeof(int),1,f1);}
636 }
637 }
638 }else{
639
640 /* 6-node 2d element */
641
642 if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
643 if(strcmp1(output,"asc")==0){
644 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
645 m1,i+1,p8,p0,imat,m2);
646 for(j=0;j<6;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+15+j]);
647 fprintf(f1,"\n");
648 }else{
649 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
650 iw=(int)ip8;fwrite(&iw,sizeof(int),1,f1);
651 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
652 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
653 for(j=0;j<6;j++){iw=(int)kon[indexe+15+j];
654 fwrite(&iw,sizeof(int),1,f1);}
655 }
656 }
657 }else if(strcmp1(&lakon[8*i+3],"6")==0){
658 if((strcmp1(&lakon[8*i+6]," ")==0)||
659 (strcmp1(&filab[4],"E")==0)){
660
661 /* 6-node wedge element */
662
663 if(strcmp1(output,"asc")==0){
664 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
665 m1,i+1,p2,p0,imat,m2);
666 for(j=0;j<6;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
667 fprintf(f1,"\n");
668 }else{
669 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
670 iw=(int)ip2;fwrite(&iw,sizeof(int),1,f1);
671 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
672 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
673 for(j=0;j<6;j++){iw=(int)kon[indexe+j];
674 fwrite(&iw,sizeof(int),1,f1);}
675 }
676 }else{
677
678 /* 3-node 2d element */
679
680 if(strcmp1(&lakon[8*i+6],"A")==0) iaxial=1;
681 if(strcmp1(output,"asc")==0){
682 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
683 m1,i+1,p7,p0,imat,m2);
684 for(j=0;j<3;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+6+j]);
685 fprintf(f1,"\n");
686 }else{
687 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
688 iw=(int)ip7;fwrite(&iw,sizeof(int),1,f1);
689 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
690 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
691 for(j=0;j<3;j++){iw=(int)kon[indexe+6+j];
692 fwrite(&iw,sizeof(int),1,f1);}
693 }
694 }
695 }else if((strcmp1(&lakon[8*i],"D")==0)&&(ithermal[1]>1)){
696 if(kon[indexe]==0){
697
698 /* 2-node 1d element (network entry element) */
699
700 if(strcmp1(output,"asc")==0){
701 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p11,p0,imat);
702 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,
703 kon[indexe+1],kon[indexe+2]);
704 }else{
705 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
706 iw=(int)ip11;fwrite(&iw,sizeof(int),1,f1);
707 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
708 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
709 iw=(int)kon[indexe+1];fwrite(&iw,sizeof(int),1,f1);
710 iw=(int)kon[indexe+2];fwrite(&iw,sizeof(int),1,f1);
711 }
712 }else if(kon[indexe+2]==0){
713
714 /* 2-node 1d element (network exit element) */
715
716 if(strcmp1(output,"asc")==0){
717 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p11,p0,imat);
718 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe],
719 kon[indexe+1]);
720 }else{
721 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
722 iw=(int)ip11;fwrite(&iw,sizeof(int),1,f1);
723 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
724 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
725 iw=(int)kon[indexe];fwrite(&iw,sizeof(int),1,f1);
726 iw=(int)kon[indexe+1];fwrite(&iw,sizeof(int),1,f1);
727 }
728 }else{
729
730 /* 3-node 1d element (genuine network element) */
731
732 if(strcmp1(output,"asc")==0){
733 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p12,p0,imat);
734 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe],
735 kon[indexe+2],kon[indexe+1]);
736 }else{
737 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
738 iw=(int)ip12;fwrite(&iw,sizeof(int),1,f1);
739 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
740 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
741 iw=(int)kon[indexe];fwrite(&iw,sizeof(int),1,f1);
742 iw=(int)kon[indexe+2];fwrite(&iw,sizeof(int),1,f1);
743 iw=(int)kon[indexe+1];fwrite(&iw,sizeof(int),1,f1);
744 }
745 }
746 }else if((strcmp1(&lakon[8*i],"E")==0)&&
747 ((strcmp1(&lakon[8*i+6],"A")==0)||
748 (strcmp1(&lakon[8*i+6],"2")==0))){
749
750 /* 2-node 1d element (spring element) */
751
752 if(strcmp1(output,"asc")==0){
753 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n",m1,i+1,p11,p0,imat);
754 fprintf(f1,"%3s%10" ITGFORMAT "%10" ITGFORMAT "\n",m2,kon[indexe],kon[indexe+1]);
755 }else{
756 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
757 iw=(int)ip11;fwrite(&iw,sizeof(int),1,f1);
758 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
759 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
760 iw=(int)kon[indexe];fwrite(&iw,sizeof(int),1,f1);
761 iw=(int)kon[indexe+1];fwrite(&iw,sizeof(int),1,f1);
762 }
763 // }else if(strcmp1(&lakon[8*i],"MASS")==0){
764
765 // /* MASS: store nothing */
766
767 /* user elements */
768
769 }else if(strcmp1(&lakon[8*i],"US45")==0){
770 if(strcmp1(output,"asc")==0){
771 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
772 m1,i+1,p9,p0,imat,m2);
773 for(j=0;j<4;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
774 fprintf(f1,"\n");
775 }else{
776 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
777 iw=(int)ip9;fwrite(&iw,sizeof(int),1,f1);
778 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
779 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
780 for(j=0;j<4;j++){iw=(int)kon[indexe+j];
781 fwrite(&iw,sizeof(int),1,f1);}
782 }
783
784 }else if(strcmp1(&lakon[8*i],"US3")==0){
785 if(strcmp1(output,"asc")==0){
786 fprintf(f1,"%3s%10" ITGFORMAT "%5s%5s%5" ITGFORMAT "\n%3s",
787 m1,i+1,p7,p0,imat,m2);
788 for(j=0;j<3;j++)fprintf(f1,"%10" ITGFORMAT "",kon[indexe+j]);
789 fprintf(f1,"\n");
790 }else{
791 iw=(int)(i+1);fwrite(&iw,sizeof(int),1,f1);
792 iw=(int)ip7;fwrite(&iw,sizeof(int),1,f1);
793 iw=(int)ip0;fwrite(&iw,sizeof(int),1,f1);
794 iw=(int)imat;fwrite(&iw,sizeof(int),1,f1);
795 for(j=0;j<3;j++){iw=(int)kon[indexe+j];
796 fwrite(&iw,sizeof(int),1,f1);}
797 }
798 }else{
799
800 /* not treated element type: may lead to an inconsistency
801 in the element count and element output, which may
802 cause a crash while reading a binary output file */
803
804 FORTRAN(writeelem,(&i,lakon));
805
806 }
807 }
808 if(strcmp1(output,"asc")==0)fprintf(f1,"%3s\n",m3);
809
810 if(*nmethod==0){fclose(f1);return;}
811 }
812
813 /* for cyclic symmetry frequency calculations only results for
814 even numbers (= odd modes, numbering starts at 0) are stored */
815
816 if((*nmethod==2)&&(((*mode/2)*2!=*mode)&&(*noddiam>=0))){fclose(f1);return;}
817
818 /* storing the displacements in the nodes */
819
820 if((*nmethod!=5)||(*mode==-1)){
821 if((strcmp1(filab,"U ")==0)&&(*ithermal!=2)){
822 iselect=1;
823
824 frdset(filab,set,&iset,istartset,iendset,ialset,
825 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
826 ngraph);
827
828 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
829 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
830
831 if(mi[1]==3){
832
833 fprintf(f1," -4 DISP 4 1\n");
834 fprintf(f1," -5 D1 1 2 1 0\n");
835 fprintf(f1," -5 D2 1 2 2 0\n");
836 fprintf(f1," -5 D3 1 2 3 0\n");
837 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
838
839 frdvector(v,&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
840 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
841
842 }else if((mi[1]>3)&&(mi[1]<7)){
843
844 fprintf(f1," -4 DISP %1" ITGFORMAT " 1\n",mi[1]+1);
845 fprintf(f1," -5 D1 1 2 1 0\n");
846 fprintf(f1," -5 D2 1 2 2 0\n");
847 fprintf(f1," -5 D3 1 2 3 0\n");
848 for(j=4;j<=mi[1];j++){
849 fprintf(f1," -5 D%1" ITGFORMAT " 1 1 0 0\n",j);
850 }
851 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
852
853 frdgeneralvector(v,&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
854 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
855 }else{
856 printf("*WARNING in frd:\n");
857 printf(" for output purposes only 4, 5 or 6\n");
858 printf(" degrees of freedom are allowed\n");
859 printf(" for generalized vectors;\n");
860 printf(" actual degrees of freedom = %" ITGFORMAT "\n",mi[1]);
861 printf(" output request ist not performed;\n");
862 }
863 }
864 }
865
866 /* storing the imaginary part of displacements in the nodes
867 for the odd modes of cyclic symmetry calculations */
868
869 if(*noddiam>=0){
870 if((strcmp1(filab,"U ")==0)&&(*ithermal!=2)){
871
872 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
873 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
874
875 fprintf(f1," -4 DISPI 4 1\n");
876 fprintf(f1," -5 D1 1 2 1 0\n");
877 fprintf(f1," -5 D2 1 2 2 0\n");
878 fprintf(f1," -5 D3 1 2 3 0\n");
879 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
880
881 frdvector(&v[*nk*mt],&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
882 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
883 }
884 }
885
886 /* storing the imaginary part of displacements in the nodes
887 for steady state calculations */
888
889 if((*nmethod==5)&&(*mode==0)){
890 if((strcmp1(filab,"U ")==0)&&(*ithermal!=2)){
891 iselect=1;
892
893 frdset(filab,set,&iset,istartset,iendset,ialset,
894 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
895 ngraph);
896
897 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
898 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
899
900 fprintf(f1," -4 DISPI 4 1\n");
901 fprintf(f1," -5 D1 1 2 1 0\n");
902 fprintf(f1," -5 D2 1 2 2 0\n");
903 fprintf(f1," -5 D3 1 2 3 0\n");
904 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
905
906 frdvector(v,&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
907 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
908 }
909 }
910
911 /* storing the velocities in the nodes */
912
913 if((strcmp1(&filab[1740],"V ")==0)&&(*ithermal!=2)){
914 iselect=1;
915
916 frdset(&filab[1740],set,&iset,istartset,iendset,ialset,
917 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
918 ngraph);
919
920 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
921 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
922
923 fprintf(f1," -4 VELO 4 1\n");
924 fprintf(f1," -5 V1 1 2 1 0\n");
925 fprintf(f1," -5 V2 1 2 2 0\n");
926 fprintf(f1," -5 V3 1 2 3 0\n");
927 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
928
929 frdvector(veold,&iset,ntrans,&filab[1740],&nkcoords,inum,m1,inotr,
930 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
931 }
932
933 /* storing the temperatures in the nodes */
934
935 if(strcmp1(&filab[87],"NT ")==0){
936 iselect=0;
937
938 frdset(&filab[87],set,&iset,istartset,iendset,ialset,
939 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
940 ngraph);
941
942 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
943 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
944
945 fprintf(f1," -4 NDTEMP 1 1\n");
946 fprintf(f1," -5 T 1 1 0 0\n");
947
948 if(*ithermal<=1){
949 frdselect(t1,t1,&iset,&nkcoords,inum,m1,istartset,iendset,
950 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
951 nfieldscalar,&iselect,m2,f1,output,m3);
952 }else{
953 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
954 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
955 nfieldvector0,&iselect,m2,f1,output,m3);
956 }
957 }
958
959 /* storing the electrical potential in the nodes */
960
961 if((strcmp1(&filab[3654],"POT ")==0)&&(*ithermal==2)){
962 iselect=0;
963
964 frdset(&filab[3654],set,&iset,istartset,iendset,ialset,
965 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
966 ngraph);
967
968 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
969 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
970
971 fprintf(f1," -4 ELPOT 1 1\n");
972 fprintf(f1," -5 V 1 1 0 0\n");
973
974 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
975 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
976 nfieldvector0,&iselect,m2,f1,output,m3);
977 }
978
979 /* storing the stresses in the nodes */
980
981 if((*nmethod!=5)||(*mode==-1)){
982 if((strcmp1(&filab[174],"S ")==0)&&(*ithermal!=2)){
983 iselect=1;
984
985 frdset(&filab[174],set,&iset,istartset,iendset,ialset,
986 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
987 ngraph);
988
989 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
990 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
991
992 fprintf(f1," -4 STRESS 6 1\n");
993 fprintf(f1," -5 SXX 1 4 1 1\n");
994 fprintf(f1," -5 SYY 1 4 2 2\n");
995 fprintf(f1," -5 SZZ 1 4 3 3\n");
996 fprintf(f1," -5 SXY 1 4 1 2\n");
997 fprintf(f1," -5 SYZ 1 4 2 3\n");
998 fprintf(f1," -5 SZX 1 4 3 1\n");
999
1000 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1001 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1002 nfieldtensor,&iselect,m2,f1,output,m3);
1003
1004 }
1005 }
1006
1007 if((*nmethod!=5)||(*mode==-1)){
1008 if((strcmp1(&filab[4350],"SNEG")==0)&&(*ithermal!=2)){
1009 iselect=1;
1010
1011 frdset(&filab[4350],set,&iset,istartset,iendset,ialset,
1012 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1013 ngraph);
1014
1015 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1016 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1017
1018 fprintf(f1," -4 STRNEG 6 1\n");
1019 fprintf(f1," -5 SXX 1 4 1 1\n");
1020 fprintf(f1," -5 SYY 1 4 2 2\n");
1021 fprintf(f1," -5 SZZ 1 4 3 3\n");
1022 fprintf(f1," -5 SXY 1 4 1 2\n");
1023 fprintf(f1," -5 SYZ 1 4 2 3\n");
1024 fprintf(f1," -5 SZX 1 4 3 1\n");
1025
1026 iflag=-1;
1027 numfield=6;
1028 if((norien>0)&&(strcmp1(&filab[4355],"L")==0)){
1029 iorienloc=1;
1030 }else{
1031 iorienloc=0;
1032 }
1033 NNEW(inumshell,ITG,*nk);
1034 FORTRAN(extrapolateshell,(stx,stn,ipkon,inumshell,kon,lakon,&numfield,
1035 nk,ne,mi,&numfield,orab,ielorien,co,
1036 &iorienloc,&filab[4354],ielmat,
1037 thicke,ielprop,prop,&iflag));
1038 SFREE(inumshell);
1039
1040 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1041 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1042 nfieldtensor,&iselect,m2,f1,output,m3);
1043
1044 }
1045 }
1046
1047 if((*nmethod!=5)||(*mode==-1)){
1048 if((strcmp1(&filab[4437],"SMID")==0)&&(*ithermal!=2)){
1049 iselect=1;
1050
1051 frdset(&filab[4437],set,&iset,istartset,iendset,ialset,
1052 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1053 ngraph);
1054
1055 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1056 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1057
1058 fprintf(f1," -4 STRMID 6 1\n");
1059 fprintf(f1," -5 SXX 1 4 1 1\n");
1060 fprintf(f1," -5 SYY 1 4 2 2\n");
1061 fprintf(f1," -5 SZZ 1 4 3 3\n");
1062 fprintf(f1," -5 SXY 1 4 1 2\n");
1063 fprintf(f1," -5 SYZ 1 4 2 3\n");
1064 fprintf(f1," -5 SZX 1 4 3 1\n");
1065
1066 iflag=0;
1067 numfield=6;
1068 if((norien>0)&&(strcmp1(&filab[4442],"L")==0)){
1069 iorienloc=1;
1070 }else{
1071 iorienloc=0;
1072 }
1073 NNEW(inumshell,ITG,*nk);
1074 FORTRAN(extrapolateshell,(stx,stn,ipkon,inumshell,kon,lakon,&numfield,
1075 nk,ne,mi,&numfield,orab,ielorien,co,
1076 &iorienloc,&filab[4441],ielmat,
1077 thicke,ielprop,prop,&iflag));
1078 SFREE(inumshell);
1079
1080 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1081 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1082 nfieldtensor,&iselect,m2,f1,output,m3);
1083
1084 }
1085 }
1086
1087 if((*nmethod!=5)||(*mode==-1)){
1088 if((strcmp1(&filab[4524],"SPOS")==0)&&(*ithermal!=2)){
1089 iselect=1;
1090
1091 frdset(&filab[4524],set,&iset,istartset,iendset,ialset,
1092 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1093 ngraph);
1094
1095 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1096 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1097
1098 fprintf(f1," -4 STRPOS 6 1\n");
1099 fprintf(f1," -5 SXX 1 4 1 1\n");
1100 fprintf(f1," -5 SYY 1 4 2 2\n");
1101 fprintf(f1," -5 SZZ 1 4 3 3\n");
1102 fprintf(f1," -5 SXY 1 4 1 2\n");
1103 fprintf(f1," -5 SYZ 1 4 2 3\n");
1104 fprintf(f1," -5 SZX 1 4 3 1\n");
1105
1106 iflag=1;
1107 numfield=6;
1108 if((norien>0)&&(strcmp1(&filab[4529],"L")==0)){
1109 iorienloc=1;
1110 }else{
1111 iorienloc=0;
1112 }
1113 NNEW(inumshell,ITG,*nk);
1114 FORTRAN(extrapolateshell,(stx,stn,ipkon,inumshell,kon,lakon,&numfield,
1115 nk,ne,mi,&numfield,orab,ielorien,co,
1116 &iorienloc,&filab[4528],ielmat,
1117 thicke,ielprop,prop,&iflag));
1118 SFREE(inumshell);
1119
1120 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1121 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1122 nfieldtensor,&iselect,m2,f1,output,m3);
1123
1124 }
1125 }
1126
1127 /* storing the imaginary part of the stresses in the nodes
1128 for the odd modes of cyclic symmetry calculations */
1129
1130 if(*noddiam>=0){
1131 if((strcmp1(&filab[174],"S ")==0)&&(*ithermal!=2)){
1132
1133 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1134 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1135
1136 fprintf(f1," -4 STRESSI 6 1\n");
1137 fprintf(f1," -5 SXX 1 4 1 1\n");
1138 fprintf(f1," -5 SYY 1 4 2 2\n");
1139 fprintf(f1," -5 SZZ 1 4 3 3\n");
1140 fprintf(f1," -5 SXY 1 4 1 2\n");
1141 fprintf(f1," -5 SYZ 1 4 2 3\n");
1142 fprintf(f1," -5 SZX 1 4 3 1\n");
1143
1144 frdselect(&stn[6**nk],stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1145 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1146 nfieldtensor,&iselect,m2,f1,output,m3);
1147
1148 }
1149 }
1150
1151 /* storing the imaginary part of the stresses in the nodes
1152 for steady state calculations */
1153
1154 if((*nmethod==5)&&(*mode==0)){
1155 if((strcmp1(&filab[174],"S ")==0)&&(*ithermal!=2)){
1156 iselect=1;
1157
1158 frdset(&filab[174],set,&iset,istartset,iendset,ialset,
1159 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1160 ngraph);
1161
1162 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1163 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1164
1165 fprintf(f1," -4 STRESSI 6 1\n");
1166 fprintf(f1," -5 SXX 1 4 1 1\n");
1167 fprintf(f1," -5 SYY 1 4 2 2\n");
1168 fprintf(f1," -5 SZZ 1 4 3 3\n");
1169 fprintf(f1," -5 SXY 1 4 1 2\n");
1170 fprintf(f1," -5 SYZ 1 4 2 3\n");
1171 fprintf(f1," -5 SZX 1 4 3 1\n");
1172
1173 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1174 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1175 nfieldtensor,&iselect,m2,f1,output,m3);
1176
1177 }
1178 }
1179
1180 /* storing the electromagnetic field E in the nodes */
1181
1182 if((strcmp1(&filab[3741],"EMFE")==0)&&(*ithermal!=2)){
1183 iselect=1;
1184
1185 frdset(&filab[3741],set,&iset,istartset,iendset,ialset,
1186 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1187 ngraph);
1188
1189 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1190 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1191
1192 fprintf(f1," -4 EMFE 4 1\n");
1193 fprintf(f1," -5 E1 1 2 1 0\n");
1194 fprintf(f1," -5 E2 1 2 2 0\n");
1195 fprintf(f1," -5 E3 1 2 3 0\n");
1196 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1197
1198 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1199 ialset,ngraph,&ncompvector,ifieldvector,icompvector,
1200 nfieldtensor,&iselect,m2,f1,output,m3);
1201
1202 if(*nmethod==2){
1203 iselect=1;
1204
1205 frdset(&filab[3741],set,&iset,istartset,iendset,ialset,
1206 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1207 ngraph);
1208
1209 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1210 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1211
1212 fprintf(f1," -4 EMFEI 4 1\n");
1213 fprintf(f1," -5 E1 1 2 1 0\n");
1214 fprintf(f1," -5 E2 1 2 2 0\n");
1215 fprintf(f1," -5 E3 1 2 3 0\n");
1216 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1217
1218 frdselect(&stn[6**nk],stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1219 ialset,ngraph,&ncompvector,ifieldvector,icompvector,
1220 nfieldtensor,&iselect,m2,f1,output,m3);
1221 }
1222 }
1223
1224 /* storing the electromagnetic field B in the nodes */
1225
1226 if((strcmp1(&filab[3828],"EMFB")==0)&&(*ithermal!=2)){
1227 iselect=1;
1228
1229 frdset(&filab[3828],set,&iset,istartset,iendset,ialset,
1230 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1231 ngraph);
1232
1233 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1234 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1235
1236 fprintf(f1," -4 EMFB 4 1\n");
1237 fprintf(f1," -5 B1 1 2 1 0\n");
1238 fprintf(f1," -5 B2 1 2 2 0\n");
1239 fprintf(f1," -5 B3 1 2 3 0\n");
1240 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1241
1242 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1243 ialset,ngraph,&ncompvector,ifieldvector,icompvectorlast,
1244 nfieldtensor,&iselect,m2,f1,output,m3);
1245
1246 if(*nmethod==2){
1247 iselect=1;
1248
1249 frdset(&filab[3828],set,&iset,istartset,iendset,ialset,
1250 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1251 ngraph);
1252
1253 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1254 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1255
1256 fprintf(f1," -4 EMFBI 4 1\n");
1257 fprintf(f1," -5 B1 1 2 1 0\n");
1258 fprintf(f1," -5 B2 1 2 2 0\n");
1259 fprintf(f1," -5 B3 1 2 3 0\n");
1260 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1261
1262 frdselect(&stn[6**nk],stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1263 ialset,ngraph,&ncompvector,ifieldvector,icompvectorlast,
1264 nfieldtensor,&iselect,m2,f1,output,m3);
1265 }
1266
1267 }
1268
1269 /* storing the total strains in the nodes */
1270
1271 if((*nmethod!=5)||(*mode==-1)){
1272 if((strcmp1(&filab[261],"E ")==0)&&(*ithermal!=2)){
1273 iselect=1;
1274
1275 frdset(&filab[261],set,&iset,istartset,iendset,ialset,
1276 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1277 ngraph);
1278
1279 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1280 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1281
1282 fprintf(f1," -4 TOSTRAIN 6 1\n");
1283 fprintf(f1," -5 EXX 1 4 1 1\n");
1284 fprintf(f1," -5 EYY 1 4 2 2\n");
1285 fprintf(f1," -5 EZZ 1 4 3 3\n");
1286 fprintf(f1," -5 EXY 1 4 1 2\n");
1287 fprintf(f1," -5 EYZ 1 4 2 3\n");
1288 fprintf(f1," -5 EZX 1 4 3 1\n");
1289
1290 frdselect(een,een,&iset,&nkcoords,inum,m1,istartset,iendset,
1291 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1292 nfieldtensor,&iselect,m2,f1,output,m3);
1293
1294 }
1295 }
1296
1297 /* storing the imaginary part of the total strains in the nodes
1298 for the odd modes of cyclic symmetry calculations */
1299
1300 if(*noddiam>=0){
1301 if((strcmp1(&filab[261],"E ")==0)&&(*ithermal!=2)){
1302
1303 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1304 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1305
1306 fprintf(f1," -4 TOSTRAII 6 1\n");
1307 fprintf(f1," -5 EXX 1 4 1 1\n");
1308 fprintf(f1," -5 EYY 1 4 2 2\n");
1309 fprintf(f1," -5 EZZ 1 4 3 3\n");
1310 fprintf(f1," -5 EXY 1 4 1 2\n");
1311 fprintf(f1," -5 EYZ 1 4 2 3\n");
1312 fprintf(f1," -5 EZX 1 4 3 1\n");
1313
1314 frdselect(&een[6**nk],een,&iset,&nkcoords,inum,m1,istartset,iendset,
1315 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1316 nfieldtensor,&iselect,m2,f1,output,m3);
1317
1318 }
1319 }
1320
1321 /* storing the imaginary part of the total strains in the nodes
1322 for steady state calculations */
1323
1324 if((*nmethod==5)&&(*mode==0)){
1325 if((strcmp1(&filab[261],"E ")==0)&&(*ithermal!=2)){
1326 iselect=1;
1327
1328 frdset(&filab[261],set,&iset,istartset,iendset,ialset,
1329 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1330 ngraph);
1331
1332 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1333 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1334
1335 fprintf(f1," -4 TOSTRAII 6 1\n");
1336 fprintf(f1," -5 EXX 1 4 1 1\n");
1337 fprintf(f1," -5 EYY 1 4 2 2\n");
1338 fprintf(f1," -5 EZZ 1 4 3 3\n");
1339 fprintf(f1," -5 EXY 1 4 1 2\n");
1340 fprintf(f1," -5 EYZ 1 4 2 3\n");
1341 fprintf(f1," -5 EZX 1 4 3 1\n");
1342
1343 frdselect(een,een,&iset,&nkcoords,inum,m1,istartset,iendset,
1344 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1345 nfieldtensor,&iselect,m2,f1,output,m3);
1346
1347 }
1348 }
1349
1350 /* storing the mechanical strains in the nodes */
1351
1352 if((*nmethod!=5)||(*mode==-1)){
1353 if((strcmp1(&filab[2697],"ME ")==0)&&(*ithermal!=2)){
1354 iselect=1;
1355
1356 frdset(&filab[2697],set,&iset,istartset,iendset,ialset,
1357 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1358 ngraph);
1359
1360 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1361 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1362
1363 fprintf(f1," -4 MESTRAIN 6 1\n");
1364 fprintf(f1," -5 MEXX 1 4 1 1\n");
1365 fprintf(f1," -5 MEYY 1 4 2 2\n");
1366 fprintf(f1," -5 MEZZ 1 4 3 3\n");
1367 fprintf(f1," -5 MEXY 1 4 1 2\n");
1368 fprintf(f1," -5 MEYZ 1 4 2 3\n");
1369 fprintf(f1," -5 MEZX 1 4 3 1\n");
1370
1371
1372 frdselect(emn,emn,&iset,&nkcoords,inum,m1,istartset,iendset,
1373 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1374 nfieldtensor,&iselect,m2,f1,output,m3);
1375
1376 }
1377 }
1378
1379 /* storing the imaginary part of the mechanical strains in the nodes
1380 for the odd modes of cyclic symmetry calculations */
1381
1382 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1383 if((strcmp1(&filab[2697],"ME ")==0)&&(*ithermal!=2)){
1384
1385 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1386 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1387
1388 fprintf(f1," -4 MESTRAII 6 1\n");
1389 fprintf(f1," -5 MEXX 1 4 1 1\n");
1390 fprintf(f1," -5 MEYY 1 4 2 2\n");
1391 fprintf(f1," -5 MEZZ 1 4 3 3\n");
1392 fprintf(f1," -5 MEXY 1 4 1 2\n");
1393 fprintf(f1," -5 MEYZ 1 4 2 3\n");
1394 fprintf(f1," -5 MEZX 1 4 3 1\n");
1395
1396 frdselect(&emn[6**nk],een,&iset,&nkcoords,inum,m1,istartset,iendset,
1397 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1398 nfieldtensor,&iselect,m2,f1,output,m3);
1399
1400 }
1401 }
1402
1403 /* storing the thermal strains in the nodes */
1404
1405 if((*nmethod!=5)||(*mode==-1)){
1406 if((strcmp1(&filab[4698],"THE ")==0)&&(*ithermal!=2)){
1407 iselect=1;
1408
1409 frdset(&filab[4698],set,&iset,istartset,iendset,ialset,
1410 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1411 ngraph);
1412
1413 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1414 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1415
1416 fprintf(f1," -4 THSTRAIN 6 1\n");
1417 fprintf(f1," -5 THXX 1 4 1 1\n");
1418 fprintf(f1," -5 THYY 1 4 2 2\n");
1419 fprintf(f1," -5 THZZ 1 4 3 3\n");
1420 fprintf(f1," -5 THXY 1 4 1 2\n");
1421 fprintf(f1," -5 THYZ 1 4 2 3\n");
1422 fprintf(f1," -5 THZX 1 4 3 1\n");
1423
1424 NNEW(ethn,double,6**nk);
1425 for(i=0;i<6**nk;i++){
1426 ethn[i]=een[i]-emn[i];
1427 }
1428
1429 frdselect(ethn,ethn,&iset,&nkcoords,inum,m1,istartset,iendset,
1430 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1431 nfieldtensor,&iselect,m2,f1,output,m3);
1432
1433 SFREE(ethn);
1434
1435 }
1436 }
1437
1438 /* storing the forces in the nodes */
1439
1440 if((*nmethod!=5)||(*mode==-1)){
1441 if((strcmp1(&filab[348],"RF ")==0)&&(*ithermal!=2)){
1442 iselect=1;
1443
1444 frdset(&filab[348],set,&iset,istartset,iendset,ialset,
1445 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1446 ngraph);
1447
1448 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1449 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1450
1451 if(mi[1]==3){
1452
1453 fprintf(f1," -4 FORC 4 1\n");
1454 fprintf(f1," -5 F1 1 2 1 0\n");
1455 fprintf(f1," -5 F2 1 2 2 0\n");
1456 fprintf(f1," -5 F3 1 2 3 0\n");
1457 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1458
1459 if((iaxial==1)&&(strcmp1(&filab[352],"I")==0)){for(i=0;i<*nk;i++){fn[1+i*mt]*=180.;fn[2+i*mt]*=180.;fn[3+i*mt]*=180.;}}
1460 frdvector(fn,&iset,ntrans,&filab[348],&nkcoords,inum,m1,inotr,
1461 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
1462 if((iaxial==1)&&(strcmp1(&filab[352],"I")==0)){for(i=0;i<*nk;i++){fn[1+i*mt]/=180.;fn[2+i*mt]/=180.;fn[3+i*mt]/=180.;}}
1463
1464 }else if((mi[1]>3)&&(mi[1]<7)){
1465
1466 fprintf(f1," -4 FORC %1" ITGFORMAT " 1\n",mi[1]+1);
1467 fprintf(f1," -5 F1 1 2 1 0\n");
1468 fprintf(f1," -5 F2 1 2 2 0\n");
1469 fprintf(f1," -5 F3 1 2 3 0\n");
1470 for(j=4;j<=mi[1];j++){
1471 fprintf(f1," -5 F%1" ITGFORMAT " 1 1 0 0\n",j);
1472 }
1473 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1474
1475 frdgeneralvector(fn,&iset,ntrans,&filab[348],&nkcoords,inum,m1,inotr,
1476 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
1477 }else{
1478 printf("*WARNING in frd:\n");
1479 printf(" for output purposes only 4, 5 or 6\n");
1480 printf(" degrees of freedom are allowed\n");
1481 printf(" for generalized vectors;\n");
1482 printf(" actual degrees of freedom = %" ITGFORMAT "\n",mi[1]);
1483 printf(" output request ist not performed;\n");
1484 }
1485 }
1486 }
1487
1488 /* storing the imaginary part of the forces in the nodes
1489 for the odd modes of cyclic symmetry calculations */
1490
1491 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1492 if((strcmp1(&filab[348],"RF ")==0)&&(*ithermal!=2)){
1493
1494 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1495 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1496
1497 fprintf(f1," -4 FORCI 4 1\n");
1498 fprintf(f1," -5 F1 1 2 1 0\n");
1499 fprintf(f1," -5 F2 1 2 2 0\n");
1500 fprintf(f1," -5 F3 1 2 3 0\n");
1501 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1502
1503 frdvector(&fn[*nk*mt],&iset,ntrans,filab,&nkcoords,inum,m1,inotr,
1504 trab,co,istartset,iendset,ialset,mi,ngraph,f1,output,m3);
1505 }
1506 }
1507
1508 /* storing the equivalent plastic strains in the nodes */
1509
1510 if((strcmp1(&filab[435],"PEEQ")==0)&&(*ithermal!=2)){
1511 iselect=1;
1512
1513 frdset(&filab[435],set,&iset,istartset,iendset,ialset,
1514 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1515 ngraph);
1516
1517 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1518 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1519
1520 fprintf(f1," -4 PE 1 1\n");
1521 fprintf(f1," -5 PE 1 1 0 0\n");
1522
1523 frdselect(epn,epn,&iset,&nkcoords,inum,m1,istartset,iendset,
1524 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1525 nfieldscalar,&iselect,m2,f1,output,m3);
1526
1527 }
1528
1529 /* storing the energy in the nodes */
1530
1531 if((*nmethod!=5)||(*mode==-1)){
1532 if((strcmp1(&filab[522],"ENER")==0)&&(*ithermal!=2)){
1533 iselect=1;
1534
1535 frdset(&filab[522],set,&iset,istartset,iendset,ialset,
1536 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1537 ngraph);
1538
1539 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1540 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1541
1542 fprintf(f1," -4 ENER 1 1\n");
1543 fprintf(f1," -5 ENER 1 1 0 0\n");
1544
1545 frdselect(enern,enern,&iset,&nkcoords,inum,m1,istartset,iendset,
1546 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1547 nfieldscalar,&iselect,m2,f1,output,m3);
1548
1549 }
1550 }
1551
1552 /* storing the contact displacements and stresses at the slave nodes */
1553
1554 /* node-to-face penalty */
1555
1556 if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar!=1)&&(*ithermal!=2)&&((*nmethod!=2)&&(*nmethod!=13))){
1557
1558 for(i=*ne-1;i>=0;i--){
1559 if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
1560 break;
1561 }
1562 noutloc=*ne-i-1;
1563
1564 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1565 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1566
1567 fprintf(f1," -4 CONTACT 6 1\n");
1568 fprintf(f1," -5 COPEN 1 4 1 1\n");
1569 fprintf(f1," -5 CSLIP1 1 4 2 2\n");
1570 fprintf(f1," -5 CSLIP2 1 4 3 3\n");
1571 fprintf(f1," -5 CPRESS 1 4 1 2\n");
1572 fprintf(f1," -5 CSHEAR1 1 4 2 3\n");
1573 fprintf(f1," -5 CSHEAR2 1 4 3 1\n");
1574
1575 for(i=*ne-1;i>=0;i--){
1576 if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
1577 break;
1578 strcpy1(text,&lakon[8*i+7],1);
1579 nope=atoi(text)+1;
1580 nodes=kon[ipkon[i]+nope-1];
1581 if(strcmp1(output,"asc")==0){
1582 fprintf(f1,"%3s%10" ITGFORMAT "",m1,nodes);
1583 for(j=0;j<6;j++)fprintf(f1,"%12.5E",(float)stx[6*mi[0]*i+j]);
1584 }else{
1585 iw=(int)(nodes);fwrite(&iw,sizeof(int),1,f1);
1586 if(strcmp1(output,"bin")==0){
1587 for(j=0;j<6;j++){
1588 fl=(float)stx[6*mi[0]*i+j];
1589 fwrite(&fl,sizeof(float),1,f1);
1590 }
1591 }else{
1592 for(j=0;j<6;j++){
1593 fwrite(&stx[6*mi[0]*i+j],sizeof(double),1,f1);
1594 }
1595 }
1596 }
1597 if(strcmp1(output,"asc")==0)fprintf(f1,"\n");
1598 }
1599
1600 if(strcmp1(output,"asc")==0)fprintf(f1,"%3s\n",m3);
1601 }
1602
1603 /* face-to-face penalty */
1604
1605 if((*nmethod!=5)||(*mode==-1)){
1606 if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==1)&&(*ithermal!=2)){
1607 iselect=1;
1608
1609 frdset(&filab[2175],set,&iset,istartset,iendset,ialset,
1610 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1611 ngraph);
1612
1613 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1614 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1615 fprintf(f1," -4 CONTACT 6 1\n");
1616 fprintf(f1," -5 COPEN 1 4 1 1\n");
1617 fprintf(f1," -5 CSLIP1 1 4 2 2\n");
1618 fprintf(f1," -5 CSLIP2 1 4 3 3\n");
1619 fprintf(f1," -5 CPRESS 1 4 1 2\n");
1620 fprintf(f1," -5 CSHEAR1 1 4 2 3\n");
1621 fprintf(f1," -5 CSHEAR2 1 4 3 1\n");
1622
1623 frdselect(cdn,cdn,&iset,&nkcoords,inum,m1,istartset,iendset,
1624 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1625 nfieldtensor,&iselect,m2,f1,output,m3);
1626
1627 }
1628 }
1629
1630 /* storing imaginary part of the differential contact displacements
1631 and the contact stresses for the odd modes of cyclic symmetry
1632 calculations */
1633
1634 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1635 if((strcmp1(&filab[2175],"CONT")==0)&&(*mortar==1)&&(*ithermal!=2)){
1636 iselect=1;
1637
1638 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1639 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1640 fprintf(f1," -4 CONTACTI 6 1\n");
1641 fprintf(f1," -5 COPEN 1 4 1 1\n");
1642 fprintf(f1," -5 CSLIP1 1 4 2 2\n");
1643 fprintf(f1," -5 CSLIP2 1 4 3 3\n");
1644 fprintf(f1," -5 CPRESS 1 4 1 2\n");
1645 fprintf(f1," -5 CSHEAR1 1 4 2 3\n");
1646 fprintf(f1," -5 CSHEAR2 1 4 3 1\n");
1647
1648 frdselect(&cdn[6**nk],cdn,&iset,&nkcoords,inum,m1,istartset,iendset,
1649 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1650 nfieldtensor,&iselect,m2,f1,output,m3);
1651
1652 }
1653 }
1654 /* storing the contact energy at the slave nodes */
1655
1656 if((strcmp1(&filab[2262],"CELS")==0)&&(*ithermal!=2)){
1657
1658 for(i=*ne-1;i>=0;i--){
1659 if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
1660 break;
1661 }
1662 noutloc=*ne-i-1;
1663
1664 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1665 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1666
1667 fprintf(f1," -4 CELS 1 1\n");
1668 fprintf(f1," -5 CELS 1 1 0 0\n");
1669
1670 for(i=*ne-1;i>=0;i--){
1671 if((strcmp1(&lakon[8*i+1],"S")!=0)||(strcmp1(&lakon[8*i+6],"C")!=0))
1672 break;
1673 nope=atoi(&lakon[8*i+7])+1;
1674 nodes=kon[ipkon[i]+nope-1];
1675 if(strcmp1(output,"asc")==0){
1676 fprintf(f1,"%3s%10" ITGFORMAT "%12.5E\n",m1,nodes,(float)ener[i*mi[0]]);
1677 }else{
1678 iw=(int)(nodes);fwrite(&iw,sizeof(int),1,f1);
1679 if(strcmp1(output,"bin")==0){
1680 fl=(float)ener[i*mi[0]];
1681 fwrite(&fl,sizeof(float),1,f1);
1682 }else{
1683 fwrite(&ener[i*mi[0]],sizeof(double),1,f1);
1684 }
1685 }
1686 }
1687
1688 if(strcmp1(output,"asc")==0)fprintf(f1,"%3s\n",m3);
1689 }
1690
1691 /* storing the internal state variables in the nodes */
1692
1693 if(strcmp1(&filab[609],"SDV ")==0){
1694 iselect=1;
1695
1696 frdset(&filab[609],set,&iset,istartset,iendset,ialset,
1697 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1698 ngraph);
1699
1700 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1701 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1702
1703 fprintf(f1," -4 SDV %3" ITGFORMAT " 1\n",*nstate_);
1704 for(j=1;j<=*nstate_;j++){
1705 fprintf(f1," -5 SDV%-3" ITGFORMAT " 1 1 0 0\n",j);
1706 }
1707
1708 for(i=0;i<*nstate_;i++){
1709 ifieldstate[i]=1;icompstate[i]=i;
1710 }
1711 nfield[0]=*nstate_;
1712
1713 frdselect(xstaten,xstaten,&iset,&nkcoords,inum,m1,istartset,iendset,
1714 ialset,ngraph,nstate_,ifieldstate,icompstate,
1715 nfield,&iselect,m2,f1,output,m3);
1716
1717 }
1718
1719 /* storing the heat flux in the nodes
1720 the heat flux has been extrapolated from the integration points
1721 in subroutine extrapolate.f, taking into account whether the
1722 results are requested in the global system or in a local system.
1723 Therefore, subroutine frdvector cannot be used, since it assumes
1724 the values are stored in the global system */
1725
1726 if((strcmp1(&filab[696],"HFL ")==0)&&(*ithermal>1)){
1727 iselect=1;
1728
1729 frdset(&filab[696],set,&iset,istartset,iendset,ialset,
1730 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1731 ngraph);
1732
1733 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1734 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1735
1736 fprintf(f1," -4 FLUX 4 1\n");
1737 fprintf(f1," -5 F1 1 2 1 0\n");
1738 fprintf(f1," -5 F2 1 2 2 0\n");
1739 fprintf(f1," -5 F3 1 2 3 0\n");
1740 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1741
1742 frdselect(qfn,qfn,&iset,&nkcoords,inum,m1,istartset,iendset,
1743 ialset,ngraph,&ncompvector,ifieldvector,icompvector,
1744 nfieldvector1,&iselect,m2,f1,output,m3);
1745
1746 }
1747
1748 /* storing the electrical current in the nodes
1749 (cf. heat flux HFL above) */
1750
1751 if((strcmp1(&filab[3567],"ECD ")==0)&&(*ithermal==2)){
1752 iselect=1;
1753
1754 frdset(&filab[3567],set,&iset,istartset,iendset,ialset,
1755 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1756 ngraph);
1757
1758 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1759 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1760
1761 fprintf(f1," -4 CURR 4 1\n");
1762 fprintf(f1," -5 j1 1 2 1 0\n");
1763 fprintf(f1," -5 j2 1 2 2 0\n");
1764 fprintf(f1," -5 j3 1 2 3 0\n");
1765 fprintf(f1," -5 ALL 1 2 0 0 1ALL\n");
1766
1767 frdselect(qfn,qfn,&iset,&nkcoords,inum,m1,istartset,iendset,
1768 ialset,ngraph,&ncompvector,ifieldvector,icompvector,
1769 nfieldvector1,&iselect,m2,f1,output,m3);
1770
1771 }
1772
1773 /* storing the heat generation in the nodes */
1774
1775 if((strcmp1(&filab[783],"RFL ")==0)&&(*ithermal>1)){
1776 iselect=1;
1777
1778 frdset(&filab[783],set,&iset,istartset,iendset,ialset,
1779 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1780 ngraph);
1781
1782 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1783 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1784
1785 fprintf(f1," -4 RFL 1 1\n");
1786 fprintf(f1," -5 RFL 1 1 0 0\n");
1787
1788 frdselect(fn,fn,&iset,&nkcoords,inum,m1,istartset,iendset,
1789 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1790 nfieldvector0,&iselect,m2,f1,output,m3);
1791
1792 }
1793
1794 /* storing the Zienkiewicz-Zhu improved stresses in the nodes */
1795
1796 if((*nmethod!=5)||(*mode==-1)){
1797 if((strcmp1(&filab[1044],"ZZS")==0)&&(*ithermal!=2)){
1798
1799 FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
1800 stx,&mi[0]));
1801
1802 iselect=1;
1803
1804 frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
1805 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1806 ngraph);
1807
1808 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1809 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1810
1811 fprintf(f1," -4 ZZSTR 6 1\n");
1812 fprintf(f1," -5 SXX 1 4 1 1\n");
1813 fprintf(f1," -5 SYY 1 4 2 2\n");
1814 fprintf(f1," -5 SZZ 1 4 3 3\n");
1815 fprintf(f1," -5 SXY 1 4 1 2\n");
1816 fprintf(f1," -5 SYZ 1 4 2 3\n");
1817 fprintf(f1," -5 SZX 1 4 3 1\n");
1818
1819 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1820 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1821 nfieldtensor,&iselect,m2,f1,output,m3);
1822
1823 }
1824 }
1825
1826 /* storing the imaginary part of the Zienkiewicz-Zhu
1827 improved stresses in the nodes
1828 for the odd modes of cyclic symmetry calculations */
1829
1830 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1831 if((strcmp1(&filab[1044],"ZZS")==0)&&(*ithermal!=2)){
1832
1833 FORTRAN(zienzhu,(co,nk,kon,ipkon,lakon,ne0,stn,ipneigh,neigh,
1834 &stx[6*mi[0]**ne],&mi[0]));
1835
1836 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1837 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1838
1839 fprintf(f1," -4 ZZSTRI 6 1\n");
1840 fprintf(f1," -5 SXX 1 4 1 1\n");
1841 fprintf(f1," -5 SYY 1 4 2 2\n");
1842 fprintf(f1," -5 SZZ 1 4 3 3\n");
1843 fprintf(f1," -5 SXY 1 4 1 2\n");
1844 fprintf(f1," -5 SYZ 1 4 2 3\n");
1845 fprintf(f1," -5 SZX 1 4 3 1\n");
1846
1847 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1848 ialset,ngraph,&ncomptensor,ifieldtensor,icomptensor,
1849 nfieldtensor,&iselect,m2,f1,output,m3);
1850
1851 }
1852 }
1853
1854 /* storing the error estimator in the nodes */
1855
1856 if((*nmethod!=5)||(*mode==-1)){
1857 if((strcmp1(&filab[1044],"ERR")==0)&&(*ithermal!=2)){
1858
1859 NNEW(errn,double,6**nk);
1860
1861 nterms=6;
1862 FORTRAN(errorestimator,(stx,errn,ipkon,kon,lakon,nk,ne,
1863 mi,ielmat,&nterms,inum,co,v,&filab[1048],
1864 ielprop,prop));
1865
1866 iselect=1;
1867
1868 frdset(&filab[1044],set,&iset,istartset,iendset,ialset,
1869 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1870 ngraph);
1871
1872 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1873 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1874
1875 fprintf(f1," -4 ERROR 1 1\n");
1876 fprintf(f1," -5 STR(%%) 1 1 0 0\n");
1877
1878 ncomp=1;
1879 ifield[0]=1;
1880 icomp[0]=0;
1881
1882 frdselect(errn,errn,&iset,&nkcoords,inum,m1,istartset,iendset,
1883 ialset,ngraph,&ncomp,ifield,icomp,
1884 nfieldtensor,&iselect,m2,f1,output,m3);
1885
1886 }
1887 }
1888
1889 /* storing the imaginary part of the error estimator in the nodes
1890 for the odd modes of cyclic symmetry calculations */
1891
1892 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1893 if((strcmp1(&filab[1044],"ERR")==0)&&(*ithermal!=2)){
1894
1895 nterms=6;
1896 FORTRAN(errorestimator,(&stx[6*mi[0]**ne],stn,ipkon,kon,lakon,nk,ne,
1897 mi,ielmat,&nterms,inum,co,v,&filab[1048],
1898 ielprop,prop));
1899
1900 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1901 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1902
1903 fprintf(f1," -4 ERRORI 1 1\n");
1904 fprintf(f1," -5 STR(%%) 1 1 0 0\n");
1905
1906 ncomp=1;
1907 ifield[0]=1;
1908 icomp[0]=0;
1909
1910 frdselect(stn,stn,&iset,&nkcoords,inum,m1,istartset,iendset,
1911 ialset,ngraph,&ncomp,ifield,icomp,
1912 nfieldtensor,&iselect,m2,f1,output,m3);
1913
1914 }
1915 }
1916
1917 /* storing the thermal error estimator in the nodes */
1918
1919 if((*nmethod!=5)||(*mode==-1)){
1920 if((strcmp1(&filab[2784],"HER")==0)&&(*ithermal>1)){
1921
1922 nterms=3;
1923 FORTRAN(errorestimator,(qfx,qfn,ipkon,kon,lakon,nk,ne,
1924 mi,ielmat,&nterms,inum,co,v,&filab[2788],
1925 ielprop,prop));
1926
1927 iselect=1;
1928
1929 frdset(&filab[2784],set,&iset,istartset,iendset,ialset,
1930 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1931 ngraph);
1932
1933 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1934 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1935
1936 fprintf(f1," -4 HERROR 1 1\n");
1937 fprintf(f1," -5 TEM(%%) 1 1 0 0\n");
1938
1939 ncomp=1;
1940 ifield[0]=1;
1941 icomp[0]=0;
1942
1943 frdselect(qfn,qfn,&iset,&nkcoords,inum,m1,istartset,iendset,
1944 ialset,ngraph,&ncomp,ifield,icomp,
1945 nfieldvector1,&iselect,m2,f1,output,m3);
1946
1947 }
1948 }
1949
1950 /* storing the imaginary part of the thermal error estimator in the nodes
1951 for the odd modes of cyclic symmetry calculations */
1952
1953 if((*noddiam>=0)||((*nmethod==5)&&(*mode==0))){
1954 if((strcmp1(&filab[2784],"HER")==0)&&(*ithermal>1)){
1955
1956 nterms=3;
1957 FORTRAN(errorestimator,(&qfx[3*mi[0]**ne],qfn,ipkon,kon,lakon,nk,ne,
1958 mi,ielmat,&nterms,inum,co,v,&filab[2788],
1959 ielprop,prop));
1960
1961 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1962 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1963
1964 fprintf(f1," -4 HERRORI 1 1\n");
1965 fprintf(f1," -5 TEM(%%) 1 1 0 0\n");
1966
1967 ncomp=1;
1968 ifield[0]=1;
1969 icomp[0]=0;
1970
1971 frdselect(qfn,qfn,&iset,&nkcoords,inum,m1,istartset,iendset,
1972 ialset,ngraph,&ncomp,ifield,icomp,
1973 nfieldtensor,&iselect,m2,f1,output,m3);
1974
1975 }
1976 }
1977
1978 /* storing the total temperatures in the network nodes */
1979
1980 if((strcmp1(&filab[1131],"TT ")==0)&&(*ithermal>1)){
1981
1982 iselect=-1;
1983 frdset(&filab[1131],set,&iset,istartset,iendset,ialset,
1984 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
1985 ngraph);
1986
1987 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
1988 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
1989
1990 fprintf(f1," -4 TOTEMP 1 1\n");
1991 fprintf(f1," -5 TT 1 1 0 0\n");
1992
1993 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
1994 ialset,ngraph,&ncompscalar,ifieldscalar,icompscalar,
1995 nfieldvector0,&iselect,m2,f1,output,m3);
1996
1997 }
1998
1999 /* storing the mass flow in the network nodes */
2000
2001 if((strcmp1(&filab[1218],"MF ")==0)&&(*ithermal>1)){
2002
2003 iselect=-1;
2004 frdset(&filab[1218],set,&iset,istartset,iendset,ialset,
2005 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2006 ngraph);
2007
2008 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2009 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2010
2011 fprintf(f1," -4 MAFLOW 1 1\n");
2012 fprintf(f1," -5 MF 1 1 0 0\n");
2013
2014 icomp[0]=1;
2015 if((iaxial==1)&&(strcmp1(&filab[1222],"I")==0)){for(i=0;i<*nk;i++)v[1+i*mt]*=180.;}
2016 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2017 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2018 nfieldvector0,&iselect,m2,f1,output,m3);
2019 if((iaxial==1)&&(strcmp1(&filab[1222],"I")==0)){for(i=0;i<*nk;i++)v[1+i*mt]/=180.;}
2020
2021 }
2022
2023 /* storing the total pressure in the network nodes */
2024
2025 if((strcmp1(&filab[1305],"PT ")==0)&&(*ithermal>1)){
2026
2027 iselect=-1;
2028 frdset(&filab[1305],set,&iset,istartset,iendset,ialset,
2029 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2030 ngraph);
2031
2032 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2033 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2034
2035 fprintf(f1," -4 TOPRES 1 1\n");
2036 fprintf(f1," -5 PT 1 1 0 0\n");
2037
2038 icomp[0]=2;
2039 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2040 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2041 nfieldvector0,&iselect,m2,f1,output,m3);
2042
2043 }
2044
2045 /* storing the static pressure in the liquid network nodes */
2046
2047 if((strcmp1(&filab[1827],"PS ")==0)&&(*ithermal>1)){
2048
2049 iselect=-1;
2050 frdset(&filab[1827],set,&iset,istartset,iendset,ialset,
2051 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2052 ngraph);
2053
2054 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2055 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2056
2057 fprintf(f1," -4 STPRES 1 1\n");
2058 fprintf(f1," -5 PS 1 1 0 0\n");
2059
2060 icomp[0]=2;
2061 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2062 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2063 nfieldvector0,&iselect,m2,f1,output,m3);
2064
2065 }
2066
2067 /* storing the liquid depth in the channel nodes */
2068
2069 if((strcmp1(&filab[2349],"DEPT")==0)&&(*ithermal>1)){
2070
2071 iselect=-1;
2072 frdset(&filab[2349],set,&iset,istartset,iendset,ialset,
2073 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2074 ngraph);
2075
2076 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2077 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2078
2079 fprintf(f1," -4 DEPTH 1 1\n");
2080 fprintf(f1," -5 DEPTH 1 1 0 0\n");
2081
2082 icomp[0]=2;
2083 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2084 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2085 nfieldvector0,&iselect,m2,f1,output,m3);
2086
2087 }
2088
2089 /* storing the critical depth in the channel nodes */
2090
2091 if((strcmp1(&filab[2436],"HCRI")==0)&&(*ithermal>1)){
2092
2093 iselect=-1;
2094 frdset(&filab[2436],set,&iset,istartset,iendset,ialset,
2095 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2096 ngraph);
2097
2098 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2099 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2100
2101 fprintf(f1," -4 HCRIT 1 1\n");
2102 fprintf(f1," -5 HCRIT 1 1 0 0\n");
2103
2104 icomp[0]=3;
2105 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2106 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2107 nfieldvector0,&iselect,m2,f1,output,m3);
2108
2109 }
2110
2111 /* storing the static temperature in the network nodes */
2112
2113 if((strcmp1(&filab[1392],"TS ")==0)&&(*ithermal>1)){
2114
2115 iselect=-1;
2116 frdset(&filab[1392],set,&iset,istartset,iendset,ialset,
2117 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2118 ngraph);
2119
2120 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2121 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2122
2123 fprintf(f1," -4 STTEMP 1 1\n");
2124 fprintf(f1," -5 TS 1 1 0 0\n");
2125
2126 icomp[0]=3;
2127 frdselect(v,v,&iset,&nkcoords,inum,m1,istartset,iendset,
2128 ialset,ngraph,&ncompscalar,ifieldscalar,icomp,
2129 nfieldvector0,&iselect,m2,f1,output,m3);
2130
2131 }
2132
2133 /* mesh refinement */
2134
2135 if(strcmp1(&filab[4089],"RM")==0){
2136 refinemesh(nk,ne,co,ipkon,kon,v,veold,stn,een,emn,epn,enern,
2137 qfn,errn,filab,mi,lakon,jobnamec,istartset,iendset,
2138 ialset,set,nset,matname,ithermal,output,nmat);
2139 }
2140
2141 /* remove auxiliary field for the error estimator at the nodes */
2142
2143 if((*nmethod!=5)||(*mode==-1)){
2144 if((strcmp1(&filab[1044],"ERR")==0)&&(*ithermal!=2)){
2145 SFREE(errn);
2146 }
2147 }
2148
2149 /* the remaining lines only apply to frequency calculations
2150 with cyclic symmetry, complex frequency and steady state calculations */
2151
2152 if((*nmethod!=2)&&(*nmethod!=13)&&(*nmethod!=5)&&(*nmethod!=6)&&(*nmethod!=7)){fclose(f1);return;}
2153 if((*nmethod==5)&&(*mode==-1)){fclose(f1);return;}
2154
2155 /* storing the displacements in the nodes (magnitude, phase) */
2156
2157 if((strcmp1(&filab[870],"PU ")==0)&&(*ithermal!=2)){
2158 iselect=1;
2159
2160 frdset(&filab[870],set,&iset,istartset,iendset,ialset,
2161 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2162 ngraph);
2163
2164 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2165 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2166
2167 fprintf(f1," -4 PDISP 6 1\n");
2168 fprintf(f1," -5 MAG1 1 12 1 0\n");
2169 fprintf(f1," -5 MAG2 1 12 2 0\n");
2170 fprintf(f1," -5 MAG3 1 12 3 0\n");
2171 fprintf(f1," -5 PHA1 1 12 4 0\n");
2172 fprintf(f1," -5 PHA2 1 12 5 0\n");
2173 fprintf(f1," -5 PHA3 1 12 6 0\n");
2174
2175 frdselect(vr,vi,&iset,&nkcoords,inum,m1,istartset,iendset,
2176 ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
2177 nfieldvectph,&iselect,m2,f1,output,m3);
2178
2179 }
2180
2181 /* storing the temperatures in the nodes (magnitude, phase) */
2182
2183 if((strcmp1(&filab[957],"PNT ")==0)&&(*ithermal>1)){
2184 iselect=1;
2185
2186 frdset(&filab[957],set,&iset,istartset,iendset,ialset,
2187 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2188 ngraph);
2189
2190 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2191 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2192
2193 fprintf(f1," -4 PNDTEMP 2 1\n");
2194 fprintf(f1," -5 MAG1 1 1 1 0\n");
2195 fprintf(f1," -5 PHA1 1 1 2 0\n");
2196
2197 frdselect(vr,vi,&iset,&nkcoords,inum,m1,istartset,iendset,
2198 ialset,ngraph,&ncompscalph,ifieldscalph,icompscalph,
2199 nfieldscalph,&iselect,m2,f1,output,m3);
2200
2201 }
2202
2203 /* storing the stresses in the nodes (magnitude, phase) */
2204
2205 if((strcmp1(&filab[1479],"PHS ")==0)&&(*ithermal!=2)){
2206 iselect=1;
2207
2208 frdset(&filab[1479],set,&iset,istartset,iendset,ialset,
2209 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2210 ngraph);
2211
2212 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2213 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2214
2215 fprintf(f1," -4 PSTRESS 12 1\n");
2216 fprintf(f1," -5 MAGXX 1 14 1 1\n");
2217 fprintf(f1," -5 MAGYY 1 14 2 2\n");
2218 fprintf(f1," -5 MAGZZ 1 14 3 3\n");
2219 fprintf(f1," -5 MAGXY 1 14 1 2\n");
2220 fprintf(f1," -5 MAGYZ 1 14 2 3\n");
2221 fprintf(f1," -5 MAGZX 1 14 3 1\n");
2222 fprintf(f1," -5 PHAXX 1 14 1 1\n");
2223 fprintf(f1," -5 PHAYY 1 14 2 2\n");
2224 fprintf(f1," -5 PHAZZ 1 14 3 3\n");
2225 fprintf(f1," -5 PHAXY 1 14 1 2\n");
2226 fprintf(f1," -5 PHAYZ 1 14 2 3\n");
2227 fprintf(f1," -5 PHAZX 1 14 3 1\n");
2228
2229 frdselect(stnr,stni,&iset,&nkcoords,inum,m1,istartset,iendset,
2230 ialset,ngraph,&ncomptensph,ifieldtensph,icomptensph,
2231 nfieldtensph,&iselect,m2,f1,output,m3);
2232
2233 }
2234
2235 /* storing the differential contact displacements and
2236 the contact stresses in the nodes (magnitude, phase)
2237 only for face-to-face penalty contact */
2238
2239 if((strcmp1(&filab[3915],"PCON")==0)&&(*ithermal!=2)&&(*mortar==1)){
2240 iselect=1;
2241
2242 frdset(&filab[3915],set,&iset,istartset,iendset,ialset,
2243 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2244 ngraph);
2245
2246 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2247 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2248
2249 fprintf(f1," -4 PCONTAC 12 1\n");
2250 fprintf(f1," -5 MAGO 1 14 1 1\n");
2251 fprintf(f1," -5 MAGSL1 1 14 2 2\n");
2252 fprintf(f1," -5 MAGSL2 1 14 3 3\n");
2253 fprintf(f1," -5 MAGP 1 14 1 2\n");
2254 fprintf(f1," -5 MAGSH1 1 14 2 3\n");
2255 fprintf(f1," -5 MAGSH2 1 14 3 1\n");
2256 fprintf(f1," -5 PHAO 1 14 1 1\n");
2257 fprintf(f1," -5 PHASL1 1 14 2 2\n");
2258 fprintf(f1," -5 PHASL2 1 14 3 3\n");
2259 fprintf(f1," -5 PHAP 1 14 1 2\n");
2260 fprintf(f1," -5 PHASH1 1 14 2 3\n");
2261 fprintf(f1," -5 PHASH2 1 14 3 1\n");
2262
2263 frdselect(cdnr,cdni,&iset,&nkcoords,inum,m1,istartset,iendset,
2264 ialset,ngraph,&ncomptensph,ifieldtensph,icomptensph,
2265 nfieldtensph,&iselect,m2,f1,output,m3);
2266
2267 }
2268
2269 /* storing the forces in the nodes (magnitude, phase) */
2270
2271 if((strcmp1(&filab[2610],"PRF ")==0)&&(*ithermal!=2)){
2272 iselect=1;
2273
2274 frdset(&filab[2610],set,&iset,istartset,iendset,ialset,
2275 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2276 ngraph);
2277
2278 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2279 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2280
2281 fprintf(f1," -4 PFORC 6 1\n");
2282 fprintf(f1," -5 MAG1 1 12 1 0\n");
2283 fprintf(f1," -5 MAG2 1 12 2 0\n");
2284 fprintf(f1," -5 MAG3 1 12 3 0\n");
2285 fprintf(f1," -5 PHA1 1 12 4 0\n");
2286 fprintf(f1," -5 PHA2 1 12 5 0\n");
2287 fprintf(f1," -5 PHA3 1 12 6 0\n");
2288
2289 frdselect(fnr,fni,&iset,&nkcoords,inum,m1,istartset,iendset,
2290 ialset,ngraph,&ncompvectph,ifieldvectph,icompvectph,
2291 nfieldvectph,&iselect,m2,f1,output,m3);
2292
2293 }
2294
2295 /* the remaining parts are for frequency calculations with cyclic symmetry only */
2296
2297 if((*nmethod!=2)&&(*nmethod!=13)){fclose(f1);return;}
2298
2299 /* storing the maximum displacements of the nodes in the base sector
2300 (components, magnitude) */
2301
2302 if((strcmp1(&filab[1566],"MAXU")==0)&&(*ithermal!=2)){
2303 iselect=1;
2304
2305 frdset(&filab[1566],set,&iset,istartset,iendset,ialset,
2306 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2307 ngraph);
2308
2309 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2310 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2311
2312 fprintf(f1," -4 MDISP 4 1\n");
2313 fprintf(f1," -5 DX 1 4 1 0\n");
2314 fprintf(f1," -5 DY 1 4 2 0\n");
2315 fprintf(f1," -5 DZ 1 4 3 0\n");
2316 fprintf(f1," -5 ANG 1 4 4 0\n");
2317
2318 ncomp=4;
2319 ifield[0]=1;icomp[0]=1;
2320 ifield[1]=1;icomp[1]=2;
2321 ifield[2]=1;icomp[2]=3;
2322 ifield[3]=1;icomp[3]=0;
2323 nfield[0]=4;nfield[1]=4;
2324
2325 frdselect(vmax,vmax,&iset,&nkcoords,inum,m1,istartset,iendset,
2326 ialset,ngraph,&ncomp,ifield,icomp,
2327 nfield,&iselect,m2,f1,output,m3);
2328
2329 }
2330
2331 /* storing the worst principal stress at the nodes
2332 in the basis sector (components, magnitude)
2333
2334 the worst principal stress is the maximum of the
2335 absolute value of all principal stresses, times
2336 its original sign */
2337
2338 if((strcmp1(&filab[1653],"MAXS")==0)&&(*ithermal!=2)){
2339 iselect=1;
2340
2341 frdset(&filab[1653],set,&iset,istartset,iendset,ialset,
2342 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2343 ngraph);
2344
2345 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2346 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2347
2348 fprintf(f1," -4 MSTRESS 7 1\n");
2349 fprintf(f1," -5 SXX 1 4 1 1\n");
2350 fprintf(f1," -5 SYY 1 4 2 2\n");
2351 fprintf(f1," -5 SZZ 1 4 3 3\n");
2352 fprintf(f1," -5 SXY 1 4 1 2\n");
2353 fprintf(f1," -5 SYZ 1 4 2 3\n");
2354 fprintf(f1," -5 SZX 1 4 3 1\n");
2355 fprintf(f1," -5 MAG 1 4 0 0\n");
2356
2357 ncomp=7;
2358 ifield[0]=1;icomp[0]=1;
2359 ifield[1]=1;icomp[1]=2;
2360 ifield[2]=1;icomp[2]=3;
2361 ifield[3]=1;icomp[3]=4;
2362 ifield[4]=1;icomp[4]=6;
2363 ifield[5]=1;icomp[5]=5;
2364 ifield[6]=1;icomp[6]=0;
2365 nfield[0]=7;nfield[1]=7;
2366
2367 frdselect(stnmax,stnmax,&iset,&nkcoords,inum,m1,istartset,iendset,
2368 ialset,ngraph,&ncomp,ifield,icomp,
2369 nfield,&iselect,m2,f1,output,m3);
2370
2371 }
2372
2373 /* storing the worst principal strain at the nodes
2374 in the basis sector (components, magnitude)
2375
2376 the worst principal strain is the maximum of the
2377 absolute value of all principal strains, times
2378 its original sign */
2379
2380 if((strcmp1(&filab[2523],"MAXE")==0)&&(*ithermal!=2)){
2381 iselect=1;
2382
2383 frdset(&filab[2523],set,&iset,istartset,iendset,ialset,
2384 inum,&noutloc,&nout,nset,&noutmin,&noutplus,&iselect,
2385 ngraph);
2386
2387 frdheader(&icounter,&oner,time,&pi,noddiam,cs,&null,mode,
2388 &noutloc,description,kode,nmethod,f1,output,istep,iinc);
2389
2390 fprintf(f1," -4 MSTRAIN 7 1\n");
2391 fprintf(f1," -5 EXX 1 4 1 1\n");
2392 fprintf(f1," -5 EYY 1 4 2 2\n");
2393 fprintf(f1," -5 EZZ 1 4 3 3\n");
2394 fprintf(f1," -5 EXY 1 4 1 2\n");
2395 fprintf(f1," -5 EYZ 1 4 2 3\n");
2396 fprintf(f1," -5 EZX 1 4 3 1\n");
2397 fprintf(f1," -5 MAG 1 4 0 0\n");
2398
2399 ncomp=7;
2400 ifield[0]=1;icomp[0]=1;
2401 ifield[1]=1;icomp[1]=2;
2402 ifield[2]=1;icomp[2]=3;
2403 ifield[3]=1;icomp[3]=4;
2404 ifield[4]=1;icomp[4]=6;
2405 ifield[5]=1;icomp[5]=5;
2406 ifield[6]=1;icomp[6]=0;
2407 nfield[0]=7;nfield[1]=7;
2408
2409 frdselect(eenmax,eenmax,&iset,&nkcoords,inum,m1,istartset,iendset,
2410 ialset,ngraph,&ncomp,ifield,icomp,
2411 nfield,&iselect,m2,f1,output,m3);
2412
2413 }
2414
2415 fclose(f1);
2416 return;
2417
2418 }
2419