1 /***********************************************************************
2 *
3 * ***** *** ***
4 * * * * * *
5 * * *** ***
6 * * * * * *
7 * ***** *** ***
8 *
9 * A FREE Finite Elements Analysis Program in ANSI C for the Windows & UNIX OS.
10 *
11 * Composed and edited and copyright by
12 * Professor Dr.-Ing. Frank Rieg, University of Bayreuth, Germany
13 *
14 * eMail:
15 * frank.rieg@uni-bayreuth.de
16 * dr.frank.rieg@t-online.de
17 *
18 * V15.0 November 18, 2015
19 *
20 * Z88 should compile and run under any UNIX OS and Windows.
21 *
22 * This program is free software; you can redistribute it and/or modify
23 * it under the terms of the GNU General Public License as published by
24 * the Free Software Foundation; either version 2, or (at your option)
25 * any later version.
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU General Public License for more details.
31 *
32 * You should have received a copy of the GNU General Public License
33 * along with this program; see the file COPYING. If not, write to
34 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
35 ***********************************************************************/
36 /***********************************************************************
37 * Compiler-Unit m4.c enthaelt:
38 * sshe88.for : 6-knoten scheibe
39 * srin88.for : 3-knoten torus
40 * 2.1.2010 Rieg
41 ***********************************************************************/
42
43 /***********************************************************************
44 * Fuer UNIX
45 ***********************************************************************/
46 #ifdef FR_UNIX
47 #include <z88r.h>
48 #include <stdio.h> /* FILE,fprintf */
49 #endif
50
51 /***********************************************************************
52 * Fuer Windows
53 ***********************************************************************/
54 #ifdef FR_WIN
55 #include <z88r.h>
56 #include <stdio.h> /* FILE,fprintf */
57 #endif
58
59 /***********************************************************************
60 * Fuer Windows & GTK+
61 ***********************************************************************/
62 #ifdef FR_GTKWIN
63 #include <z88r.h>
64 #include <stdio.h>
65 #endif
66
67 /***********************************************************************
68 * Formate
69 ***********************************************************************/
70 #define NL "\n"
71
72 #ifdef FR_XDOUB
73 #define P10E " %+#10.2lE"
74 #define P11E " %+#11.3lE"
75 #endif
76
77 #ifdef FR_XQUAD
78 #define P10E " %+#10.2LE"
79 #define P11E " %+#11.3LE"
80 #endif
81
82 /***********************************************************************
83 * Functions
84 ***********************************************************************/
85 FR_DOUBLE sheigh(FR_DOUBLE sig[]);
86 FR_DOUBLE torgh(FR_DOUBLE sig[]);
87 FR_DOUBLE sheinh(FR_DOUBLE sig[]);
88 FR_DOUBLE tornh(FR_DOUBLE sig[]);
89 FR_DOUBLE sheish(FR_DOUBLE sig[]);
90 FR_DOUBLE torsh(FR_DOUBLE sig[]);
91
92 /***********************************************************************
93 * hier beginnt Function sshe88
94 ***********************************************************************/
sshe88(void)95 int sshe88(void)
96 {
97 extern FILE *fo3,*fo5;
98
99 extern FR_DOUBLEAY gmw;
100 extern FR_DOUBLEAY sigvku;
101
102 extern FR_DOUBLE ul[];
103 extern FR_DOUBLE xk[],yk[];
104
105 extern FR_DOUBLE emode,rnuee;
106
107 extern FR_INT4 kflag,isflag,kc,jpri,ifnili,ngau;
108
109 FR_DOUBLE sig[4];
110
111 FR_DOUBLE rr,phirad,phideg,sinphi,cosphi;
112 FR_DOUBLE qsin,qcos,sincos;
113 FR_DOUBLE x21,x31,y21,y31,det,uxsi,vxsi,ueta,veta,xs,ys;
114 FR_DOUBLE epsx,epsy,gamxy,sigr,sigt,taurt,sigv;
115
116 /***********************************************************************
117 * Hilfswerte berechnen
118 ***********************************************************************/
119 x21= xk[2] - xk[1];
120 x31= xk[3] - xk[1];
121 y21= yk[2] - yk[1];
122 y31= yk[3] - yk[1];
123
124 /***********************************************************************
125 * Jacobi-Determinate
126 ***********************************************************************/
127 det= x21 * y31 - x31 * y21;
128
129 /***********************************************************************
130 * Partielle Ableitungen im Einheitsdreieck glch.2.169,s.123
131 ***********************************************************************/
132 uxsi= (ul[3] - ul[1] + 4. * (ul[9] - ul[11])) / 3.;
133 vxsi= (ul[4] - ul[2] + 4. * (ul[10] - ul[12])) / 3.;
134 ueta= (ul[5] - ul[1] + 4. * (ul[9] - ul[7] )) / 3.;
135 veta= (ul[6] - ul[2] + 4. * (ul[10] - ul[8] )) / 3.;
136
137 /***********************************************************************
138 * Dehnungen lt. Glch. 2.167, s.123
139 ***********************************************************************/
140 epsx= (y31 * uxsi - y21 * ueta) / det;
141 epsy= (x21 * veta - x31 * vxsi) / det;
142 gamxy= (x21 * ueta - x31 * uxsi + y31 * vxsi - y21 * veta) / det;
143
144 /***********************************************************************
145 * Spannungen lt. Glch. 2.166, s.123
146 ***********************************************************************/
147 sig[1]= emode * (epsx + rnuee * epsy) / (1. - rnuee * rnuee);
148 sig[2]= emode * (rnuee * epsx + epsy) / (1. - rnuee * rnuee);
149 sig[3]= emode * gamxy / (2. * (1. + rnuee));
150
151 /***********************************************************************
152 * Entscheiden, ob mit oder ohne Polarkoordinaten
153 ***********************************************************************/
154 xs= (xk[1] + xk[2] + xk[3]) / 3.;
155 ys= (yk[1] + yk[2] + yk[3]) / 3.;
156
157 if(isflag == 0) /* ohne VglS */
158 {
159 if(kflag == 1) /* mit sigr */
160 {
161 rr= FR_SQRT(xs * xs + ys * ys);
162 if(xs == 0.) xs= 1e-10;
163 phirad= FR_ATAN(ys/xs);
164 phideg= phirad*57.29578;
165
166 sinphi= FR_SIN(phirad);
167 cosphi= FR_COS(phirad);
168 qsin= sinphi*sinphi;
169 qcos= cosphi*cosphi;
170 sincos= sinphi*cosphi;
171
172 sigr= sig[1] * qcos + sig[2] * qsin + 2. * sig[3] * sincos;
173 sigt= sig[1] * qsin + sig[2] * qcos - 2. * sig[3] * sincos;
174 taurt= (sig[2] - sig[1]) * sincos + sig[3] * (qcos-qsin);
175
176 if(ifnili == 0)
177 {
178 fprintf(fo3,NL P10E P10E P11E P11E P11E P10E P10E P11E P11E P11E
179 ,xs,ys,sig[1],sig[2],sig[3],rr,phideg,sigr,sigt,taurt);
180 }
181 } /* e if mit sigr */
182 else
183 { /* ohne sigr */
184 if(ifnili == 0)
185 {
186 fprintf(fo3,NL P10E P10E P11E P11E P11E
187 ,xs,ys,sig[1],sig[2],sig[3]);
188 } /* e if jpri */
189 } /* e if ohne sigr */
190 } /* e if ohne GEH */
191
192 if(isflag == 1 || isflag == 2 || isflag == 3) /* mit GEH */
193 {
194 if(isflag == 1) sigv= sheigh(sig);
195 if(isflag == 2) sigv= sheinh(sig);
196 if(isflag == 3) sigv= sheish(sig);
197
198 ngau++;
199 sigvku[ngau]+= sigv;
200
201 if(kflag == 1) /* mit sigr */
202 {
203 rr= FR_SQRT(xs * xs + ys * ys);
204 if(xs == 0.) xs= 1e-10;
205 phirad= FR_ATAN(ys/xs);
206 phideg= phirad*57.29578;
207
208 sinphi= FR_SIN(phirad);
209 cosphi= FR_COS(phirad);
210 qsin= sinphi*sinphi;
211 qcos= cosphi*cosphi;
212 sincos= sinphi*cosphi;
213
214 sigr= sig[1] * qcos + sig[2] * qsin + 2. * sig[3] * sincos;
215 sigt= sig[1] * qsin + sig[2] * qcos - 2. * sig[3] * sincos;
216 taurt= (sig[2] - sig[1]) * sincos + sig[3] * (qcos-qsin);
217
218 if(ifnili == 0)
219 {
220 fprintf(fo3,NL P10E P10E P11E P11E P11E P10E P10E P11E P11E P11E P11E
221 ,xs,ys,sig[1],sig[2],sig[3],rr,phideg,sigr,sigt,taurt,sigv);
222 }
223 } /* e if mit sigr */
224 else
225 { /* ohne sigr */
226 if(ifnili == 0)
227 {
228 fprintf(fo3,P10E P10E P11E P11E P11E P11E
229 ,xs,ys,sig[1],sig[2],sig[3],sigv);
230 }
231 } /* e if ohne sigr */
232
233 if(jpri == 1)
234 {
235 fprintf(fo5,NL P11E P11E P11E,xs,ys,sigvku[ngau]);
236 }
237 gmw[kc]= sigv;
238 } /* e if mit GEH */
239
240 return(0);
241 }
242
243 /***********************************************************************
244 * hier beginnt Function srin88
245 ***********************************************************************/
srin88(void)246 int srin88(void)
247 {
248 extern FILE *fo3,*fo5;
249
250 extern FR_DOUBLEAY gmw;
251 extern FR_DOUBLEAY sigvku;
252
253 extern FR_DOUBLE ul[];
254 extern FR_DOUBLE xk[],yk[];
255
256 extern FR_DOUBLE emode,rnuee;
257
258 extern FR_INT4 isflag,kc,jpri,ifnili,ngau;
259
260 FR_DOUBLE a[4],b[4],c[4],d[4];
261 FR_DOUBLE hdb[25]; /* 4 x 6 +1 */
262 FR_DOUBLE sig[5],sigmit[5];
263
264 FR_DOUBLE f0,f1,f2,flaec,tsig1,tsig2,tsig3,tsig4,sigv,xx,yy,xs,ys;
265
266 FR_INT4 k,i,j;
267
268 static FR_INT4 izyk[6]= { 0,
269 1,2,3,1,2 };
270
271 /***********************************************************************
272 * Flaeche des Dreiecks
273 ***********************************************************************/
274 flaec= 0.5 * ( xk[1] * (yk[2] - yk[3]) + xk[2] * (yk[3] - yk[1]) +
275 + xk[3] * (yk[1] - yk[2]) );
276
277 /***********************************************************************
278 * Berechnung f0, f1 & f2
279 ***********************************************************************/
280 f0= emode*(1.-rnuee)/(1.+rnuee)/(1.-2.*rnuee);
281 f1= rnuee/(1.-rnuee);
282 f2= (1.-2.*rnuee)/2./(1.-rnuee);
283
284 /***********************************************************************
285 * Berechnung von a(i), b(i) ,c(i) & d(i)
286 ***********************************************************************/
287 for(i = 1;i <= 4;i++)
288 sigmit[i]= 0.;
289
290 for(i = 1;i <= 3;i++) /* ueber 3 Knoten */
291 { /* 70 */
292 xx= xk[i];
293 yy= yk[i];
294 for(k = 1;k <= 3;k++)
295 {
296 a[k]= xk[izyk[k+1]] * yk[izyk[k+2]] - xk[izyk[k+2]] * yk[izyk[k+1]];
297 b[k]= yk[izyk[k+1]] - yk[izyk[k+2]];
298 c[k]= xk[izyk[k+2]] - xk[izyk[k+1]];
299 d[k]= a[k] / xx + b[k] + c[k] * yy / xx;
300 } /* e 30 */
301
302 /***********************************************************************
303 * Matrix hdb = d*b berechnen
304 ***********************************************************************/
305 hdb[1 ]= f1 * b[1] + f1 * d[1];
306 hdb[7 ]= b[1] + f1 * d[1];
307 hdb[13]= f1 * b[1] + d[1];
308 hdb[19]= f2 * c[1];
309
310 hdb[2 ]= c[1];
311 hdb[8 ]= f1 * c[1];
312 hdb[14]= f1 * c[1];
313 hdb[20]= f2 * b[1];
314
315 hdb[3 ]= f1 * b[2] + f1 * d[2];
316 hdb[9 ]= b[2] + f1 * d[2];
317 hdb[15]= f1 * b[2] + d[2];
318 hdb[21]= f2 * c[2];
319
320 hdb[4 ]= c[2];
321 hdb[10]= f1 * c[2];
322 hdb[16]= f1 * c[2];
323 hdb[22]= f2 * b[2];
324
325 hdb[5 ]= f1 * b[3] + f1 * d[3];
326 hdb[11]= b[3] + f1 * d[3];
327 hdb[17]= f1 * b[3] + d[3];
328 hdb[23]= f2 * c[3];
329
330 hdb[6 ]= c[3];
331 hdb[12]= f1 * c[3];
332 hdb[18]= f1 * c[3];
333 hdb[24]= f2 * b[3];
334
335 /***********************************************************************
336 * nun Spannungen berechnen
337 ***********************************************************************/
338 for(j = 1;j <= 4;j++)
339 {
340 sig[j]= 0.;
341 for(k = 1;k <= 6;k++)
342 {
343 sig[j]= sig[j] + hdb[(j-1)*6 + k] * ul[k] * f0 / 2. / flaec;
344 }
345 sigmit[j]= sigmit[j] + sig[j];
346 }
347 } /* e70 */
348
349 /*---------------------------------------------------------------------
350 * mittlere Spannungen im Elementschwerpunkt berechnen
351 *--------------------------------------------------------------------*/
352 for(i = 1;i <= 4;i++)
353 sigmit[i]= sigmit[i] * 0.33333333; /* durch 3, weil 3 Knoten */
354
355 /*---------------------------------------------------------------------
356 * bis jetzt gilt:
357 * sigmit[1] = ueber die Knoten gemittelte SIGZZ
358 * sigmit[2] = ueber die Knoten gemittelte SIGRR
359 * sigmit[3] = ueber die Knoten gemittelte SIGTE
360 * sigmit[4] = ueber die Knoten gemittelte TAURZ
361 * ab jetzt werden die sigmit vertauscht: fuer torgh und Kompatibilitaet
362 * zu sqsh88 und scsh88
363 *--------------------------------------------------------------------*/
364 tsig1= sigmit[1];
365 tsig2= sigmit[2];
366 tsig3= sigmit[3];
367 tsig4= sigmit[4];
368 sigmit[1]= tsig2;
369 sigmit[2]= tsig1;
370 sigmit[3]= tsig4;
371 sigmit[4]= tsig3;
372
373 /*---------------------------------------------------------------------
374 * ab jetzt gilt:
375 * sigmit[1] = ueber die Knoten gemittelte SIGRR
376 * sigmit[2] = ueber die Knoten gemittelte SIGZZ
377 * sigmit[3] = ueber die Knoten gemittelte TAURZ
378 * sigmit[4] = ueber die Knoten gemittelte SIGTE
379 *--------------------------------------------------------------------*/
380
381 /***********************************************************************
382 * Beschreiben der Files Z88O3.TXT & Z88O5.TXT
383 ***********************************************************************/
384 xs= 0.33333333 * (xk[1] + xk[2] + xk[3]);
385 ys= 0.33333333 * (yk[1] + yk[2] + yk[3]);
386
387 if(isflag == 0 && ifnili == 0) /* ohne VglS */
388 {
389 fprintf(fo3,NL P10E P10E P11E P11E P11E P11E
390 ,xs,ys,sigmit[1],sigmit[2],sigmit[3],sigmit[4]);
391 }
392
393 if(isflag == 1 || isflag == 2 || isflag == 3) /* mit GEH */
394 {
395 if(isflag == 1) sigv= torgh(sigmit);
396 if(isflag == 2) sigv= tornh(sigmit);
397 if(isflag == 3) sigv= torsh(sigmit);
398
399 ngau++;
400 sigvku[ngau]+= sigv;
401
402 if(ifnili == 0)
403 {
404 fprintf(fo3,NL P10E P10E P11E P11E P11E P11E P11E
405 ,xs,ys,sigmit[1],sigmit[2],sigmit[3],sigmit[4],sigv);
406 }
407 if(jpri == 1)
408 {
409 fprintf(fo5,NL P11E P11E P11E,xs,ys,sigvku[ngau]);
410 }
411 gmw[kc]= sigv;
412 } /* e if mit GEH */
413 return(0);
414 }
415