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 * Z88 should compile and run under any Windows OS and UNIX OS and
22 * GTK+.
23 *
24 * This program is free software; you can redistribute it and/or modify
25 * it under the terms of the GNU General Public License as published by
26 * the Free Software Foundation; either version 2, or (at your option)
27 * any later version.
28 *
29 * This program is distributed in the hope that it will be useful,
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 * GNU General Public License for more details.
33 *
34 * You should have received a copy of the GNU General Public License
35 * along with this program; see the file COPYING. If not, write to
36 * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
37 ***********************************************************************/
38 /***********************************************************************
39 * diese Compilerunit umfasst: lqua88 - Elementsteifigkeitsroutine
40 * lb88 - Berechnung der Matrix b
41 * diese Compilerunit enthaelt Routinen, die gedanklich an FORTRAN-
42 * Quellen von H.J.Bathe, MIT, Cambridge, MA, USA angelehnt sind.
43 * 2.1.2010 Rieg
44 ***********************************************************************/
45
46 /***********************************************************************
47 * Fuer UNIX
48 ***********************************************************************/
49 #ifdef FR_UNIX
50 #include <z88r.h>
51 #endif
52
53 /***********************************************************************
54 * Fuer Windows
55 ***********************************************************************/
56 #ifdef FR_WIN
57 #include <z88r.h>
58 #endif
59
60 /***********************************************************************
61 * Fuer Windows & GTK+
62 ***********************************************************************/
63 #ifdef FR_GTKWIN
64 #include <z88r.h>
65 #endif
66
67 /***********************************************************************
68 * Functions
69 ***********************************************************************/
70 int lb88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s,FR_DOUBLE *t);
71
72 /***********************************************************************
73 * hier beginnt Function lqua88
74 ***********************************************************************/
lqua88(void)75 int lqua88(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= 24,i,lx,ly,lz,j,k,l;
92
93 int iret;
94
95 /*----------------------------------------------------------------------
96 * Gauss-Legendre Stuetzstellen
97 *---------------------------------------------------------------------*/
98 static FR_DOUBLE xg[17]= { 0.,
99 0., -.5773502691896, -.7745966692415, -.8611363115941,
100 0., +.5773502691896, 0., -.3399810435849,
101 0., 0., +.7745966692415, +.3399810435849,
102 0., 0., 0., +.8611363115941 };
103
104 /*----------------------------------------------------------------------
105 * Gauss-Legendre Integrationsgewichte
106 *---------------------------------------------------------------------*/
107 static FR_DOUBLE wgt[17]= { 0.,
108 2., 1., +.5555555555556, +.3478548451375,
109 0., 1., +.8888888888889, +.6521451548625,
110 0., 0., +.5555555555556, +.6521451548625,
111 0., 0., 0., +.3478548451375 };
112
113 /*----------------------------------------------------------------------
114 * xk und yk umspeichern
115 *---------------------------------------------------------------------*/
116 for(i = 1;i <= 8;i++)
117 {
118 xx[i] = xk[i];
119 xx[8 +i]= yk[i];
120 xx[16+i]= zk[i];
121 }
122
123 /*----------------------------------------------------------------------
124 * Materialkonstanten
125 *---------------------------------------------------------------------*/
126 f0= emode*(1.-rnuee) / ((1.+rnuee)*(1.-2.*rnuee));
127
128 f1= rnuee/(1.-rnuee) * f0;
129 f2= (1.-2.*rnuee) / (2.*(1.-rnuee)) * f0;
130
131 /*----------------------------------------------------------------------
132 * Elastizitaetsmatrix aufstellen
133 *---------------------------------------------------------------------*/
134 for(i = 1;i <= 36;i++)
135 d[i]= 0.;
136
137 d[1] = f0;
138 d[7] = f1;
139 d[13]= f1;
140 d[2] = f1;
141 d[8] = f0;
142 d[14]= f1;
143 d[3] = f1;
144 d[9] = f1;
145 d[15]= f0;
146 d[22]= f2;
147 d[29]= f2;
148 d[36]= f2;
149
150 /*----------------------------------------------------------------------
151 * Elementsteifigkeitsmatrix aufstellen
152 *---------------------------------------------------------------------*/
153 for(i = 1;i <= 576;i++)
154 se[i]= 0.;
155
156 for(lx = 1;lx <= intore;lx++) /* 120 */
157 {
158 r= xg[(lx-1)*4 + intore];
159 for(ly = 1;ly <= intore;ly++) /* 110 */
160 {
161 s= xg[(ly-1)*4 + intore];
162 for(lz = 1;lz <= intore;lz++) /* 110 */
163 {
164 t= xg[(lz-1)*4 + intore];
165
166 /*======================================================================
167 * Matrix b der partiellen Ableitungen & Jacobi Determinante holen
168 *=====================================================================*/
169 iret= lb88(&det,&r,&s,&t);
170 if(iret != 0) return(iret);
171
172 wt= wgt[(lx-1)*4 + intore] *
173 wgt[(ly-1)*4 + intore] *
174 wgt[(lz-1)*4 + intore] * det;
175
176 for(j = 1;j <= 24;j++) /* 90 */
177 {
178 for(k = 1;k <= 6;k++) /* 60 */
179 {
180 db[k]= 0.;
181 for(l = 1;l <= 6;l++) /* 50 */
182 {
183 db[k]= db[k] + d[(k-1)*6 + l] * b[(l-1)*24 + j];
184 } /* e 50 */
185 } /* e 60 */
186
187 for(i = j;i <= 24;i++) /* 80 */
188 {
189 stiff= 0.;
190 for(l = 1;l <= 6;l++) /* 70 */
191 {
192 stiff+= b[(l-1)*24 + i] * db[l];
193 } /* e 70 */
194 se[i+ne*(j-1)]= se[i+ne*(j-1)] + stiff * wt;
195 } /* e 80 */
196 } /* e 90 */
197 } /* e 100 */
198 } /* e 110 */
199 } /* e 120 */
200
201 for(j = 1;j <= 24;j++)
202 { /* 140 */
203 for(i = j;i <= 24;i++) /* 130 */
204 {
205 se[j+ne*(i-1)]= se[i+ne*(j-1)];
206 } /* e 130 */
207 } /* e 140 */
208
209 return(0);
210 }
211
212 /***********************************************************************
213 * hier beginnt Function lb88
214 ***********************************************************************/
lb88(FR_DOUBLE * det,FR_DOUBLE * r,FR_DOUBLE * s,FR_DOUBLE * t)215 int lb88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s,FR_DOUBLE *t)
216 {
217 /*---------------------------------------------------------------------
218 * xx geht rein, unveraendert (ex)
219 * b geht raus, neu (ex)
220 * det geht raus, neu
221 * r,s,t gehen rein, unveraendert
222 *--------------------------------------------------------------------*/
223
224 extern FR_DOUBLE b[],xx[],p[];
225
226 FR_DOUBLE xj[10], xji[10]; /* ist 3x3 +1 */
227
228 FR_DOUBLE dum,rs,rt,st;
229
230 FR_INT4 i,j,k,k3;
231
232 /*----------------------------------------------------------------------
233 * Faktoren der Ableitungen der Formfunktionen belegen
234 *---------------------------------------------------------------------*/
235 rs= (*r) * (*s) ;
236 rt= (*r) * (*t) ;
237 st= (*s) * (*t) ;
238
239 /*----------------------------------------------------------------------
240 * Partielle Ableitung der Formfunktionen nach r
241 *---------------------------------------------------------------------*/
242 p[1]= .125 *(+1. + (*s) + (*t) + st);
243 p[2]= .125 *(-1. - (*s) - (*t) - st);
244 p[3]= .125 *(-1. + (*s) - (*t) + st);
245 p[4]= .125 *(+1. - (*s) + (*t) - st);
246 p[5]= .125 *(+1. + (*s) - (*t) - st);
247 p[6]= .125 *(-1. - (*s) + (*t) + st);
248 p[7]= .125 *(-1. + (*s) + (*t) - st);
249 p[8]= .125 *(+1. - (*s) - (*t) + st);
250
251 /*----------------------------------------------------------------------
252 * Partielle Ableitung der Formfunktionen nach s
253 *---------------------------------------------------------------------*/
254 p[9] = .125 *(+1. + (*r) + (*t) + rt);
255 p[10]= .125 *(+1. - (*r) + (*t) - rt);
256 p[11]= .125 *(-1. + (*r) - (*t) + rt);
257 p[12]= .125 *(-1. - (*r) - (*t) - rt);
258 p[13]= .125 *(+1. + (*r) - (*t) - rt);
259 p[14]= .125 *(+1. - (*r) - (*t) + rt);
260 p[15]= .125 *(-1. + (*r) + (*t) - rt);
261 p[16]= .125 *(-1. - (*r) + (*t) + rt);
262
263 /*----------------------------------------------------------------------
264 * Partielle Ableitung der Formfunktionen nach t
265 *---------------------------------------------------------------------*/
266 p[17]= .125 *(+1. + (*r) + (*s) + rs);
267 p[18]= .125 *(+1. - (*r) + (*s) - rs);
268 p[19]= .125 *(+1. - (*r) - (*s) + rs);
269 p[20]= .125 *(+1. + (*r) - (*s) - rs);
270 p[21]= .125 *(-1. - (*r) - (*s) - rs);
271 p[22]= .125 *(-1. + (*r) - (*s) + rs);
272 p[23]= .125 *(-1. + (*r) + (*s) - rs);
273 p[24]= .125 *(-1. - (*r) + (*s) + rs);
274
275 /*----------------------------------------------------------------------
276 * Jacobi-Matrix am Punkt (r,s,t) entwickeln
277 *---------------------------------------------------------------------*/
278 for(i = 1;i <= 3;i++)
279 {
280 for(j = 1;j <= 3;j++)
281 {
282 dum= 0.;
283 for(k = 1;k <= 8;k++)
284 {
285 dum+= p[(i-1)*8 + k] * xx[(j-1)*8 + k];
286 }
287 xj[(i-1)*3 + j]= dum;
288 }
289 }
290
291 /*----------------------------------------------------------------------
292 * Jacobi-Determinante am Punkt (r,s,t) entwickeln
293 *---------------------------------------------------------------------*/
294 (*det)= (xj[1] * xj[5] * xj[9]) - (xj[1] * xj[6] * xj[8]) +
295 (xj[2] * xj[6] * xj[7]) - (xj[2] * xj[4] * xj[9]) +
296 (xj[3] * xj[4] * xj[8]) - (xj[3] * xj[5] * xj[7]);
297
298 if((*det) < 0.00000001)
299 return(AL_JACNEG);
300
301 /*----------------------------------------------------------------------
302 * Berechnung der inversen Jacobi-Matrix
303 *---------------------------------------------------------------------*/
304 dum= 1./(*det);
305
306 xji[1]= (xj[5] * xj[9] - xj[8] * xj[6])*dum;
307 xji[2]= -(xj[2] * xj[9] - xj[8] * xj[3])*dum;
308 xji[3]= (xj[2] * xj[6] - xj[5] * xj[3])*dum;
309 xji[4]= -(xj[4] * xj[9] - xj[7] * xj[6])*dum;
310 xji[5]= (xj[1] * xj[9] - xj[7] * xj[3])*dum;
311 xji[6]= -(xj[1] * xj[6] - xj[4] * xj[3])*dum;
312 xji[7]= (xj[4] * xj[8] - xj[7] * xj[5])*dum;
313 xji[8]= -(xj[1] * xj[8] - xj[7] * xj[2])*dum;
314 xji[9]= (xj[1] * xj[5] - xj[4] * xj[2])*dum;
315
316 /*----------------------------------------------------------------------
317 * Entwickeln der Matrix b
318 *---------------------------------------------------------------------*/
319 for(i = 1;i <= 144;i++)
320 b[i]= 0.;
321
322 k3= 0;
323
324 for(k = 1;k <= 8;k++)
325 {
326 k3+= 3;
327
328 for(i = 1;i <= 3;i++)
329 {
330 b[ k3-2]= b[ k3-2] + xji[ i] * p[(i-1)*8 + k];
331 b[24 + k3-1]= b[24 + k3-1] + xji[3 + i] * p[(i-1)*8 + k];
332 b[48 + k3 ]= b[48 + k3 ] + xji[6 + i] * p[(i-1)*8 + k];
333 }
334 b[72 + k3-2]= b[24 + k3-1];
335 b[72 + k3-1]= b[ k3-2];
336
337 b[96 + k3-1]= b[48+ k3 ];
338 b[96 + k3 ]= b[24 + k3-1];
339
340 b[120 + k3-2]= b[48 +k3 ];
341 b[120 + k3 ]= b[ k3-2];
342 }
343
344 return(0);
345 }
346
347