1 /************************************************************************
2 *
3 * ***** *** ***
4 * * * * * *
5 * * *** ***
6 * * * * * *
7 * ***** *** ***
8 *
9 * A FREE Finite Elements Analysis Program in ANSI C for the Windows
10 * and UNIX OS.
11 *
12 * Composed and edited and copyright by
13 * Professor Dr.-Ing. Frank Rieg, University of Bayreuth, Germany
14 *
15 * eMail:
16 * frank.rieg@uni-bayreuth.de
17 * dr.frank.rieg@t-online.de
18 *
19 * V15.0 November 18, 2015
20 *
21 * This program is free software; you can redistribute it and/or modify
22 * it under the terms of the GNU General Public License as published by
23 * the Free Software Foundation; either version 2, or (at your option)
24 * any later version.
25 *
26 * This program is distributed in the hope that it will be useful,
27 * but WITHOUT ANY WARRANTY; without even the implied warranty of
28 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
29 * GNU General Public License for more details.
30 *
31 * You should have received a copy of the GNU General Public License
32 * along with this program; see the file COPYING. If not, write to
33 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
34 ***********************************************************************/
35 /***********************************************************************
36 * diese Compilerunit umfasst: shad88 - Elementsteifigkeitsroutine
37 * dd88 - Berechnung der Matrix b
38 * diese Compilerunit enthaelt Routinen, die gedanklich an FORTRAN-
39 * Quellen von H.J.Bathe, MIT, Cambridge, MA, USA angelehnt sind.
40 * 17.7.2011 Rieg
41 ***********************************************************************/
42
43 /***********************************************************************
44 * Fuer UNIX
45 ***********************************************************************/
46 #ifdef FR_UNIX
47 #include <z88r.h>
48 #include <stdio.h>
49 #endif
50
51 /***********************************************************************
52 * Fuer Windows
53 ***********************************************************************/
54 #ifdef FR_WIN
55 #include <z88r.h>
56 #include <stdio.h>
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 * Functions
69 ***********************************************************************/
70 int dd88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s,FR_DOUBLE *t);
71
72 /***********************************************************************
73 * hier beginnt Function shad88
74 ***********************************************************************/
shad88(void)75 int shad88(void)
76 {
77 extern FR_DOUBLEAY se;
78
79 extern FR_DOUBLE xk[],yk[],zk[];
80 extern FR_DOUBLE b[],xx[],d[];
81
82 extern FR_DOUBLE emode,rnuee;
83
84 extern FR_INT4 intore;
85
86 FR_DOUBLE db[7];
87
88 FR_DOUBLE r,s,t,det,wt,stiff;
89 FR_DOUBLE f0,f1,f2;
90
91 FR_INT4 ne= 36,i,igauss,lz,j,k,l,intor_t;
92
93 int iret;
94
95 /*----------------------------------------------------------------------
96 * Gauss-Legendre Stuetzstellen fuer r
97 *---------------------------------------------------------------------*/
98 static FR_DOUBLE rg[40]=
99 {
100 0.,0.,0.,0.,0.,0.,0., /* Elemente 0 - 6 leer */
101 0.1666666666667, /* intore = 3, 1.Ele Start bei i=7 */
102 0.6666666666667,
103 0.1666666666667,
104 0.,0.,0.,0.,0., /* Elemente 10-14 leer */
105 0.1012865073235, /* intore = 7, 1.Ele Start bei i=15 */
106 0.7974269853531,
107 0.1012865073235,
108 0.4701420641051,
109 0.4701420641051,
110 0.0597158717898,
111 0.3333333333333,
112 0.,0.,0.,0.,0., /* Elemente 22-26 leer */
113 0.0651301029022, /* intore =13, 1.Ele Start bei i=27 */
114 0.8697397941956,
115 0.0651301029022,
116 0.3128654960049,
117 0.6384441885698,
118 0.0486903154253,
119 0.6384441885698,
120 0.3128654960049,
121 0.0486903154253,
122 0.2603459660790,
123 0.4793080678419,
124 0.2603459660790,
125 0.3333333333333
126 };
127
128 /*----------------------------------------------------------------------
129 * Gauss-Legendre Stuetzstellen fuer s
130 *---------------------------------------------------------------------*/
131 static FR_DOUBLE sg[40]=
132 {
133 0.,0.,0.,0.,0.,0.,0., /* Elemente 0 - 6 leer */
134 0.1666666666667, /* intore = 3, 1.Ele Start bei i=7 */
135 0.1666666666667,
136 0.6666666666667,
137 0.,0.,0.,0.,0., /* Elemente 10-14 leer */
138 0.1012865073235, /* intore = 7, 1.Ele Start bei i=15 */
139 0.1012865073235,
140 0.7974269853531,
141 0.0597158717898,
142 0.4701420641051,
143 0.4701420641051,
144 0.3333333333333,
145 0.,0.,0.,0.,0., /* Elemente 22-26 leer */
146 0.0651301029022, /* intore =13, 1.Ele Start bei i=27 */
147 0.0651301029022,
148 0.8697397941956,
149 0.0486903154253,
150 0.3128654960049,
151 0.6384441885698,
152 0.0486903154253,
153 0.6384441885698,
154 0.3128654960049,
155 0.2603459660790,
156 0.2603459660790,
157 0.4793080678419,
158 0.3333333333333
159 };
160
161 /*----------------------------------------------------------------------
162 * Gauss-Legendre Integrationsgewichte
163 *---------------------------------------------------------------------*/
164 static FR_DOUBLE wg[40]=
165 {
166 0.,0.,0.,0.,0.,0.,0., /* Elemente 0 - 6 leer */
167 0.3333333333333, /* intore = 3, 1.Ele Start bei i=7 */
168 0.3333333333333,
169 0.3333333333333,
170 0.,0.,0.,0.,0., /* Elemente 10-14 leer */
171 0.1259391805448, /* intore = 7, 1.Ele Start bei i=15 */
172 0.1259391805448,
173 0.1259391805448,
174 0.1323941527885,
175 0.1323941527885,
176 0.1323941527885,
177 0.225,
178 0.,0.,0.,0.,0., /* Elemente 22-26 leer */
179 0.0533472356088, /* intore =13, 1.Ele Start bei i=27 */
180 0.0533472356088,
181 0.0533472356088,
182 0.0771137608903,
183 0.0771137608903,
184 0.0771137608903,
185 0.0771137608903,
186 0.0771137608903,
187 0.0771137608903,
188 0.1756152574332,
189 0.1756152574332,
190 0.1756152574332,
191 -0.1495700444677
192 };
193
194 /*----------------------------------------------------------------------
195 * Gauss-Legendre Stuetzstellen
196 *---------------------------------------------------------------------*/
197 static FR_DOUBLE xg[17]= { 0.,
198 0., -.5773502691896, -.7745966692415, -.8611363115941,
199 0., +.5773502691896, 0., -.3399810435849,
200 0., 0., +.7745966692415, +.3399810435849,
201 0., 0., 0., +.8611363115941 };
202
203 /*----------------------------------------------------------------------
204 * Gauss-Legendre Integrationsgewichte
205 *---------------------------------------------------------------------*/
206 static FR_DOUBLE wgt[17]= { 0.,
207 2., 1., +.5555555555556, +.3478548451375,
208 0., 1., +.8888888888889, +.6521451548625,
209 0., 0., +.5555555555556, +.6521451548625,
210 0., 0., 0., +.3478548451375 };
211
212 /*----------------------------------------------------------------------
213 * xk, yk und zk umspeichern
214 *---------------------------------------------------------------------*/
215 for(i = 1;i <= 12;i++)
216 {
217 xx[i] = xk[i];
218 xx[12+i]= yk[i];
219 xx[24+i]= zk[i];
220 }
221
222 /*----------------------------------------------------------------------
223 * Materialkonstanten
224 *---------------------------------------------------------------------*/
225 f0= emode*(1.-rnuee) / ((1.+rnuee)*(1.-2.*rnuee));
226
227 f1= rnuee/(1.-rnuee) * f0;
228 f2= (1.-2.*rnuee) / (2.*(1.-rnuee)) * f0;
229
230 /*----------------------------------------------------------------------
231 * Elastizitaetsmatrix aufstellen
232 *---------------------------------------------------------------------*/
233 for(i = 1;i <= 36;i++)
234 d[i]= 0.;
235
236 d[1] = f0;
237 d[7] = f1;
238 d[13]= f1;
239 d[2] = f1;
240 d[8] = f0;
241 d[14]= f1;
242 d[3] = f1;
243 d[9] = f1;
244 d[15]= f0;
245 d[22]= f2;
246 d[29]= f2;
247 d[36]= f2;
248
249 /*----------------------------------------------------------------------
250 * Elementsteifigkeitsmatrix aufstellen
251 *---------------------------------------------------------------------*/
252 intor_t= 2;
253
254 for(i = 1;i <= 1296;i++) /* 36 x 36 */
255 se[i]= 0.;
256
257 for(igauss = 1;igauss <= intore;igauss++)
258 {
259 r= rg[igauss+2*intore];
260 s= sg[igauss+2*intore];
261 for(lz = 1;lz <= intor_t;lz++)
262 {
263 t= xg[(lz-1)*4 + intor_t];
264
265 /*======================================================================
266 * Matrix b der partiellen Ableitungen & Jacobi Determinante holen
267 *=====================================================================*/
268 iret= dd88(&det,&r,&s,&t);
269 if(iret != 0) return(iret);
270
271 wt= wg[igauss+2*intore]*0.5 * wgt[(lz-1)*4+intor_t] * det;
272
273 for(j = 1;j <= 36;j++)
274 {
275 for(k = 1;k <= 6;k++)
276 {
277 db[k]= 0.;
278 for(l = 1;l <= 6;l++)
279 {
280 db[k]= db[k] + d[(k-1)*6 + l] * b[(l-1)*36 + j];
281 }
282 }
283
284 for(i = j;i <= 36;i++)
285 {
286 stiff= 0.;
287 for(l = 1;l <= 6;l++)
288 {
289 stiff+= b[(l-1)*36 + i] * db[l];
290 }
291 se[i+ne*(j-1)]= se[i+ne*(j-1)] + stiff * wt;
292 }
293 }
294 }
295 }
296
297 for(j = 1;j <= 36;j++)
298 {
299 for(i = j;i <= 36;i++)
300 {
301 se[j+ne*(i-1)]= se[i+ne*(j-1)];
302 }
303 }
304
305 return(0);
306 }
307
308 /***********************************************************************
309 * hier beginnt Function dd88
310 ***********************************************************************/
dd88(FR_DOUBLE * det,FR_DOUBLE * r,FR_DOUBLE * s,FR_DOUBLE * t)311 int dd88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s,FR_DOUBLE *t)
312 {
313 /*---------------------------------------------------------------------
314 * xx geht rein, unveraendert (ex)
315 * b geht raus, neu (ex)
316 * det geht raus, neu
317 * r,s,t gehen rein, unveraendert
318 *--------------------------------------------------------------------*/
319
320 extern FR_DOUBLE b[],xx[],h[],p[];
321
322 FR_DOUBLE xj[10], xji[10]; /* ist 3x3 +1 */
323
324 FR_DOUBLE rr2,ss2,r4,r3,s4,s3,rs4,dum,ept,emt;
325
326 FR_INT4 i,j,k,k3;
327
328 /*----------------------------------------------------------------------
329 * Klammern der Formfunktionen belegen
330 *---------------------------------------------------------------------*/
331 rr2= 2. * (*r) * (*r);
332 ss2= 2. * (*s) * (*s);
333 r4 = 4. * (*r);
334 r3 = 3. * (*r);
335 s4 = 4. * (*s);
336 s3 = 3. * (*s);
337 rs4= 4. * (*r) * (*s);
338 ept= 0.5* (1. + (*t));
339 emt= 0.5* (1. - (*t));
340
341 /*----------------------------------------------------------------------
342 * Formfunktionen obere und untere Ebene
343 *---------------------------------------------------------------------*/
344 h[ 1]= (rr2 + ss2 + rs4 - r3 - s3 + 1.) * ept;
345 h[ 2]= (rr2 - (*r)) * ept;
346 h[ 3]= (ss2 - (*s)) * ept;
347 h[ 4]= (r4 - 2*rr2 - rs4) * ept;
348 h[ 5]= rs4 * ept;
349 h[ 6]= (s4 - 2*ss2 - rs4) * ept;
350
351 h[ 7]= (rr2 + ss2 + rs4 - r3 - s3 + 1.) * emt;
352 h[ 8]= (rr2 - (*r)) * emt;
353 h[ 9]= (ss2 - (*s)) * emt;
354 h[10]= (r4 - 2*rr2 - rs4) * emt;
355 h[11]= rs4 * emt;
356 h[12]= (s4 - 2*ss2 - rs4) * emt;
357
358 /*----------------------------------------------------------------------
359 * Partielle Ableitung der Formfunktionen nach r
360 *---------------------------------------------------------------------*/
361 p[ 1]= (r4 + s4 - 3.) * ept;
362 p[ 2]= (r4 - 1.) * ept;
363 p[ 3]= 0.;
364 p[ 4]= (4. - 8*(*r) -s4) * ept;
365 p[ 5]= s4 * ept;
366 p[ 6]= -s4 * ept;
367
368 p[ 7]= (r4 + s4 - 3.) * emt;
369 p[ 8]= (r4 - 1.) * emt;
370 p[ 9]= 0.;
371 p[10]= (4. - 8*(*r) -s4) * emt;
372 p[11]= s4 * emt;
373 p[12]= -s4 * emt;
374
375 /*----------------------------------------------------------------------
376 * Partielle Ableitung der Formfunktionen nach s
377 *---------------------------------------------------------------------*/
378 p[13] = (s4 + r4 - 3.) * ept;
379 p[14] = 0.;
380 p[15] = (s4 - 1.) * ept;
381 p[16]= -r4 * ept;
382 p[17]= r4 * ept;
383 p[18]= (4. - r4 - 8*(*s)) * ept;
384
385 p[19] = (s4 + r4 - 3.) * emt;
386 p[20] = 0.;
387 p[21] = (s4 - 1.) * emt;
388 p[22]= -r4 * emt;
389 p[23]= r4 * emt;
390 p[24]= (4. - r4 - 8*(*s)) * emt;
391
392 /*----------------------------------------------------------------------
393 * Partielle Ableitung der Formfunktionen nach t
394 *---------------------------------------------------------------------*/
395 p[25]= (rr2 + ss2 + rs4 - r3 - s3 + 1.) * 0.5;
396 p[26]= (rr2 - (*r)) * 0.5;
397 p[27]= (ss2 - (*s)) * 0.5;
398 p[28]= (r4 - 2*rr2 - rs4) * 0.5;
399 p[29]= rs4 * 0.5;
400 p[30]= (s4 - 2*ss2 - rs4) * 0.5;
401
402 p[31]= -(rr2 + ss2 + rs4 - r3 - s3 + 1.) * 0.5;
403 p[32]= -(rr2 - (*r)) * 0.5;
404 p[33]= -(ss2 - (*s)) * 0.5;
405 p[34]= -(r4 - 2*rr2 - rs4) * 0.5;
406 p[35]= -rs4 * 0.5;
407 p[36]= -(s4 - 2*ss2 - rs4) * 0.5;
408
409 /*----------------------------------------------------------------------
410 * Jacobi-Matrix am Punkt (r,s,t) entwickeln
411 *---------------------------------------------------------------------*/
412 for(i = 1;i <= 3;i++)
413 {
414 for(j = 1;j <= 3;j++)
415 {
416 dum= 0.;
417 for(k = 1;k <= 12;k++)
418 {
419 dum+= p[(i-1)*12 + k] * xx[(j-1)*12 + k];
420 }
421 xj[(i-1)*3 + j]= dum;
422 }
423 }
424
425 /*----------------------------------------------------------------------
426 * Jacobi-Determinante am Punkt (r,s,t) entwickeln
427 *---------------------------------------------------------------------*/
428 (*det)= (xj[1] * xj[5] * xj[9]) - (xj[1] * xj[6] * xj[8]) +
429 (xj[2] * xj[6] * xj[7]) - (xj[2] * xj[4] * xj[9]) +
430 (xj[3] * xj[4] * xj[8]) - (xj[3] * xj[5] * xj[7]);
431
432 if(FR_FABS(*det) < 0.00000000001)
433 {
434 return(AL_JACNEG);
435 }
436
437 /*----------------------------------------------------------------------
438 * Berechnung der inversen Jacobi-Matrix
439 *---------------------------------------------------------------------*/
440 dum= 1./(*det);
441
442 xji[1]= (xj[5] * xj[9] - xj[8] * xj[6])*dum;
443 xji[2]= -(xj[2] * xj[9] - xj[8] * xj[3])*dum;
444 xji[3]= (xj[2] * xj[6] - xj[5] * xj[3])*dum;
445 xji[4]= -(xj[4] * xj[9] - xj[7] * xj[6])*dum;
446 xji[5]= (xj[1] * xj[9] - xj[7] * xj[3])*dum;
447 xji[6]= -(xj[1] * xj[6] - xj[4] * xj[3])*dum;
448 xji[7]= (xj[4] * xj[8] - xj[7] * xj[5])*dum;
449 xji[8]= -(xj[1] * xj[8] - xj[7] * xj[2])*dum;
450 xji[9]= (xj[1] * xj[5] - xj[4] * xj[2])*dum;
451
452 /*----------------------------------------------------------------------
453 * Entwickeln der Matrix b
454 *---------------------------------------------------------------------*/
455 for(i = 1;i <= 216;i++) /* 12 Knoten * 3 FG * 6 */
456 b[i]= 0.;
457
458 k3= 0;
459
460 for(k = 1;k <= 12;k++)
461 {
462 k3+= 3;
463
464 for(i = 1;i <= 3;i++)
465 {
466 b[ k3-2]+= xji[ i] * p[(i-1)*12 + k];
467 b[36 + k3-1]+= xji[3 + i] * p[(i-1)*12 + k];
468 b[72 + k3 ]+= xji[6 + i] * p[(i-1)*12 + k];
469 }
470 b[108 + k3-2]= b[36 + k3-1];
471 b[108 + k3-1]= b[ k3-2];
472
473 b[144 + k3-1]= b[72 + k3 ];
474 b[144 + k3 ]= b[36 + k3-1];
475
476 b[180 + k3-2]= b[72 + k3 ];
477 b[180 + k3 ]= b[ k3-2];
478 }
479
480 return(0);
481 }
482
483